build_nat.c

Go to the documentation of this file.
00001 
00019 #include <string.h>
00020 #include <stdlib.h>
00021 #include <stdio.h>
00022 #include <grass/glocale.h>
00023 #include <grass/gis.h>
00024 #include <grass/Vect.h>
00025 
00037 int Vect_build_line_area(struct Map_info *Map, int iline, int side)
00038 {
00039     int j, area, isle, n_lines, line, type, direction;
00040     static int first = 1;
00041     long offset;
00042     struct Plus_head *plus;
00043     P_LINE *BLine;
00044     static struct line_pnts *Points, *APoints;
00045     plus_t *lines;
00046     double area_size;
00047 
00048     plus = &(Map->plus);
00049 
00050     G_debug(3, "Vect_build_line_area() line = %d, side = %d", iline, side);
00051 
00052     if (first) {
00053         Points = Vect_new_line_struct();
00054         APoints = Vect_new_line_struct();
00055         first = 0;
00056     }
00057 
00058     area = dig_line_get_area(plus, iline, side);
00059     if (area != 0) {
00060         G_debug(3, "  area/isle = %d -> skip", area);
00061         return 0;
00062     }
00063 
00064     n_lines = dig_build_area_with_line(plus, iline, side, &lines);
00065     G_debug(3, "  n_lines = %d", n_lines);
00066     if (n_lines < 1) {
00067         return 0;
00068     }                           /* area was not built */
00069 
00070     /* Area or island ? */
00071     Vect_reset_line(APoints);
00072     for (j = 0; j < n_lines; j++) {
00073         line = abs(lines[j]);
00074         BLine = plus->Line[line];
00075         offset = BLine->offset;
00076         G_debug(3, "  line[%d] = %d, offset = %ld", j, line, offset);
00077         type = Vect_read_line(Map, Points, NULL, line);
00078         if (lines[j] > 0)
00079             direction = GV_FORWARD;
00080         else
00081             direction = GV_BACKWARD;
00082         Vect_append_points(APoints, Points, direction);
00083     }
00084 
00085     dig_find_area_poly(APoints, &area_size);
00086     G_debug(3, "  area/isle size = %f", area_size);
00087 
00088     if (area_size > 0) {        /* area */
00089         /* add area structure to plus */
00090         area = dig_add_area(plus, n_lines, lines);
00091         if (area == -1) {       /* error */
00092             Vect_close(Map);
00093             G_fatal_error(_("Unable to add area (map closed, topo saved)"));
00094         }
00095         G_debug(3, "  -> area %d", area);
00096         return area;
00097     }
00098     else if (area_size < 0) {   /* island */
00099         isle = dig_add_isle(plus, n_lines, lines);
00100         if (isle == -1) {       /* error */
00101             Vect_close(Map);
00102             G_fatal_error(_("Unable to add isle (map closed, topo saved)"));
00103         }
00104         G_debug(3, "  -> isle %d", isle);
00105         return -isle;
00106     }
00107     else {
00108         /* TODO: What to do with such areas? Should be areas/isles of size 0 stored,
00109          *        so that may be found and cleaned by some utility
00110          *  Note: it would be useful for vertical closed polygons, but such would be added twice
00111          *        as area */
00112         G_warning(_("Area of size = 0.0 ignored"));
00113     }
00114     return 0;
00115 }
00116 
00126 int Vect_isle_find_area(struct Map_info *Map, int isle)
00127 {
00128     int j, line, node, sel_area, first, area, poly;
00129     static int first_call = 1;
00130     struct Plus_head *plus;
00131     P_LINE *Line;
00132     P_NODE *Node;
00133     P_ISLE *Isle;
00134     P_AREA *Area;
00135     double size, cur_size;
00136     BOUND_BOX box, abox;
00137     static struct ilist *List;
00138     static struct line_pnts *APoints;
00139 
00140     /* Note: We should check all isle points (at least) because if topology is not clean
00141      * and two areas overlap, isle which is not completely within area may be attached,
00142      * but it would take long time */
00143 
00144     G_debug(3, "Vect_isle_find_area () island = %d", isle);
00145     plus = &(Map->plus);
00146 
00147     if (plus->Isle[isle] == NULL) {
00148         G_warning(_("Request to find area outside nonexistent isle"));
00149         return 0;
00150     }
00151 
00152     if (first_call) {
00153         List = Vect_new_list();
00154         APoints = Vect_new_line_struct();
00155         first_call = 0;
00156     }
00157 
00158     Isle = plus->Isle[isle];
00159     line = abs(Isle->lines[0]);
00160     Line = plus->Line[line];
00161     node = Line->N1;
00162     Node = plus->Node[node];
00163 
00164     /* select areas by box */
00165     box.E = Node->x;
00166     box.W = Node->x;
00167     box.N = Node->y;
00168     box.S = Node->y;
00169     box.T = PORT_DOUBLE_MAX;
00170     box.B = -PORT_DOUBLE_MAX;
00171     Vect_select_areas_by_box(Map, &box, List);
00172     G_debug(3, "%d areas overlap island boundary point", List->n_values);
00173 
00174     sel_area = 0;
00175     cur_size = -1;
00176     first = 1;
00177     Vect_get_isle_box(Map, isle, &box);
00178     for (j = 0; j < List->n_values; j++) {
00179         area = List->value[j];
00180         G_debug(3, "area = %d", area);
00181 
00182         Area = plus->Area[area];
00183 
00184         /* Before other tests, simply exclude those areas inside isolated isles formed by one boundary */
00185         if (abs(Isle->lines[0]) == abs(Area->lines[0])) {
00186             G_debug(3, "  area inside isolated isle");
00187             continue;
00188         }
00189 
00190         /* Check box */
00191         /* Note: If build is run on large files of areas imported from nontopo format (shapefile)
00192          * attaching of isles takes very  long time because each area is also isle and select by
00193          * box all overlapping areas selects all areas with box overlapping first node. 
00194          * Then reading coordinates for all those areas would take a long time -> check first 
00195          * if isle's box is completely within area box */
00196         Vect_get_area_box(Map, area, &abox);
00197         if (box.E > abox.E || box.W < abox.W || box.N > abox.N ||
00198             box.S < abox.S) {
00199             G_debug(3, "  isle not completely inside area box");
00200             continue;
00201         }
00202 
00203         poly = Vect_point_in_area_outer_ring(Node->x, Node->y, Map, area);
00204         G_debug(3, "  poly = %d", poly);
00205 
00206         if (poly == 1) {        /* pint in area, but node is not part of area inside isle (would be poly == 2) */
00207             /* In rare case island is inside more areas in that case we have to calculate area
00208              * of outer ring and take the smaller */
00209             if (sel_area == 0) {        /* first */
00210                 sel_area = area;
00211             }
00212             else {              /* is not first */
00213                 if (cur_size < 0) {     /* second area */
00214                     /* This is slow, but should not be called often */
00215                     Vect_get_area_points(Map, sel_area, APoints);
00216                     G_begin_polygon_area_calculations();
00217                     cur_size =
00218                         G_area_of_polygon(APoints->x, APoints->y,
00219                                           APoints->n_points);
00220                     G_debug(3, "  first area size = %f (n points = %d)",
00221                             cur_size, APoints->n_points);
00222 
00223                 }
00224 
00225                 Vect_get_area_points(Map, area, APoints);
00226                 size =
00227                     G_area_of_polygon(APoints->x, APoints->y,
00228                                       APoints->n_points);
00229                 G_debug(3, "  area size = %f (n points = %d)", cur_size,
00230                         APoints->n_points);
00231 
00232                 if (size < cur_size) {
00233                     sel_area = area;
00234                     cur_size = size;
00235                 }
00236             }
00237             G_debug(3, "sel_area = %d cur_size = %f", sel_area, cur_size);
00238         }
00239     }
00240     if (sel_area > 0) {
00241         G_debug(3, "Island %d in area %d", isle, sel_area);
00242     }
00243     else {
00244         G_debug(3, "Island %d is not in area", isle);
00245     }
00246 
00247     return sel_area;
00248 }
00249 
00258 int Vect_attach_isle(struct Map_info *Map, int isle)
00259 {
00260     int sel_area;
00261     P_ISLE *Isle;
00262     struct Plus_head *plus;
00263 
00264     /* Note!: If topology is not clean and areas overlap, one island may fall to more areas
00265      *  (partially or fully). Before isle is attached to area it must be check if it is not attached yet */
00266     G_debug(3, "Vect_attach_isle (): isle = %d", isle);
00267 
00268     plus = &(Map->plus);
00269 
00270     sel_area = Vect_isle_find_area(Map, isle);
00271     G_debug(3, "      isle = %d -> area outside = %d", isle, sel_area);
00272     if (sel_area > 0) {
00273         Isle = plus->Isle[isle];
00274         if (Isle->area > 0) {
00275             G_debug(3,
00276                     "Attempt to attach isle %d to more areas (=>topology is not clean)",
00277                     isle);
00278         }
00279         else {
00280             Isle->area = sel_area;
00281             dig_area_add_isle(plus, sel_area, isle);
00282         }
00283     }
00284     return 0;
00285 }
00286 
00295 int Vect_attach_isles(struct Map_info *Map, BOUND_BOX * box)
00296 {
00297     int i, isle;
00298     static int first = 1;
00299     static struct ilist *List;
00300     struct Plus_head *plus;
00301 
00302     G_debug(3, "Vect_attach_isles ()");
00303 
00304     plus = &(Map->plus);
00305 
00306     if (first) {
00307         List = Vect_new_list();
00308         first = 0;
00309     }
00310 
00311     Vect_select_isles_by_box(Map, box, List);
00312     G_debug(3, "  number of isles to attach = %d", List->n_values);
00313 
00314     for (i = 0; i < List->n_values; i++) {
00315         isle = List->value[i];
00316         Vect_attach_isle(Map, isle);
00317     }
00318     return 0;
00319 }
00320 
00329 int Vect_attach_centroids(struct Map_info *Map, BOUND_BOX * box)
00330 {
00331     int i, sel_area, centr;
00332     static int first = 1;
00333     static struct ilist *List;
00334     P_AREA *Area;
00335     P_LINE *Line;
00336     struct Plus_head *plus;
00337 
00338     G_debug(3, "Vect_attach_centroids ()");
00339 
00340     plus = &(Map->plus);
00341 
00342     if (first) {
00343         List = Vect_new_list();
00344         first = 0;
00345     }
00346 
00347     /* Warning: If map is updated on level2, it may happen that previously correct island 
00348      * becomes incorrect. In that case, centroid of area forming the island is reattached 
00349      * to outer area, because island polygon is not excluded. 
00350      *
00351      * +-----------+     +-----------+
00352      * |   1       |     |   1       |
00353      * | +---+---+ |     | +---+---+ |     
00354      * | | 2 | 3 | |     | | 2 |     |   
00355      * | | x |   | |  -> | | x |     |  
00356      * | |   |   | |     | |   |     | 
00357      * | +---+---+ |     | +---+---+ |
00358      * |           |     |           |
00359      * +-----------+     +-----------+
00360      * centroid is       centroid is
00361      * attached to 2     reattached to 1
00362      *
00363      * Because of this, when the centroid is reattached to another area, it is always necessary
00364      * to check if original area exist, unregister centroid from previous area.
00365      * To simplify code, this is implemented so that centroid is always firs unregistered 
00366      * and if new area is found, it is registered again.
00367      *
00368      * This problem can be avoided altogether if properly attached centroids
00369      * are skipped
00370      * MM 2009
00371      */
00372 
00373     Vect_select_lines_by_box(Map, box, GV_CENTROID, List);
00374     G_debug(3, "  number of centroids to reattach = %d", List->n_values);
00375     for (i = 0; i < List->n_values; i++) {
00376         int orig_area;
00377 
00378         centr = List->value[i];
00379         Line = plus->Line[centr];
00380 
00381         /* only attach unregistered and duplicate centroids because 
00382          * 1) all properly attached centroids are properly attached, really! Don't touch.
00383          * 2) Vect_find_area() below does not always return the correct area
00384          * 3) it's faster
00385          */
00386         if (Line->left > 0)
00387             continue; 
00388 
00389         orig_area = Line->left;
00390 
00391         sel_area = Vect_find_area(Map, Line->E, Line->N);
00392         G_debug(3, "  centroid %d is in area %d", centr, sel_area);
00393         if (sel_area > 0) {
00394             Area = plus->Area[sel_area];
00395             if (Area->centroid == 0) {  /* first centroid */
00396                 G_debug(3, "  first centroid -> attach to area");
00397                 Area->centroid = centr;
00398                 Line->left = sel_area;
00399                 
00400                 if (sel_area != orig_area && plus->do_uplist)
00401                     dig_line_add_updated(plus, centr);
00402             }
00403             else if (Area->centroid != centr) { /* duplicate centroid */
00404                 /* Note: it cannot happen that Area->centroid == centr, because the centroid
00405                  * was not registered or a duplicate */
00406                 G_debug(3, "  duplicate centroid -> do not attach to area");
00407                 Line->left = -sel_area;
00408 
00409                 if (-sel_area != orig_area && plus->do_uplist)
00410                     dig_line_add_updated(plus, centr);
00411             }
00412         }
00413 
00414         if (sel_area != orig_area && plus->do_uplist)
00415             dig_line_add_updated(plus, centr);
00416     }
00417 
00418     return 0;
00419 }
00420 
00430 int Vect_build_nat(struct Map_info *Map, int build)
00431 {
00432     struct Plus_head *plus;
00433     int i, s, type, lineid;
00434     long offset;
00435     int side, line, area;
00436     struct line_pnts *Points, *APoints;
00437     struct line_cats *Cats;
00438     P_LINE *Line;
00439     P_NODE *Node;
00440     P_AREA *Area;
00441     BOUND_BOX box;
00442     struct ilist *List;
00443 
00444     G_debug(3, "Vect_build_nat() build = %d", build);
00445 
00446     plus = &(Map->plus);
00447 
00448     if (build == plus->built)
00449         return 1;               /* Do nothing */
00450 
00451     /* Check if upgrade or downgrade */
00452     if (build < plus->built) {  /* lower level request */
00453 
00454         /* release old sources (this also initializes structures and numbers of elements) */
00455         if (plus->built >= GV_BUILD_CENTROIDS && build < GV_BUILD_CENTROIDS) {
00456             /* reset info about areas stored for centroids */
00457             int nlines = Vect_get_num_lines(Map);
00458 
00459             for (line = 1; line <= nlines; line++) {
00460                 Line = plus->Line[line];
00461                 if (Line && Line->type == GV_CENTROID)
00462                     Line->left = 0;
00463             }
00464             dig_free_plus_areas(plus);
00465             dig_spidx_free_areas(plus);
00466             dig_free_plus_isles(plus);
00467             dig_spidx_free_isles(plus);
00468         }
00469 
00470 
00471         if (plus->built >= GV_BUILD_AREAS && build < GV_BUILD_AREAS) {
00472             /* reset info about areas stored for lines */
00473             int nlines = Vect_get_num_lines(Map);
00474 
00475             for (line = 1; line <= nlines; line++) {
00476                 Line = plus->Line[line];
00477                 if (Line && Line->type == GV_BOUNDARY) {
00478                     Line->left = 0;
00479                     Line->right = 0;
00480                 }
00481             }
00482             dig_free_plus_areas(plus);
00483             dig_spidx_free_areas(plus);
00484             dig_free_plus_isles(plus);
00485             dig_spidx_free_isles(plus);
00486         }
00487         if (plus->built >= GV_BUILD_BASE && build < GV_BUILD_BASE) {
00488             dig_free_plus_nodes(plus);
00489             dig_spidx_free_nodes(plus);
00490             dig_free_plus_lines(plus);
00491             dig_spidx_free_lines(plus);
00492         }
00493 
00494         plus->built = build;
00495         return 1;
00496     }
00497 
00498     Points = Vect_new_line_struct();
00499     APoints = Vect_new_line_struct();
00500     Cats = Vect_new_cats_struct();
00501     List = Vect_new_list();
00502 
00503     if (plus->built < GV_BUILD_BASE) {
00504         int npoints, format;
00505 
00506         format = G_info_format();
00507 
00508         /* 
00509          *  We shall go through all primitives in coor file and 
00510          *  add new node for each end point to nodes structure
00511          *  if the node with the same coordinates doesn't exist yet.
00512          */
00513 
00514         /* register lines, create nodes */
00515         Vect_rewind(Map);
00516         G_message(_("Registering primitives..."));
00517         i = 1;
00518         npoints = 0;
00519         while (1) {
00520             /* register line */
00521             type = Vect_read_next_line(Map, Points, Cats);
00522 
00523             /* Note: check for dead lines is not needed, because they are skipped by V1_read_next_line_nat() */
00524             if (type == -1) {
00525                 G_warning(_("Unable to read vector map"));
00526                 return 0;
00527             }
00528             else if (type == -2) {
00529                 break;
00530             }
00531 
00532             npoints += Points->n_points;
00533 
00534             offset = Map->head.last_offset;
00535 
00536             G_debug(3, "Register line: offset = %ld", offset);
00537             lineid = dig_add_line(plus, type, Points, offset);
00538             dig_line_box(Points, &box);
00539             if (lineid == 1)
00540                 Vect_box_copy(&(plus->box), &box);
00541             else
00542                 Vect_box_extend(&(plus->box), &box);
00543 
00544             /* Add all categories to category index */
00545             if (build == GV_BUILD_ALL) {
00546                 int c;
00547 
00548                 for (c = 0; c < Cats->n_cats; c++) {
00549                     dig_cidx_add_cat(plus, Cats->field[c], Cats->cat[c],
00550                                      lineid, type);
00551                 }
00552                 if (Cats->n_cats == 0)  /* add field 0, cat 0 */
00553                     dig_cidx_add_cat(plus, 0, 0, lineid, type);
00554             }
00555 
00556             if (G_verbose() > G_verbose_min() && i % 1000 == 0) {
00557                 if (format == G_INFO_FORMAT_PLAIN)
00558                     fprintf(stderr, "%d..", i);
00559                 else
00560                     fprintf(stderr, "%7d\b\b\b\b\b\b\b", i);
00561             }
00562             
00563             i++;
00564         }
00565         
00566         if ( (G_verbose() > G_verbose_min() ) && format != G_INFO_FORMAT_PLAIN )
00567             fprintf(stderr, "\r");
00568 
00569         G_message(_("%d primitives registered"), plus->n_lines);
00570         G_message(_("%d vertices registered"), npoints);
00571 
00572         plus->built = GV_BUILD_BASE;
00573     }
00574 
00575     if (build < GV_BUILD_AREAS)
00576         return 1;
00577 
00578     if (plus->built < GV_BUILD_AREAS) {
00579         /* Build areas */
00580         /* Go through all bundaries and try to build area for both sides */
00581         G_message(_("Building areas..."));
00582         for (i = 1; i <= plus->n_lines; i++) {
00583             G_percent(i, plus->n_lines, 1);
00584 
00585             /* build */
00586             if (plus->Line[i] == NULL) {
00587                 continue;
00588             }                   /* dead line */
00589             Line = plus->Line[i];
00590             if (Line->type != GV_BOUNDARY) {
00591                 continue;
00592             }
00593 
00594             for (s = 0; s < 2; s++) {
00595                 if (s == 0)
00596                     side = GV_LEFT;
00597                 else
00598                     side = GV_RIGHT;
00599 
00600                 G_debug(3, "Build area for line = %d, side = %d", i, side);
00601                 Vect_build_line_area(Map, i, side);
00602             }
00603         }
00604         G_message(_("%d areas built"), plus->n_areas);
00605         G_message(_("%d isles built"), plus->n_isles);
00606         plus->built = GV_BUILD_AREAS;
00607     }
00608 
00609     if (build < GV_BUILD_ATTACH_ISLES)
00610         return 1;
00611 
00612     /* Attach isles to areas */
00613     if (plus->built < GV_BUILD_ATTACH_ISLES) {
00614         G_message(_("Attaching islands..."));
00615         for (i = 1; i <= plus->n_isles; i++) {
00616             G_percent(i, plus->n_isles, 1);
00617             Vect_attach_isle(Map, i);
00618         }
00619         plus->built = GV_BUILD_ATTACH_ISLES;
00620     }
00621 
00622     if (build < GV_BUILD_CENTROIDS)
00623         return 1;
00624 
00625     /* Attach centroids to areas */
00626     if (plus->built < GV_BUILD_CENTROIDS) {
00627         int nlines;
00628 
00629         G_message(_("Attaching centroids..."));
00630 
00631         nlines = Vect_get_num_lines(Map);
00632         for (line = 1; line <= nlines; line++) {
00633             G_percent(line, nlines, 1);
00634 
00635             Line = plus->Line[line];
00636             if (!Line)
00637                 continue;       /* Dead */
00638 
00639             if (Line->type != GV_CENTROID)
00640                 continue;
00641 
00642             Node = plus->Node[Line->N1];
00643 
00644             area = Vect_find_area(Map, Node->x, Node->y);
00645 
00646             if (area > 0) {
00647                 G_debug(3, "Centroid (line=%d) in area %d", line, area);
00648 
00649                 Area = plus->Area[area];
00650 
00651                 if (Area->centroid == 0) {      /* first */
00652                     Area->centroid = line;
00653                     Line->left = area;
00654                 }
00655                 else {          /* duplicate */
00656                     Line->left = -area;
00657                 }
00658             }
00659         }
00660         plus->built = GV_BUILD_CENTROIDS;
00661     }
00662 
00663     /* Add areas to category index */
00664     for (area = 1; area <= plus->n_areas; area++) {
00665         int c;
00666 
00667         if (plus->Area[area] == NULL)
00668             continue;
00669 
00670         if (plus->Area[area]->centroid > 0) {
00671             Vect_read_line(Map, NULL, Cats, plus->Area[area]->centroid);
00672 
00673             for (c = 0; c < Cats->n_cats; c++) {
00674                 dig_cidx_add_cat(plus, Cats->field[c], Cats->cat[c], area,
00675                                  GV_AREA);
00676             }
00677         }
00678 
00679         if (plus->Area[area]->centroid == 0 || Cats->n_cats == 0)       /* no centroid or no cats */
00680             dig_cidx_add_cat(plus, 0, 0, area, GV_AREA);
00681     }
00682 
00683     return 1;
00684 }