build_nat.c

Go to the documentation of this file.
00001 /*
00002 ****************************************************************************
00003 *
00004 * MODULE:       Vector library 
00005 *               
00006 * AUTHOR(S):    Original author CERL, probably Dave Gerdes or Mike Higgins.
00007 *               Update to GRASS 5.7 Radim Blazek and David D. Gray.
00008 *
00009 * PURPOSE:      Higher level functions for reading/writing/manipulating vectors.
00010 *
00011 * COPYRIGHT:    (C) 2001 by the GRASS Development Team
00012 *
00013 *               This program is free software under the GNU General Public
00014 *               License (>=v2). Read the file COPYING that comes with GRASS
00015 *               for details.
00016 *
00017 *****************************************************************************/
00018 #include <string.h>
00019 #include <stdlib.h>
00020 #include <stdio.h>
00021 #include <grass/glocale.h>
00022 #include <grass/gis.h>
00023 #include <grass/Vect.h>
00024 
00025 extern FILE *Msgout;
00026 extern int prnmsg ( char *msg, ...) ;
00027 
00034 int
00035 Vect_build_line_area ( struct Map_info *Map, int iline, int side )
00036 {
00037     int    j, area, isle, n_lines, line, type, direction;
00038     static int first = 1;
00039     long   offset;
00040     struct Plus_head *plus ;
00041     P_LINE *BLine;
00042     static struct line_pnts *Points, *APoints;
00043     plus_t *lines;
00044     double area_size;
00045 
00046     plus = &(Map->plus);
00047 
00048     G_debug ( 3, "Vect_build_line_area() line = %d, side = %d", iline, side );
00049     
00050     if ( first ) {
00051         Points = Vect_new_line_struct ();
00052         APoints = Vect_new_line_struct ();
00053         first = 0;
00054     }
00055 
00056     area = dig_line_get_area (plus, iline, side); 
00057     if ( area != 0 ) {
00058         G_debug ( 3, "  area/isle = %d -> skip", area );
00059         return 0;
00060     }
00061     
00062     n_lines = dig_build_area_with_line (plus, iline, side, &lines);
00063     G_debug ( 3, "  n_lines = %d", n_lines );
00064     if ( n_lines < 1 ) { return 0; } /* area was not built */
00065 
00066     /* Area or island ? */
00067     Vect_reset_line ( APoints ); 
00068     for (j = 0; j < n_lines; j++){
00069         line = abs(lines[j]);
00070         BLine = plus->Line[line];
00071         offset = BLine->offset;
00072         G_debug ( 3, "  line[%d] = %d, offset = %d", j, line, offset );
00073         type = Vect_read_line (Map, Points, NULL, line );
00074         if ( lines[j] > 0 ) direction = GV_FORWARD; 
00075         else direction = GV_BACKWARD;
00076         Vect_append_points ( APoints, Points, direction);
00077     }
00078     
00079     dig_find_area_poly (APoints, &area_size);
00080     G_debug ( 3, "  area/isle size = %f", area_size );
00081     
00082     if (area_size > 0) {  /* area */
00083         /* add area structure to plus */
00084         area = dig_add_area (plus, n_lines, lines);
00085         if ( area == -1 ) { /* error */
00086             Vect_close ( Map );
00087             G_fatal_error ( _("Cannot add area (map closed, topo saved)") );
00088         }
00089         G_debug ( 3, "  -> area %d", area );
00090         return area;
00091     } else if (area_size < 0) { /* island */
00092         isle = dig_add_isle (plus, n_lines, lines);
00093         if ( isle == -1 ) { /* error */
00094             Vect_close ( Map );
00095             G_fatal_error ( _("Cannot add isle (map closed, topo saved)") );
00096         }
00097         G_debug ( 3, "  -> isle %d", isle );
00098         return -isle;
00099     } else { 
00100         /* TODO: What to do with such areas? Should be areas/isles of size 0 stored,
00101         *        so that may be found and cleaned by some utility
00102         *  Note: it would be useful for vertical closed polygons, but such would be added twice
00103         *        as area */
00104         G_warning (_("Area of size = 0.0 ignored")); 
00105     }
00106     return 0;
00107 }
00108 
00115 int
00116 Vect_isle_find_area ( struct Map_info *Map, int isle ) 
00117 {
00118     int    j, line, node, sel_area, first, area, poly;
00119     static int first_call = 1;
00120     struct Plus_head *plus ;
00121     P_LINE *Line;
00122     P_NODE *Node;
00123     P_ISLE *Isle;
00124     P_AREA *Area;
00125     double  size, cur_size;
00126     BOUND_BOX box, abox;
00127     static struct ilist *List;
00128     static struct line_pnts *APoints;
00129 
00130     /* Note: We should check all isle points (at least) because if topology is not clean
00131      * and two areas overlap, isle which is not completely within area may be attached,
00132      * but it would take long time */
00133 
00134     G_debug ( 3, "Vect_isle_find_area () island = %d", isle );
00135     plus = &(Map->plus);
00136 
00137     if (  plus->Isle[isle] == NULL ) {
00138         G_warning ("Request to find area outside nonexisting isle");
00139         return 0;
00140     }
00141 
00142     if ( first_call ) {
00143         List = Vect_new_list ();
00144         APoints = Vect_new_line_struct ();
00145         first_call = 0;
00146     }
00147 
00148     Isle = plus->Isle[isle];
00149     line = abs(Isle->lines[0]);
00150     Line = plus->Line[line];
00151     node = Line->N1;
00152     Node = plus->Node[node];
00153     
00154     /* select areas by box */
00155     box.E = Node->x; box.W = Node->x; box.N = Node->y; box.S = Node->y;
00156     box.T = PORT_DOUBLE_MAX; box.B = -PORT_DOUBLE_MAX;
00157     Vect_select_areas_by_box (Map, &box, List);
00158     G_debug ( 3, "%d areas overlap island boundary point", List->n_values );
00159    
00160     sel_area = 0; cur_size = -1;
00161     first = 1;
00162     Vect_get_isle_box ( Map, isle, &box);
00163     for (j = 0; j < List->n_values; j++) {
00164         area = List->value[j];
00165         G_debug ( 3, "area = %d", area );
00166 
00167         Area = plus->Area[area];
00168         
00169         /* Before other tests, simply exclude those areas inside isolated isles formed by one boundary */
00170         if ( abs(Isle->lines[0]) == abs(Area->lines[0]) ) {
00171             G_debug ( 3, "  area inside isolated isle" );
00172             continue;
00173         }
00174         
00175         /* Check box */
00176         /* Note: If build is run on large files of areas imported from nontopo format (shapefile)
00177          * attaching of isles takes very  long time because each area is also isle and select by
00178          * box all overlaping areas selects all areas with box overlaping first node. 
00179          * Then reading coordinates for all those areas would take a long time -> check first 
00180          * if isle's box is completely within area box */
00181         Vect_get_area_box ( Map, area, &abox);
00182         if ( box.E > abox.E || box.W < abox.W || box.N > abox.N || box.S < abox.S ) {
00183             G_debug ( 3, "  isle not completely inside area box" );
00184             continue;
00185         }
00186         
00187         poly = Vect_point_in_area_outer_ring ( Node->x, Node->y, Map, area);
00188         G_debug ( 3, "  poly = %d", poly );
00189         
00190         if ( poly == 1) { /* pint in area, but node is not part of area inside isle (would be poly == 2) */
00191             /* In rare case island is inside more areas in that case we have to calculate area
00192              * of outer ring and take the smaller */
00193             if ( sel_area == 0 ) { /* first */
00194                 sel_area = area;
00195             } else { /* is not first */
00196                 if ( cur_size < 0 ) { /* second area */
00197                     /* This is slow, but should not be called often */
00198                     Vect_get_area_points (Map, sel_area, APoints);
00199                     G_begin_polygon_area_calculations();
00200                     cur_size = G_area_of_polygon(APoints->x, APoints->y, APoints->n_points);
00201                     G_debug ( 3, "  first area size = %f (n points = %d)", cur_size, APoints->n_points );
00202 
00203                 } 
00204 
00205                 Vect_get_area_points (Map, area, APoints);
00206                 size = G_area_of_polygon(APoints->x, APoints->y, APoints->n_points);
00207                 G_debug ( 3, "  area size = %f (n points = %d)", cur_size, APoints->n_points );
00208 
00209                 if ( size < cur_size ) {
00210                     sel_area = area;
00211                     cur_size = size;
00212                 }
00213             }
00214             G_debug ( 3, "sel_area = %d cur_size = %f", sel_area, cur_size );
00215         }
00216     }
00217     if ( sel_area > 0 ) {       
00218         G_debug ( 3, "Island %d in area %d", isle, sel_area );
00219     } else { 
00220         G_debug ( 3, "Island %d is not in area", isle );
00221     }
00222 
00223     return sel_area;
00224 }
00225 
00232 int
00233 Vect_attach_isle ( struct Map_info *Map, int isle ) 
00234 {
00235     int sel_area;
00236     P_ISLE *Isle;
00237     struct Plus_head *plus ;
00238 
00239     /* Note!: If topology is not clean and areas overlap, one island may fall to more areas
00240     *  (partially or fully). Before isle is attached to area it must be check if it is not attached yet */
00241     G_debug ( 3, "Vect_attach_isle (): isle = %d", isle);
00242 
00243     plus = &(Map->plus);
00244     
00245     sel_area = Vect_isle_find_area ( Map, isle );
00246     G_debug ( 3, "      isle = %d -> area outside = %d", isle, sel_area);
00247     if ( sel_area > 0 ) {
00248         Isle = plus->Isle[isle];
00249         if ( Isle->area > 0 ) {
00250             G_debug (3, "Attempt to attach isle %d to more areas (=>topology is not clean)");
00251         } else {
00252             Isle->area = sel_area;
00253             dig_area_add_isle ( plus, sel_area, isle );
00254         }
00255     }
00256     return 0;
00257 }
00258 
00265 int
00266 Vect_attach_isles ( struct Map_info *Map, BOUND_BOX *box )
00267 {
00268     int i, isle;
00269     static int first = 1;
00270     static struct ilist *List;
00271     struct Plus_head *plus ;
00272 
00273     G_debug ( 3, "Vect_attach_isles ()");
00274 
00275     plus = &(Map->plus);
00276 
00277     if ( first ) {
00278         List = Vect_new_list (); 
00279         first = 0;
00280     }
00281     
00282     Vect_select_isles_by_box ( Map, box, List);
00283     G_debug ( 3, "  number of isles to attach = %d", List->n_values);
00284 
00285     for (i = 0; i < List->n_values; i++) { 
00286         isle = List->value[i];
00287         Vect_attach_isle ( Map, isle );
00288     }
00289     return 0;
00290 }
00291 
00298 int
00299 Vect_attach_centroids ( struct Map_info *Map, BOUND_BOX *box ) 
00300 {
00301     int i, sel_area, centr;
00302     static int first = 1;
00303     static struct ilist *List;
00304     P_AREA *Area;
00305     P_NODE *Node;
00306     P_LINE *Line;
00307     struct Plus_head *plus ;
00308 
00309     G_debug ( 3, "Vect_attach_centroids ()");
00310     
00311     plus = &(Map->plus);
00312 
00313     if ( first ) {
00314         List = Vect_new_list (); 
00315         first = 0;
00316     }
00317 
00318     /* Warning: If map is updated on level2, it may happen that previously correct island 
00319      * becomes incorrect. In that case, centroid of area forming the island is reattached 
00320      * to outer area, because island polygon is not excluded. 
00321      *
00322      * +-----------+     +-----------+
00323      * |   1       |     |   1       |
00324      * | +---+---+ |     | +---+---+ |     
00325      * | | 2 | 3 | |     | | 2 |     |   
00326      * | | x |   | |  -> | | x |     |  
00327      * | |   |   | |     | |   |     | 
00328      * | +---+---+ |     | +---+---+ |
00329      * |           |     |           |
00330      * +-----------+     +-----------+
00331      * centroid is       centroid is
00332      * attached to 2     reattached to 1
00333      *
00334      * Because of this, when the centroid is reattached to another area, it is always necessary
00335      * to check if original area exist, unregister centroid from previous area.
00336      * To simplify code, this is implemented so that centroid is always firs unregistered 
00337      * and if new area is found, it is registered again.
00338      */
00339     
00340     Vect_select_lines_by_box ( Map, box, GV_CENTROID, List);
00341     G_debug ( 3, "  number of centroids to reattach = %d", List->n_values);
00342     for (i = 0; i < List->n_values; i++) { 
00343         int orig_area;
00344         centr = List->value[i];
00345         Line = plus->Line[centr];
00346         Node = plus->Node[Line->N1];
00347 
00348         /* Unregister centroid */
00349         orig_area = Line->left;
00350         if ( orig_area > 0 ) {
00351             if ( plus->Area[orig_area] != NULL ) {
00352                 plus->Area[orig_area]->centroid = 0;
00353             }
00354         }
00355         Line->left = 0;
00356         
00357         sel_area = Vect_find_area ( Map, Node->x, Node->y );
00358         G_debug ( 3, "  centroid %d is in area %d", centr, sel_area);
00359         if ( sel_area > 0 ) {
00360             Area = plus->Area[sel_area];
00361             if ( Area->centroid == 0 ) { /* first centroid */
00362                 G_debug ( 3, "  first centroid -> attach to area");
00363                 Area->centroid = centr;
00364                 Line->left = sel_area;
00365             } else if ( Area->centroid != centr ) {  /* duplicate centroid */
00366                 /* Note: it cannot happen that Area->centroid == centr, because the centroid
00367                  * was previously unregistered */
00368                 G_debug ( 3, "  duplicate centroid -> do not attach to area");
00369                 Line->left = -sel_area;
00370             }
00371         }
00372             
00373         if ( sel_area != orig_area && plus->do_uplist ) dig_line_add_updated ( plus, centr );
00374     }
00375         
00376     return 0;
00377 }
00378     
00385 int
00386 Vect_build_nat ( struct Map_info *Map, int build, FILE *msgout )
00387 {
00388     struct Plus_head *plus ;
00389     int    i, j, s, type, lineid;
00390     long offset; 
00391     int    side, line, area;
00392     struct line_pnts *Points, *APoints;
00393     struct line_cats *Cats;
00394     P_LINE *Line;
00395     P_NODE *Node;
00396     P_AREA *Area;
00397     BOUND_BOX box;
00398     struct ilist *List;
00399 
00400     G_debug (3, "Vect_build_nat() build = %d", build);
00401     
00402     plus = &(Map->plus);
00403     Msgout = msgout;
00404     
00405     if ( build == plus->built ) return 1; /* Do nothing */
00406     
00407     /* Check if upgrade or downgrade */
00408     if ( build < plus->built ) { /* lower level request */
00409 
00410         /* release old sources (this also initializes structures and numbers of elements) */
00411         if ( plus->built >= GV_BUILD_CENTROIDS && build < GV_BUILD_CENTROIDS) {
00412             /* reset info about areas stored for centroids */
00413             int nlines = Vect_get_num_lines (Map);
00414             for ( line = 1; line <= nlines; line++ ) {
00415                 Line = plus->Line[line];
00416                 if ( Line && Line->type == GV_CENTROID ) 
00417                     Line->left = 0;
00418             }
00419             dig_free_plus_areas (plus);
00420             dig_spidx_free_areas (plus);
00421             dig_free_plus_isles (plus);
00422             dig_spidx_free_isles (plus);
00423         }
00424 
00425 
00426         if ( plus->built >= GV_BUILD_AREAS && build < GV_BUILD_AREAS) {
00427             /* reset info about areas stored for lines */
00428             int nlines = Vect_get_num_lines (Map);
00429             for ( line = 1; line <= nlines; line++ ) {
00430                 Line = plus->Line[line];
00431                 if ( Line && Line->type == GV_BOUNDARY ) {
00432                     Line->left = 0;
00433                     Line->right = 0;
00434                 }
00435             }
00436             dig_free_plus_areas (plus);
00437             dig_spidx_free_areas (plus);
00438             dig_free_plus_isles (plus);
00439             dig_spidx_free_isles (plus);
00440         }
00441         if ( plus->built >= GV_BUILD_BASE && build < GV_BUILD_BASE) {
00442             dig_free_plus_nodes (plus);
00443             dig_spidx_free_nodes (plus);
00444             dig_free_plus_lines (plus);
00445             dig_spidx_free_lines (plus);
00446         }
00447         
00448         plus->built = build;
00449         return 1;
00450     }
00451     
00452     Points = Vect_new_line_struct ();
00453     APoints = Vect_new_line_struct ();
00454     Cats = Vect_new_cats_struct ();
00455     List = Vect_new_list ();
00456 
00457     if ( plus->built < GV_BUILD_BASE ) { 
00458         /* 
00459         *  We shall go through all primitives in coor file and 
00460         *  add new node for each end point to nodes structure
00461         *  if the node with the same coordinates doesn't exist yet.
00462         */
00463        
00464         /* register lines, create nodes */ 
00465         Vect_rewind ( Map );
00466         prnmsg (_("Registering lines: "));
00467         i = 1; j = 1;
00468         while ( 1 ) {
00469             /* register line */
00470             type = Vect_read_next_line (Map, Points, Cats);
00471             /* Note: check for dead lines is not needed, because they are skipped by V1_read_next_line_nat() */
00472             if ( type == -1 ) { 
00473                 fprintf (stderr, "\nERROR: vector file - can't read\n" );
00474                 return 0;
00475             } else if ( type == -2 ) {
00476                 break;
00477             }
00478             
00479             offset = Map->head.last_offset;
00480 
00481             G_debug ( 3, "Register line: offset = %ld", offset );
00482             lineid = dig_add_line ( plus, type, Points, offset );
00483             dig_line_box ( Points, &box );
00484             if ( lineid == 1 ) 
00485                 Vect_box_copy (&(plus->box), &box);
00486             else
00487                 Vect_box_extend (&(plus->box), &box);
00488 
00489             /* Add all categories to category index */
00490             if ( build == GV_BUILD_ALL ) {
00491                 int c; 
00492 
00493                 for ( c = 0; c < Cats->n_cats; c++ ) {
00494                     dig_cidx_add_cat (plus, Cats->field[c], Cats->cat[c], lineid, type);
00495                 }
00496                 if ( Cats->n_cats == 0 ) /* add field 0, cat 0 */
00497                     dig_cidx_add_cat (plus, 0, 0, lineid, type);
00498             }
00499             
00500             /* print progress */
00501             if ( i == 1000 ) {
00502                 prnmsg ("%7d\b\b\b\b\b\b\b", j);
00503                 i = 0;
00504             }
00505             i++; j++;
00506         }
00507         prnmsg (_("\r%d primitives registered      \n"), plus->n_lines);
00508 
00509         plus->built = GV_BUILD_BASE;
00510     }
00511 
00512     if ( build < GV_BUILD_AREAS ) return 1;
00513 
00514     if ( plus->built < GV_BUILD_AREAS ) {
00515         /* Build areas */
00516         /* Go through all bundaries and try to build area for both sides */
00517         prnmsg ("Building areas: ");
00518         for (i = 1; i <= plus->n_lines; i++) {
00519             G_percent2 ( i, plus->n_lines, 1, msgout );
00520 
00521             /* build */
00522             if ( plus->Line[i] == NULL ) { continue; } /* dead line */
00523             Line = plus->Line[i];
00524             if ( Line->type != GV_BOUNDARY ) { continue; }
00525 
00526             for (s = 0; s < 2; s++) {
00527                 if ( s == 0 ) side = GV_LEFT;
00528                 else side = GV_RIGHT;
00529                
00530                 G_debug ( 3, "Build area for line = %d, side = %d", i, side );
00531                 Vect_build_line_area ( Map, i, side );
00532             }
00533         }
00534         prnmsg (_("\r%d areas built      \n%d isles built\n"), plus->n_areas, plus->n_isles );
00535         plus->built = GV_BUILD_AREAS;
00536     }
00537     
00538     if ( build < GV_BUILD_ATTACH_ISLES ) return 1;
00539     
00540     /* Attach isles to areas */
00541     if ( plus->built < GV_BUILD_ATTACH_ISLES ) {
00542         prnmsg (_("Attaching islands: "));
00543         for (i = 1; i <= plus->n_isles; i++) {
00544             G_percent2 ( i, plus->n_isles, 1, msgout );
00545             Vect_attach_isle ( Map, i ) ;
00546         }
00547         if (i==1) prnmsg ("\n");
00548         plus->built = GV_BUILD_ATTACH_ISLES;
00549     }
00550     
00551     if ( build < GV_BUILD_CENTROIDS ) return 1;
00552 
00553     /* Attach centroids to areas */
00554     if ( plus->built < GV_BUILD_CENTROIDS ) {
00555         int nlines;
00556         
00557         prnmsg (_("Attaching centroids: "));
00558         
00559         nlines = Vect_get_num_lines (Map);
00560         for ( line = 1; line <= nlines; line++ ) {
00561             G_percent2 ( line, nlines, 1, msgout );
00562 
00563             Line = plus->Line[line];
00564             if ( !Line ) continue; /* Dead */
00565 
00566             if ( Line->type != GV_CENTROID ) continue;
00567 
00568             Node = plus->Node[Line->N1];
00569 
00570             area = Vect_find_area ( Map, Node->x, Node->y );
00571 
00572             if ( area > 0 ) {
00573                 G_debug ( 3, "Centroid (line=%d) in area %d", line, area );
00574 
00575                 Area =  plus->Area[area];
00576 
00577                 if ( Area->centroid == 0 ) { /* first */
00578                     Area->centroid = line;
00579                     Line->left = area;
00580                 } else { /* duplicate */
00581                     Line->left = -area;
00582                 }
00583             }
00584         }
00585         plus->built = GV_BUILD_CENTROIDS;
00586     }
00587 
00588     /* Add areas to category index */
00589     for ( area = 1; area <= plus->n_areas; area++ ) {
00590         int c;
00591 
00592         if ( plus->Area[area] == NULL ) continue;
00593 
00594         if ( plus->Area[area]->centroid > 0 ) { 
00595             Vect_read_line (Map, NULL, Cats, plus->Area[area]->centroid );
00596 
00597             for ( c = 0; c < Cats->n_cats; c++ ) {
00598                 dig_cidx_add_cat (plus, Cats->field[c], Cats->cat[c], area, GV_AREA);
00599             }
00600         }
00601 
00602         if ( plus->Area[area]->centroid == 0 || Cats->n_cats == 0 ) /* no centroid or no cats */
00603             dig_cidx_add_cat (plus, 0, 0, area, GV_AREA);
00604     }
00605 
00606     return 1;
00607 }
00608 
00609 

Generated on Sun Apr 6 17:32:44 2008 for GRASS by  doxygen 1.5.5