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_NODE *Node;
00336     P_LINE *Line;
00337     struct Plus_head *plus;
00338 
00339     G_debug(3, "Vect_attach_centroids ()");
00340 
00341     plus = &(Map->plus);
00342 
00343     if (first) {
00344         List = Vect_new_list();
00345         first = 0;
00346     }
00347 
00348     /* Warning: If map is updated on level2, it may happen that previously correct island 
00349      * becomes incorrect. In that case, centroid of area forming the island is reattached 
00350      * to outer area, because island polygon is not excluded. 
00351      *
00352      * +-----------+     +-----------+
00353      * |   1       |     |   1       |
00354      * | +---+---+ |     | +---+---+ |     
00355      * | | 2 | 3 | |     | | 2 |     |   
00356      * | | x |   | |  -> | | x |     |  
00357      * | |   |   | |     | |   |     | 
00358      * | +---+---+ |     | +---+---+ |
00359      * |           |     |           |
00360      * +-----------+     +-----------+
00361      * centroid is       centroid is
00362      * attached to 2     reattached to 1
00363      *
00364      * Because of this, when the centroid is reattached to another area, it is always necessary
00365      * to check if original area exist, unregister centroid from previous area.
00366      * To simplify code, this is implemented so that centroid is always firs unregistered 
00367      * and if new area is found, it is registered again.
00368      */
00369 
00370     Vect_select_lines_by_box(Map, box, GV_CENTROID, List);
00371     G_debug(3, "  number of centroids to reattach = %d", List->n_values);
00372     for (i = 0; i < List->n_values; i++) {
00373         int orig_area;
00374 
00375         centr = List->value[i];
00376         Line = plus->Line[centr];
00377         Node = plus->Node[Line->N1];
00378 
00379         /* Unregister centroid */
00380         orig_area = Line->left;
00381         if (orig_area > 0) {
00382             if (plus->Area[orig_area] != NULL) {
00383                 plus->Area[orig_area]->centroid = 0;
00384             }
00385         }
00386         Line->left = 0;
00387 
00388         sel_area = Vect_find_area(Map, Node->x, Node->y);
00389         G_debug(3, "  centroid %d is in area %d", centr, sel_area);
00390         if (sel_area > 0) {
00391             Area = plus->Area[sel_area];
00392             if (Area->centroid == 0) {  /* first centroid */
00393                 G_debug(3, "  first centroid -> attach to area");
00394                 Area->centroid = centr;
00395                 Line->left = sel_area;
00396             }
00397             else if (Area->centroid != centr) { /* duplicate centroid */
00398                 /* Note: it cannot happen that Area->centroid == centr, because the centroid
00399                  * was previously unregistered */
00400                 G_debug(3, "  duplicate centroid -> do not attach to area");
00401                 Line->left = -sel_area;
00402             }
00403         }
00404 
00405         if (sel_area != orig_area && plus->do_uplist)
00406             dig_line_add_updated(plus, centr);
00407     }
00408 
00409     return 0;
00410 }
00411 
00421 int Vect_build_nat(struct Map_info *Map, int build)
00422 {
00423     struct Plus_head *plus;
00424     int i, s, type, lineid;
00425     long offset;
00426     int side, line, area;
00427     struct line_pnts *Points, *APoints;
00428     struct line_cats *Cats;
00429     P_LINE *Line;
00430     P_NODE *Node;
00431     P_AREA *Area;
00432     BOUND_BOX box;
00433     struct ilist *List;
00434 
00435     G_debug(3, "Vect_build_nat() build = %d", build);
00436 
00437     plus = &(Map->plus);
00438 
00439     if (build == plus->built)
00440         return 1;               /* Do nothing */
00441 
00442     /* Check if upgrade or downgrade */
00443     if (build < plus->built) {  /* lower level request */
00444 
00445         /* release old sources (this also initializes structures and numbers of elements) */
00446         if (plus->built >= GV_BUILD_CENTROIDS && build < GV_BUILD_CENTROIDS) {
00447             /* reset info about areas stored for centroids */
00448             int nlines = Vect_get_num_lines(Map);
00449 
00450             for (line = 1; line <= nlines; line++) {
00451                 Line = plus->Line[line];
00452                 if (Line && Line->type == GV_CENTROID)
00453                     Line->left = 0;
00454             }
00455             dig_free_plus_areas(plus);
00456             dig_spidx_free_areas(plus);
00457             dig_free_plus_isles(plus);
00458             dig_spidx_free_isles(plus);
00459         }
00460 
00461 
00462         if (plus->built >= GV_BUILD_AREAS && build < GV_BUILD_AREAS) {
00463             /* reset info about areas stored for lines */
00464             int nlines = Vect_get_num_lines(Map);
00465 
00466             for (line = 1; line <= nlines; line++) {
00467                 Line = plus->Line[line];
00468                 if (Line && Line->type == GV_BOUNDARY) {
00469                     Line->left = 0;
00470                     Line->right = 0;
00471                 }
00472             }
00473             dig_free_plus_areas(plus);
00474             dig_spidx_free_areas(plus);
00475             dig_free_plus_isles(plus);
00476             dig_spidx_free_isles(plus);
00477         }
00478         if (plus->built >= GV_BUILD_BASE && build < GV_BUILD_BASE) {
00479             dig_free_plus_nodes(plus);
00480             dig_spidx_free_nodes(plus);
00481             dig_free_plus_lines(plus);
00482             dig_spidx_free_lines(plus);
00483         }
00484 
00485         plus->built = build;
00486         return 1;
00487     }
00488 
00489     Points = Vect_new_line_struct();
00490     APoints = Vect_new_line_struct();
00491     Cats = Vect_new_cats_struct();
00492     List = Vect_new_list();
00493 
00494     if (plus->built < GV_BUILD_BASE) {
00495         int npoints, format;
00496 
00497         format = G_info_format();
00498 
00499         /* 
00500          *  We shall go through all primitives in coor file and 
00501          *  add new node for each end point to nodes structure
00502          *  if the node with the same coordinates doesn't exist yet.
00503          */
00504 
00505         /* register lines, create nodes */
00506         Vect_rewind(Map);
00507         G_message(_("Registering primitives..."));
00508         i = 1;
00509         npoints = 0;
00510         while (1) {
00511             /* register line */
00512             type = Vect_read_next_line(Map, Points, Cats);
00513 
00514             /* Note: check for dead lines is not needed, because they are skipped by V1_read_next_line_nat() */
00515             if (type == -1) {
00516                 G_warning(_("Unable to read vector map"));
00517                 return 0;
00518             }
00519             else if (type == -2) {
00520                 break;
00521             }
00522 
00523             npoints += Points->n_points;
00524 
00525             offset = Map->head.last_offset;
00526 
00527             G_debug(3, "Register line: offset = %ld", offset);
00528             lineid = dig_add_line(plus, type, Points, offset);
00529             dig_line_box(Points, &box);
00530             if (lineid == 1)
00531                 Vect_box_copy(&(plus->box), &box);
00532             else
00533                 Vect_box_extend(&(plus->box), &box);
00534 
00535             /* Add all categories to category index */
00536             if (build == GV_BUILD_ALL) {
00537                 int c;
00538 
00539                 for (c = 0; c < Cats->n_cats; c++) {
00540                     dig_cidx_add_cat(plus, Cats->field[c], Cats->cat[c],
00541                                      lineid, type);
00542                 }
00543                 if (Cats->n_cats == 0)  /* add field 0, cat 0 */
00544                     dig_cidx_add_cat(plus, 0, 0, lineid, type);
00545             }
00546 
00547             if (G_verbose() > G_verbose_min() && i % 1000 == 0) {
00548                 if (format == G_INFO_FORMAT_PLAIN)
00549                     fprintf(stderr, "%d..", i);
00550                 else
00551                     fprintf(stderr, "%7d\b\b\b\b\b\b\b", i);
00552             }
00553             
00554             i++;
00555         }
00556         
00557         if ( (G_verbose() > G_verbose_min() ) && format != G_INFO_FORMAT_PLAIN )
00558             fprintf(stderr, "\r");
00559 
00560         G_message(_("%d primitives registered"), plus->n_lines);
00561         G_message(_("%d vertices registered"), npoints);
00562 
00563         plus->built = GV_BUILD_BASE;
00564     }
00565 
00566     if (build < GV_BUILD_AREAS)
00567         return 1;
00568 
00569     if (plus->built < GV_BUILD_AREAS) {
00570         /* Build areas */
00571         /* Go through all bundaries and try to build area for both sides */
00572         G_message(_("Building areas..."));
00573         for (i = 1; i <= plus->n_lines; i++) {
00574             G_percent(i, plus->n_lines, 1);
00575 
00576             /* build */
00577             if (plus->Line[i] == NULL) {
00578                 continue;
00579             }                   /* dead line */
00580             Line = plus->Line[i];
00581             if (Line->type != GV_BOUNDARY) {
00582                 continue;
00583             }
00584 
00585             for (s = 0; s < 2; s++) {
00586                 if (s == 0)
00587                     side = GV_LEFT;
00588                 else
00589                     side = GV_RIGHT;
00590 
00591                 G_debug(3, "Build area for line = %d, side = %d", i, side);
00592                 Vect_build_line_area(Map, i, side);
00593             }
00594         }
00595         G_message(_("%d areas built"), plus->n_areas);
00596         G_message(_("%d isles built"), plus->n_isles);
00597         plus->built = GV_BUILD_AREAS;
00598     }
00599 
00600     if (build < GV_BUILD_ATTACH_ISLES)
00601         return 1;
00602 
00603     /* Attach isles to areas */
00604     if (plus->built < GV_BUILD_ATTACH_ISLES) {
00605         G_message(_("Attaching islands..."));
00606         for (i = 1; i <= plus->n_isles; i++) {
00607             G_percent(i, plus->n_isles, 1);
00608             Vect_attach_isle(Map, i);
00609         }
00610         plus->built = GV_BUILD_ATTACH_ISLES;
00611     }
00612 
00613     if (build < GV_BUILD_CENTROIDS)
00614         return 1;
00615 
00616     /* Attach centroids to areas */
00617     if (plus->built < GV_BUILD_CENTROIDS) {
00618         int nlines;
00619 
00620         G_message(_("Attaching centroids..."));
00621 
00622         nlines = Vect_get_num_lines(Map);
00623         for (line = 1; line <= nlines; line++) {
00624             G_percent(line, nlines, 1);
00625 
00626             Line = plus->Line[line];
00627             if (!Line)
00628                 continue;       /* Dead */
00629 
00630             if (Line->type != GV_CENTROID)
00631                 continue;
00632 
00633             Node = plus->Node[Line->N1];
00634 
00635             area = Vect_find_area(Map, Node->x, Node->y);
00636 
00637             if (area > 0) {
00638                 G_debug(3, "Centroid (line=%d) in area %d", line, area);
00639 
00640                 Area = plus->Area[area];
00641 
00642                 if (Area->centroid == 0) {      /* first */
00643                     Area->centroid = line;
00644                     Line->left = area;
00645                 }
00646                 else {          /* duplicate */
00647                     Line->left = -area;
00648                 }
00649             }
00650         }
00651         plus->built = GV_BUILD_CENTROIDS;
00652     }
00653 
00654     /* Add areas to category index */
00655     for (area = 1; area <= plus->n_areas; area++) {
00656         int c;
00657 
00658         if (plus->Area[area] == NULL)
00659             continue;
00660 
00661         if (plus->Area[area]->centroid > 0) {
00662             Vect_read_line(Map, NULL, Cats, plus->Area[area]->centroid);
00663 
00664             for (c = 0; c < Cats->n_cats; c++) {
00665                 dig_cidx_add_cat(plus, Cats->field[c], Cats->cat[c], area,
00666                                  GV_AREA);
00667             }
00668         }
00669 
00670         if (plus->Area[area]->centroid == 0 || Cats->n_cats == 0)       /* no centroid or no cats */
00671             dig_cidx_add_cat(plus, 0, 0, area, GV_AREA);
00672     }
00673 
00674     return 1;
00675 }

Generated on Thu Jul 16 13:21:17 2009 for GRASS Programmer's Manual by  doxygen 1.5.6