read_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 <grass/gis.h>
00017 #include <grass/Vect.h>
00018 
00019 #ifdef HAVE_OGR
00020 #include <ogr_api.h>
00021 
00022 /*
00023 *  Recursively read feature and add all elements to points_cache and types_cache.
00024 *  ftype : if > 0 use this type (because parts of Polygon are read as wkbLineString)
00025 * 
00026 *  return 0 OK
00027 *         1 Error
00028 */
00029 static int cache_feature ( struct Map_info *Map, OGRGeometryH hGeom, int ftype ) 
00030 {
00031     int line, i, np, ng, tp;
00032     OGRwkbGeometryType type;
00033     OGRGeometryH hGeom2;
00034 
00035     G_debug (4, "cache_feature" ); 
00036 
00037     /* Alloc space */
00038     line = Map->fInfo.ogr.lines_num;
00039     if ( line == Map->fInfo.ogr.lines_alloc ) { 
00040         Map->fInfo.ogr.lines_alloc += 20;
00041         Map->fInfo.ogr.lines = (struct line_pnts **) G_realloc ( (void *)Map->fInfo.ogr.lines, 
00042                       Map->fInfo.ogr.lines_alloc * sizeof(struct line_pnts *) );
00043 
00044         Map->fInfo.ogr.lines_types = (int *) G_realloc ( Map->fInfo.ogr.lines_types, 
00045                       Map->fInfo.ogr.lines_alloc * sizeof(int) );
00046 
00047         for ( i = Map->fInfo.ogr.lines_num; i < Map->fInfo.ogr.lines_alloc; i++ ) 
00048             Map->fInfo.ogr.lines[i] = Vect_new_line_struct();
00049         
00050     }
00051     Vect_reset_line ( Map->fInfo.ogr.lines[line] );
00052     
00053     type = wkbFlatten ( OGR_G_GetGeometryType (hGeom) );
00054     
00055     switch ( type ) {
00056         case wkbPoint:
00057             G_debug (4, "Point" );
00058             Vect_append_point ( Map->fInfo.ogr.lines[line], 
00059                                 OGR_G_GetX(hGeom,0), OGR_G_GetY(hGeom,0), OGR_G_GetZ(hGeom,0) );
00060             Map->fInfo.ogr.lines_types[line] = GV_POINT;
00061             Map->fInfo.ogr.lines_num++;
00062             return 0;
00063             break;
00064             
00065         case wkbLineString:
00066             G_debug (4, "LineString" );
00067             np = OGR_G_GetPointCount(hGeom);
00068             for( i = 0; i < np; i++ ) {
00069                 Vect_append_point ( Map->fInfo.ogr.lines[line],
00070                                     OGR_G_GetX(hGeom,i), OGR_G_GetY(hGeom,i), OGR_G_GetZ(hGeom,i) );
00071             }
00072             
00073             if ( ftype > 0 ) { /* Polygon rings */
00074                 Map->fInfo.ogr.lines_types[line] = ftype;
00075             } else { 
00076                 Map->fInfo.ogr.lines_types[line] = GV_LINE;
00077             }
00078             Map->fInfo.ogr.lines_num++;
00079             return 0;
00080             break;
00081 
00082         case wkbMultiPoint:
00083         case wkbMultiLineString:
00084         case wkbPolygon:
00085         case wkbMultiPolygon:
00086         case wkbGeometryCollection:
00087             ng = OGR_G_GetGeometryCount( hGeom );
00088             G_debug (4, "%d geoms -> next level", ng );
00089             if ( type == wkbPolygon ) {
00090                 tp = GV_BOUNDARY;
00091             } else {
00092                 tp = -1;
00093             }
00094             for ( i = 0 ; i < ng; i++ ) {
00095                 hGeom2 = OGR_G_GetGeometryRef( hGeom, i );
00096                 cache_feature ( Map, hGeom2, tp );
00097             }
00098             return 0;
00099             break;
00100 
00101         default:
00102             G_warning ("OGR feature type %d not supported\n", type);
00103             return 1;
00104             break;
00105     }
00106 }
00107 
00108 /*
00109  *  Read next line from OGR layer. Skip empty features.
00110  *
00111  *  Returns  line type
00112  *           -2 No more features (EOF)
00113  *           -1 Out of memory
00114  *
00115  *  The action of this routine can be modified by:
00116  *    Vect_read_constraint_region ()
00117  *    Vect_read_constraint_type   ()
00118  *    Vect_remove_constraints     ()
00119  */
00120 int
00121 V1_read_next_line_ogr (struct Map_info *Map, struct line_pnts *line_p, struct line_cats *line_c)
00122 {
00123     int      itype;
00124     BOUND_BOX lbox, mbox;
00125     OGRFeatureH hFeature;
00126     OGRGeometryH hGeom;
00127 
00128     G_debug (3, "V1_read_next_line_ogr()");
00129 
00130     if (line_p != NULL) Vect_reset_line (line_p);
00131     if (line_c != NULL) Vect_reset_cats (line_c);
00132 
00133     if (Map->Constraint_region_flag)
00134         Vect_get_constraint_box (Map, &mbox);
00135 
00136     while (1) {
00137         /* Read feature to chache if necessary */
00138         while ( Map->fInfo.ogr.lines_next == Map->fInfo.ogr.lines_num ) {
00139                 hFeature = OGR_L_GetNextFeature(Map->fInfo.ogr.layer);
00140 
00141                 if ( hFeature == NULL ) { return -2; } /* no more features */
00142 
00143                 hGeom = OGR_F_GetGeometryRef( hFeature );
00144                 if ( hGeom == NULL ) { /* feature without geometry */ 
00145                     OGR_F_Destroy( hFeature );
00146                     continue; 
00147                 }
00148 
00149                 Map->fInfo.ogr.feature_cache_id = (int) OGR_F_GetFID ( hFeature );
00150                 if ( Map->fInfo.ogr.feature_cache_id == OGRNullFID ) {
00151                     G_warning ("OGR feature without ID." );
00152                 }
00153 
00154                 /* Cache the feature */
00155                 Map->fInfo.ogr.lines_num = 0; 
00156                 cache_feature ( Map, hGeom, -1 );
00157                 G_debug (4, "%d lines read to cache", Map->fInfo.ogr.lines_num);
00158                 OGR_F_Destroy( hFeature );
00159 
00160                 Map->fInfo.ogr.lines_next = 0; /* next to be read from cache */
00161         }
00162 
00163         /* Read next part of the feature */
00164         G_debug ( 4, "read next cached line %d", Map->fInfo.ogr.lines_next );
00165         itype = Map->fInfo.ogr.lines_types[Map->fInfo.ogr.lines_next]; 
00166 
00167         /* Constraint on Type of line 
00168          * Default is all of  Point, Line, Area and whatever else comes along
00169          */
00170         if (Map->Constraint_type_flag) {
00171             if (!(itype & Map->Constraint_type)) {
00172                 Map->fInfo.ogr.lines_next++;
00173                 continue;
00174             }
00175         }
00176 
00177         /* Constraint on specified region */
00178         if (Map->Constraint_region_flag) {
00179             Vect_line_box ( Map->fInfo.ogr.lines[Map->fInfo.ogr.lines_next], &lbox);
00180 
00181             if (!Vect_box_overlap (&lbox, &mbox)) {
00182                 Map->fInfo.ogr.lines_next++;
00183                 continue;
00184             }
00185         }
00186 
00187         if (line_p != NULL) 
00188             Vect_append_points ( line_p, Map->fInfo.ogr.lines[Map->fInfo.ogr.lines_next], GV_FORWARD ); 
00189 
00190         if (line_c != NULL && Map->fInfo.ogr.feature_cache_id != OGRNullFID ) 
00191             Vect_cat_set ( line_c, 1, Map->fInfo.ogr.feature_cache_id );
00192 
00193         Map->fInfo.ogr.lines_next++;
00194         G_debug ( 4, "next line read, type = %d", itype );
00195         return (itype);
00196     }
00197     return -2; /* not reached */
00198 }
00199 
00200 /*
00201  *  Read next line from OGR layer.
00202  *
00203  *  Returns  line type
00204  *           -2 no more features (EOF)
00205  *           -1 Out of memory
00206  */
00207 int
00208 V2_read_next_line_ogr (struct Map_info *Map, struct line_pnts *line_p, struct line_cats *line_c)
00209 {
00210     if ( Map->next_line >= Map->plus.n_lines ) return -2;
00211 
00212     return V2_read_line_ogr (Map, line_p, line_c, Map->next_line++);
00213 }
00214 
00215 /* Recursively descend to feature and read the part
00216  * return: 0 - OK
00217  *         1 - error
00218  */
00219 static int read_line ( struct Map_info *Map, OGRGeometryH hGeom, long offset, struct line_pnts *Points )
00220 {
00221     int     i, nPoints;
00222     int     eType;
00223     OGRGeometryH hGeom2;
00224 
00225     /* Read coors if hGeom is a simple element (wkbPoint, wkbLineString) otherwise
00226      * descend to geometry specified by offset[offset] */
00227 
00228     G_debug (4, "read_line() offset = %d", offset);
00229 
00230     eType = wkbFlatten (OGR_G_GetGeometryType (hGeom));
00231     G_debug (4, "OGR Geometry of type: %d", eType);
00232 
00233     switch ( eType ) {
00234         case wkbPoint:
00235             G_debug (4, "Point" );
00236             Vect_append_point ( Points, OGR_G_GetX(hGeom,0), OGR_G_GetY(hGeom,0), OGR_G_GetZ(hGeom,0) );
00237             return 0;
00238             break;
00239 
00240         case wkbLineString:
00241             G_debug (4, "LineString" );
00242             nPoints = OGR_G_GetPointCount(hGeom);
00243             for( i = 0; i < nPoints; i++ ) {
00244                 Vect_append_point ( Points, OGR_G_GetX(hGeom,i), OGR_G_GetY(hGeom,i), OGR_G_GetZ(hGeom,i) );
00245             }
00246             return 0;
00247             break;
00248 
00249         case wkbPolygon:
00250         case wkbMultiPoint:
00251         case wkbMultiLineString:
00252         case wkbMultiPolygon:
00253         case wkbGeometryCollection:
00254             G_debug (4, " more geoms -> part %d", Map->fInfo.ogr.offset[offset] );
00255             hGeom2 = OGR_G_GetGeometryRef( hGeom, Map->fInfo.ogr.offset[offset] );
00256             return ( read_line ( Map, hGeom2, offset+1, Points ) );
00257             break;
00258 
00259         default:
00260             G_warning ("OGR feature type %d not supported\n", eType);
00261             break;
00262     }
00263     return 1;
00264 }
00265 
00266 /*
00267  *  Read line from layer on given offset.
00268  *
00269  *  Returns  line type
00270  *           -2 End of table (last row)
00271  *           -1 Out of memory
00272  */
00273 int
00274 V2_read_line_ogr (struct Map_info *Map, struct line_pnts *line_p,  struct line_cats *line_c, int line)
00275 {
00276     int  node;
00277     int  offset;
00278     long FID;
00279     P_LINE *Line;
00280     P_NODE *Node;
00281     OGRGeometryH hGeom;
00282 
00283     G_debug (4, "V2_read_line_ogr() line = %d", line);
00284 
00285     if (line_p != NULL) Vect_reset_line (line_p);
00286     if (line_c != NULL) Vect_reset_cats (line_c);
00287     
00288     Line = Map->plus.Line[line];
00289     offset = (int) Line->offset;
00290 
00291     if ( Line->type == GV_CENTROID ) {
00292         G_debug (4, "Centroid");
00293         node = Line->N1;
00294         Node = Map->plus.Node[node];
00295 
00296         /* coordinates */
00297         if (line_p != NULL) {
00298             Vect_append_point (line_p, Node->x, Node->y, 0.0);
00299         }
00300 
00301         /* category */
00302         if (line_c != NULL) {
00303             /* cat = FID and offset = FID for centroid */
00304             Vect_cat_set (line_c, 1, (int) offset); 
00305         }
00306 
00307         return (GV_CENTROID);
00308     } else {
00309         FID = Map->fInfo.ogr.offset[offset];
00310         G_debug (4, "  FID = %d", FID);
00311         
00312         /* coordinates */
00313         if (line_p != NULL) {
00314             /* Read feature to cache if necessary */
00315             if ( Map->fInfo.ogr.feature_cache_id != FID ) {
00316                 G_debug(4, "Read feature (FID = %ld) to cache.", FID);
00317                 if ( Map->fInfo.ogr.feature_cache ) {
00318                     OGR_F_Destroy( Map->fInfo.ogr.feature_cache );
00319                 }
00320                 Map->fInfo.ogr.feature_cache = OGR_L_GetFeature ( Map->fInfo.ogr.layer, FID );
00321                 if ( Map->fInfo.ogr.feature_cache == NULL ) {
00322                     G_fatal_error("Cannot read feature, FID = %ld", FID);
00323                 }
00324                 Map->fInfo.ogr.feature_cache_id = FID;
00325             }
00326                 
00327             hGeom =  OGR_F_GetGeometryRef ( Map->fInfo.ogr.feature_cache );
00328             if (hGeom == NULL) {
00329                 G_fatal_error("Cannot get feature geometry, FID = %ld", FID);
00330             }
00331             
00332             read_line ( Map, hGeom, Line->offset + 1, line_p );
00333         }
00334             
00335         /* category */
00336         if (line_c != NULL) {
00337             Vect_cat_set (line_c, 1, (int)FID);
00338         }
00339 
00340         return Line->type;
00341     }
00342     
00343     return -2; /* not reached */
00344 }
00345 
00346 #endif

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