build_ogr.c

Go to the documentation of this file.
00001 /****************************************************************************
00002 *
00003 * MODULE:       Vector library 
00004 *               
00005 * AUTHOR(S):    Radim Blazek, Piero Cavalieri
00006 *
00007 * PURPOSE:      Higher level functions for reading/writing/manipulating vectors.
00008 *
00009 * COPYRIGHT:    (C) 2001 by the GRASS Development Team
00010 *
00011 *               This program is free software under the GNU General Public
00012 *               License (>=v2). Read the file COPYING that comes with GRASS
00013 *               for details.
00014 *
00015 *****************************************************************************/
00016 #include <string.h>
00017 #include <stdlib.h>
00018 #include <stdio.h>
00019 #include <grass/gis.h>
00020 #include <grass/Vect.h>
00021 #include <grass/glocale.h>
00022 
00023 /* 
00024 *   Line offset is
00025 *      - centroids   : FID
00026 *      - other types : index of the first record (which is FID) in offset array.
00027 *
00028 *   Category: FID, not all layer have FID, OGRNullFID is defined (5/2004) as -1, so FID should be only >= 0
00029 *
00030 */
00031 
00032 
00033 #ifdef HAVE_OGR
00034 #include <ogr_api.h>
00035 
00036 extern FILE *Msgout;
00037 extern int prnmsg (char *msg, ...);
00038 
00039 /* 
00040 *  This structure keeps info about geometry parts above current geometry, path to curent geometry in the 
00041 *  feature. First 'part' number however is feature Id 
00042 */
00043 typedef struct {
00044     int *part;
00045     int a_parts;
00046     int n_parts;
00047 } GEOM_PARTS;
00048 
00049 /* Init parts */
00050 static void init_parts ( GEOM_PARTS *parts ) 
00051 {
00052     parts->part = NULL;
00053     parts->a_parts = parts->n_parts = 0;
00054 }
00055 
00056 /* Reset parts */
00057 static void reset_parts ( GEOM_PARTS *parts ) 
00058 {
00059     parts->n_parts = 0;
00060 }
00061 
00062 /* Free parts */
00063 static void free_parts ( GEOM_PARTS *parts ) 
00064 {
00065     free ( parts->part );
00066     parts->a_parts = parts->n_parts = 0;
00067 }
00068 
00069 /* Add new part number to parts */
00070 static void add_part ( GEOM_PARTS *parts, int part ) 
00071 {
00072     if ( parts->a_parts == parts->n_parts ) {
00073         parts->a_parts += 10;
00074         parts->part = (int *) G_realloc ( (void*) parts->part, parts->a_parts * sizeof(int) );
00075     }
00076     parts->part[parts->n_parts] = part;
00077     parts->n_parts++;
00078 }
00079 
00080 /* Remove last part */
00081 static void del_part ( GEOM_PARTS *parts ) 
00082 {
00083     parts->n_parts--;
00084 }
00085 
00086 /* Add parts to offset */
00087 static void add_parts_to_offset ( struct Map_info *Map, GEOM_PARTS *parts, int line ) 
00088 {
00089     int i, j;
00090 
00091     if ( Map->fInfo.ogr.offset_num + parts->n_parts >= Map->fInfo.ogr.offset_alloc ) {
00092         Map->fInfo.ogr.offset_alloc += parts->n_parts + 1000;
00093         Map->fInfo.ogr.offset = (int *) G_realloc ( Map->fInfo.ogr.offset,
00094                                                     Map->fInfo.ogr.offset_alloc * sizeof(int) );
00095     }
00096     j = Map->fInfo.ogr.offset_num;
00097     for ( i = 0; i < parts->n_parts; i++ ) {
00098         G_debug (4, "add offset %d", parts->part[i] );
00099         Map->fInfo.ogr.offset[j] = parts->part[i];
00100         j++;
00101     }
00102     Map->fInfo.ogr.offset_num += parts->n_parts;
00103 }
00104 
00105 /* add line to support structures */
00106 static int add_line ( struct Map_info *Map, int type, struct line_pnts *Points, 
00107                       int FID, GEOM_PARTS *parts )
00108 {
00109     int     line;
00110     struct  Plus_head *plus;
00111     long    offset;
00112     BOUND_BOX box;
00113 
00114     plus = &(Map->plus);
00115 
00116     if ( type != GV_CENTROID ) {
00117         offset = Map->fInfo.ogr.offset_num; /* beginning in the offset array */
00118     } else { 
00119         /* TODO : could be used to statore category ? */
00120         offset = FID; /* because centroids are read from topology, not from layer */
00121     }
00122     G_debug (4, "Register line: FID = %d offset = %ld", FID, offset );
00123     line = dig_add_line (plus, type, Points, offset );
00124     G_debug (4, "Line registered with line = %d", line);
00125 
00126     /* Set box */
00127     dig_line_box (Points, &box);
00128     if (line == 1) Vect_box_copy (&(plus->box), &box);
00129     else Vect_box_extend (&(plus->box), &box);
00130 
00131     if ( type != GV_BOUNDARY ) { 
00132         dig_cidx_add_cat (plus, 1, (int)FID, line, type);
00133     } else {
00134         dig_cidx_add_cat (plus, 0, 0, line, type);
00135     }
00136 
00137     if ( type != GV_CENTROID ) /* because centroids are read from topology, not from layer */ 
00138         add_parts_to_offset ( Map, parts, line );
00139 
00140     return line;
00141 }
00142 
00143 /* Recursively add geometry to topology 
00144 * 
00145 *  
00146 *
00147 */
00148 static int add_geometry ( struct Map_info *Map, OGRGeometryH hGeom, int FID, GEOM_PARTS *parts )
00149 {
00150     struct  Plus_head *plus;
00151     int     i, ret;
00152     int     line;
00153     int     area, isle, outer_area=0;
00154     int     lines[1];
00155     static struct line_pnts **Points = NULL;
00156     static int alloc_points = 0;
00157     BOUND_BOX box;
00158     P_LINE  *Line;
00159     double   area_size, x, y;
00160     int      eType, nRings, iPart, nParts, nPoints;
00161     OGRGeometryH hGeom2, hRing;
00162 
00163     G_debug (4, "add_geometry() FID = %d", FID );
00164     plus = &(Map->plus);
00165 
00166     if ( !Points ) {
00167         alloc_points = 1;
00168         Points = (struct line_pnts **) G_malloc ( sizeof (struct line_pnts *) );
00169         Points[0] = Vect_new_line_struct ();
00170     }
00171     Vect_reset_line ( Points[0] );
00172 
00173     eType = wkbFlatten (OGR_G_GetGeometryType (hGeom));
00174     G_debug (4, "OGR type = %d", eType);
00175 
00176     switch ( eType ) {
00177         case wkbPoint:
00178             G_debug (4, "Point" );
00179             Vect_append_point ( Points[0], OGR_G_GetX(hGeom,0), OGR_G_GetY(hGeom,0), OGR_G_GetZ(hGeom,0) );
00180             add_line ( Map, GV_POINT, Points[0], FID, parts );
00181             break;
00182 
00183         case wkbLineString:
00184             G_debug (4, "LineString" );
00185             nPoints = OGR_G_GetPointCount(hGeom);
00186             for( i = 0; i < nPoints; i++ ) {
00187                 Vect_append_point ( Points[0], 
00188                                     OGR_G_GetX(hGeom,i), OGR_G_GetY(hGeom,i), OGR_G_GetZ(hGeom,i) );
00189             }
00190             add_line ( Map, GV_LINE, Points[0], FID, parts );
00191             break;
00192 
00193         case wkbPolygon:
00194             G_debug (4, "Polygon");
00195 
00196             nRings = OGR_G_GetGeometryCount (hGeom);
00197             G_debug (4, "Number of rings: %d", nRings);
00198 
00199             /* Alloc space for islands */
00200             if ( nRings >= alloc_points ) {
00201                 Points = (struct line_pnts **) G_realloc ( (void*)Points, 
00202                                                            nRings * sizeof (struct line_pnts *));
00203                 for ( i = alloc_points; i < nRings; i++) {
00204                     Points[i] = Vect_new_line_struct ();
00205                 }
00206             }
00207 
00208             for (iPart = 0; iPart < nRings; iPart++) {
00209                 hRing = OGR_G_GetGeometryRef (hGeom, iPart);
00210                 nPoints = OGR_G_GetPointCount (hRing);
00211                 G_debug (4, "  ring %d : nPoints = %d", iPart, nPoints );
00212                 
00213                 
00214                 Vect_reset_line ( Points[iPart] );
00215                 for ( i = 0; i < nPoints; i++ ) {
00216                     Vect_append_point ( Points[iPart], 
00217                             OGR_G_GetX (hRing, i), OGR_G_GetY (hRing, i), OGR_G_GetZ (hRing, i));
00218                 }
00219 
00220                 /* register line */
00221                 add_part ( parts, iPart );
00222                 line = add_line ( Map, GV_BOUNDARY, Points[iPart], FID, parts );
00223                 del_part ( parts );
00224 
00225                 /* add area (each inner ring is also area) */
00226                 dig_line_box ( Points[iPart], &box);
00227                 dig_find_area_poly (Points[iPart], &area_size);
00228 
00229                 if ( area_size > 0 ) /* clockwise */
00230                     lines[0] = line;    
00231                 else 
00232                     lines[0] = -line;   
00233                         
00234                 area = dig_add_area (plus, 1, lines);
00235                 dig_area_set_box ( plus, area, &box );
00236 
00237                 /* Each area is also isle */
00238                 lines[0] = -lines[0]; /* island is counter clockwise */ 
00239                 
00240                 isle = dig_add_isle (plus, 1, lines);
00241                 dig_isle_set_box (plus, isle, &box);
00242                 
00243                 if ( iPart == 0 ) { /* outer ring */
00244                     outer_area = area;
00245                 } else { /* inner ring */
00246                     P_ISLE  *Isle;
00247 
00248                     Isle = plus->Isle[isle];
00249                     Isle->area = outer_area;
00250 
00251                     dig_area_add_isle ( plus, outer_area, isle );
00252                 }
00253             }
00254 
00255             /* create virtual centroid */
00256             ret = Vect_get_point_in_poly_isl ( Points[0], Points+1, nRings-1, &x, &y );
00257             if (ret < -1) {
00258                 G_warning (_("Cannot calculate centroid for area %d"), outer_area );
00259             } else { 
00260                 P_AREA  *Area;
00261 
00262                 G_debug (4, "  Centroid: %f, %f", x, y );
00263                 Vect_reset_line ( Points[0] );
00264                 Vect_append_point ( Points[0], x, y, 0.0);
00265                 line = add_line ( Map, GV_CENTROID, Points[0], FID, parts );
00266                 dig_line_box (Points[0], &box);
00267                 dig_line_set_box (plus, line, &box);
00268 
00269                 Line = plus->Line[line];
00270                 Line->left = outer_area;
00271 
00272                 /* register centroid to area */
00273                 Area = plus->Area[outer_area];
00274                 Area->centroid = line;
00275             }
00276             break;
00277 
00278         case wkbMultiPoint:
00279         case wkbMultiLineString:
00280         case wkbMultiPolygon:
00281         case wkbGeometryCollection:
00282             nParts = OGR_G_GetGeometryCount( hGeom );
00283             G_debug (4, "%d geoms -> next level", nParts );
00284             for ( i = 0 ; i < nParts; i++ ) {
00285                 add_part ( parts, i );
00286                 hGeom2 = OGR_G_GetGeometryRef( hGeom, i );
00287                 add_geometry ( Map, hGeom2, FID, parts );
00288                 del_part ( parts );
00289             }
00290             break;
00291 
00292         default:
00293             G_warning (_("OGR feature type %d not supported)"), eType);
00294             break;
00295     }
00296 
00297     return 0;
00298 }
00299 
00306 int
00307 Vect_build_ogr (struct Map_info *Map, int build, FILE * msgout)
00308 {
00309     int          iFeature, count, FID;
00310     GEOM_PARTS   parts;
00311     OGRFeatureH  hFeature;
00312     OGRGeometryH hGeom;
00313 
00314     if ( build != GV_BUILD_ALL ) G_fatal_error ("Partial build for OGR is not supported.");
00315 
00316     Msgout = msgout;
00317 
00318     /* TODO move this init to better place (Vect_open_ ?), because in theory build may be reused on level2 */
00319     Map->fInfo.ogr.offset = NULL;
00320     Map->fInfo.ogr.offset_num = 0;
00321     Map->fInfo.ogr.offset_alloc = 0;
00322 
00323     /* test layer capabilities */
00324     if ( !OGR_L_TestCapability ( Map->fInfo.ogr.layer, OLCRandomRead ) ) {
00325         G_warning ("Random read is not supported by OGR for this layer, cannot build support." );
00326         return 0;
00327     }
00328     
00329     init_parts (&parts);
00330 
00331     /* Note: Do not use OGR_L_GetFeatureCount (it may scan all features)!!! */
00332     prnmsg ("Feature: ");
00333 
00334     OGR_L_ResetReading ( Map->fInfo.ogr.layer );
00335     count = iFeature = 0;
00336     while ((hFeature = OGR_L_GetNextFeature ( Map->fInfo.ogr.layer )) != NULL) {
00337         iFeature++;
00338         count++;
00339         
00340         G_debug (4, "---- Feature %d ----", iFeature);
00341 
00342         /* print progress */
00343         if ( count == 1000 ) {
00344             prnmsg ("%7d\b\b\b\b\b\b\b", iFeature);
00345             count = 0;
00346         }
00347 
00348         hGeom = OGR_F_GetGeometryRef (hFeature);
00349         if (hGeom == NULL) {
00350             G_warning (_("Feature %d without geometry ignored"), iFeature);
00351             OGR_F_Destroy( hFeature );
00352             continue;
00353         }
00354 
00355         FID = (int) OGR_F_GetFID ( hFeature );
00356         if ( FID == OGRNullFID ) {
00357             G_warning (_("OGR feature without ID ignored."));
00358             OGR_F_Destroy( hFeature );
00359             continue;
00360         }
00361         G_debug (3, "FID =  %d", FID);
00362 
00363         reset_parts (&parts);
00364         add_part ( &parts, FID );
00365         add_geometry ( Map, hGeom, FID, &parts );
00366 
00367         OGR_F_Destroy( hFeature );
00368     }                           /* while */
00369     free_parts (&parts);
00370 
00371     prnmsg ("%4d\n", iFeature);
00372 
00373     Map->plus.built = GV_BUILD_ALL;
00374     return 1;
00375 }
00376 #endif

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