read_ogr.c

Go to the documentation of this file.
00001 
00020 #include <grass/gis.h>
00021 #include <grass/Vect.h>
00022 #include <grass/glocale.h>
00023 
00024 #ifdef HAVE_OGR
00025 #include <ogr_api.h>
00026 
00039 static int cache_feature(struct Map_info *Map, OGRGeometryH hGeom, int ftype)
00040 {
00041     int line, i, np, ng, tp;
00042     OGRwkbGeometryType type;
00043     OGRGeometryH hGeom2;
00044 
00045     G_debug(4, "cache_feature");
00046 
00047     /* Alloc space */
00048     line = Map->fInfo.ogr.lines_num;
00049     if (line == Map->fInfo.ogr.lines_alloc) {
00050         Map->fInfo.ogr.lines_alloc += 20;
00051         Map->fInfo.ogr.lines =
00052             (struct line_pnts **)G_realloc((void *)Map->fInfo.ogr.lines,
00053                                            Map->fInfo.ogr.lines_alloc *
00054                                            sizeof(struct line_pnts *));
00055 
00056         Map->fInfo.ogr.lines_types =
00057             (int *)G_realloc(Map->fInfo.ogr.lines_types,
00058                              Map->fInfo.ogr.lines_alloc * sizeof(int));
00059 
00060         for (i = Map->fInfo.ogr.lines_num; i < Map->fInfo.ogr.lines_alloc;
00061              i++)
00062             Map->fInfo.ogr.lines[i] = Vect_new_line_struct();
00063 
00064     }
00065     Vect_reset_line(Map->fInfo.ogr.lines[line]);
00066 
00067     type = wkbFlatten(OGR_G_GetGeometryType(hGeom));
00068 
00069     switch (type) {
00070     case wkbPoint:
00071         G_debug(4, "Point");
00072         Vect_append_point(Map->fInfo.ogr.lines[line],
00073                           OGR_G_GetX(hGeom, 0), OGR_G_GetY(hGeom, 0),
00074                           OGR_G_GetZ(hGeom, 0));
00075         Map->fInfo.ogr.lines_types[line] = GV_POINT;
00076         Map->fInfo.ogr.lines_num++;
00077         return 0;
00078         break;
00079 
00080     case wkbLineString:
00081         G_debug(4, "LineString");
00082         np = OGR_G_GetPointCount(hGeom);
00083         for (i = 0; i < np; i++) {
00084             Vect_append_point(Map->fInfo.ogr.lines[line],
00085                               OGR_G_GetX(hGeom, i), OGR_G_GetY(hGeom, i),
00086                               OGR_G_GetZ(hGeom, i));
00087         }
00088 
00089         if (ftype > 0) {        /* Polygon rings */
00090             Map->fInfo.ogr.lines_types[line] = ftype;
00091         }
00092         else {
00093             Map->fInfo.ogr.lines_types[line] = GV_LINE;
00094         }
00095         Map->fInfo.ogr.lines_num++;
00096         return 0;
00097         break;
00098 
00099     case wkbMultiPoint:
00100     case wkbMultiLineString:
00101     case wkbPolygon:
00102     case wkbMultiPolygon:
00103     case wkbGeometryCollection:
00104         ng = OGR_G_GetGeometryCount(hGeom);
00105         G_debug(4, "%d geoms -> next level", ng);
00106         if (type == wkbPolygon) {
00107             tp = GV_BOUNDARY;
00108         }
00109         else {
00110             tp = -1;
00111         }
00112         for (i = 0; i < ng; i++) {
00113             hGeom2 = OGR_G_GetGeometryRef(hGeom, i);
00114             cache_feature(Map, hGeom2, tp);
00115         }
00116         return 0;
00117         break;
00118 
00119     default:
00120         G_warning(_("OGR feature type %d not supported"), type);
00121         return 1;
00122         break;
00123     }
00124 }
00125 
00142 int
00143 V1_read_next_line_ogr(struct Map_info *Map, struct line_pnts *line_p,
00144                       struct line_cats *line_c)
00145 {
00146     int itype;
00147     BOUND_BOX lbox, mbox;
00148     OGRFeatureH hFeature;
00149     OGRGeometryH hGeom;
00150 
00151     G_debug(3, "V1_read_next_line_ogr()");
00152 
00153     if (line_p != NULL)
00154         Vect_reset_line(line_p);
00155     if (line_c != NULL)
00156         Vect_reset_cats(line_c);
00157 
00158     if (Map->Constraint_region_flag)
00159         Vect_get_constraint_box(Map, &mbox);
00160 
00161     while (1) {
00162         /* Read feature to chache if necessary */
00163         while (Map->fInfo.ogr.lines_next == Map->fInfo.ogr.lines_num) {
00164             hFeature = OGR_L_GetNextFeature(Map->fInfo.ogr.layer);
00165 
00166             if (hFeature == NULL) {
00167                 return -2;
00168             }                   /* no more features */
00169 
00170             hGeom = OGR_F_GetGeometryRef(hFeature);
00171             if (hGeom == NULL) {        /* feature without geometry */
00172                 OGR_F_Destroy(hFeature);
00173                 continue;
00174             }
00175 
00176             Map->fInfo.ogr.feature_cache_id = (int)OGR_F_GetFID(hFeature);
00177             if (Map->fInfo.ogr.feature_cache_id == OGRNullFID) {
00178                 G_warning(_("OGR feature without ID"));
00179             }
00180 
00181             /* Cache the feature */
00182             Map->fInfo.ogr.lines_num = 0;
00183             cache_feature(Map, hGeom, -1);
00184             G_debug(4, "%d lines read to cache", Map->fInfo.ogr.lines_num);
00185             OGR_F_Destroy(hFeature);
00186 
00187             Map->fInfo.ogr.lines_next = 0;      /* next to be read from cache */
00188         }
00189 
00190         /* Read next part of the feature */
00191         G_debug(4, "read next cached line %d", Map->fInfo.ogr.lines_next);
00192         itype = Map->fInfo.ogr.lines_types[Map->fInfo.ogr.lines_next];
00193 
00194         /* Constraint on Type of line 
00195          * Default is all of  Point, Line, Area and whatever else comes along
00196          */
00197         if (Map->Constraint_type_flag) {
00198             if (!(itype & Map->Constraint_type)) {
00199                 Map->fInfo.ogr.lines_next++;
00200                 continue;
00201             }
00202         }
00203 
00204         /* Constraint on specified region */
00205         if (Map->Constraint_region_flag) {
00206             Vect_line_box(Map->fInfo.ogr.lines[Map->fInfo.ogr.lines_next],
00207                           &lbox);
00208 
00209             if (!Vect_box_overlap(&lbox, &mbox)) {
00210                 Map->fInfo.ogr.lines_next++;
00211                 continue;
00212             }
00213         }
00214 
00215         if (line_p != NULL)
00216             Vect_append_points(line_p,
00217                                Map->fInfo.ogr.lines[Map->fInfo.ogr.
00218                                                     lines_next], GV_FORWARD);
00219 
00220         if (line_c != NULL && Map->fInfo.ogr.feature_cache_id != OGRNullFID)
00221             Vect_cat_set(line_c, 1, Map->fInfo.ogr.feature_cache_id);
00222 
00223         Map->fInfo.ogr.lines_next++;
00224         G_debug(4, "next line read, type = %d", itype);
00225         return (itype);
00226     }
00227     return -2;                  /* not reached */
00228 }
00229 
00241 int
00242 V2_read_next_line_ogr(struct Map_info *Map, struct line_pnts *line_p,
00243                       struct line_cats *line_c)
00244 {
00245     if (Map->next_line >= Map->plus.n_lines)
00246         return -2;
00247 
00248     return V2_read_line_ogr(Map, line_p, line_c, Map->next_line++);
00249 }
00250 
00262 static int read_line(struct Map_info *Map, OGRGeometryH hGeom, long offset,
00263                      struct line_pnts *Points)
00264 {
00265     int i, nPoints;
00266     int eType;
00267     OGRGeometryH hGeom2;
00268 
00269     /* Read coors if hGeom is a simple element (wkbPoint, wkbLineString) otherwise
00270      * descend to geometry specified by offset[offset] */
00271 
00272     G_debug(4, "read_line() offset = %ld", offset);
00273 
00274     eType = wkbFlatten(OGR_G_GetGeometryType(hGeom));
00275     G_debug(4, "OGR Geometry of type: %d", eType);
00276 
00277     switch (eType) {
00278     case wkbPoint:
00279         G_debug(4, "Point");
00280         Vect_append_point(Points, OGR_G_GetX(hGeom, 0), OGR_G_GetY(hGeom, 0),
00281                           OGR_G_GetZ(hGeom, 0));
00282         return 0;
00283         break;
00284 
00285     case wkbLineString:
00286         G_debug(4, "LineString");
00287         nPoints = OGR_G_GetPointCount(hGeom);
00288         for (i = 0; i < nPoints; i++) {
00289             Vect_append_point(Points, OGR_G_GetX(hGeom, i),
00290                               OGR_G_GetY(hGeom, i), OGR_G_GetZ(hGeom, i));
00291         }
00292         return 0;
00293         break;
00294 
00295     case wkbPolygon:
00296     case wkbMultiPoint:
00297     case wkbMultiLineString:
00298     case wkbMultiPolygon:
00299     case wkbGeometryCollection:
00300         G_debug(4, " more geoms -> part %d", Map->fInfo.ogr.offset[offset]);
00301         hGeom2 = OGR_G_GetGeometryRef(hGeom, Map->fInfo.ogr.offset[offset]);
00302         return (read_line(Map, hGeom2, offset + 1, Points));
00303         break;
00304 
00305     default:
00306         G_warning(_("OGR feature type %d not supported"), eType);
00307         break;
00308     }
00309     return 1;
00310 }
00311 
00324 int
00325 V2_read_line_ogr(struct Map_info *Map, struct line_pnts *line_p,
00326                  struct line_cats *line_c, int line)
00327 {
00328     int node;
00329     int offset;
00330     long FID;
00331     P_LINE *Line;
00332     P_NODE *Node;
00333     OGRGeometryH hGeom;
00334 
00335     G_debug(4, "V2_read_line_ogr() line = %d", line);
00336 
00337     if (line_p != NULL)
00338         Vect_reset_line(line_p);
00339     if (line_c != NULL)
00340         Vect_reset_cats(line_c);
00341 
00342     Line = Map->plus.Line[line];
00343     offset = (int)Line->offset;
00344 
00345     if (Line->type == GV_CENTROID) {
00346         G_debug(4, "Centroid");
00347         node = Line->N1;
00348         Node = Map->plus.Node[node];
00349 
00350         /* coordinates */
00351         if (line_p != NULL) {
00352             Vect_append_point(line_p, Node->x, Node->y, 0.0);
00353         }
00354 
00355         /* category */
00356         if (line_c != NULL) {
00357             /* cat = FID and offset = FID for centroid */
00358             Vect_cat_set(line_c, 1, (int)offset);
00359         }
00360 
00361         return (GV_CENTROID);
00362     }
00363     else {
00364         FID = Map->fInfo.ogr.offset[offset];
00365         G_debug(4, "  FID = %ld", FID);
00366 
00367         /* coordinates */
00368         if (line_p != NULL) {
00369             /* Read feature to cache if necessary */
00370             if (Map->fInfo.ogr.feature_cache_id != FID) {
00371                 G_debug(4, "Read feature (FID = %ld) to cache.", FID);
00372                 if (Map->fInfo.ogr.feature_cache) {
00373                     OGR_F_Destroy(Map->fInfo.ogr.feature_cache);
00374                 }
00375                 Map->fInfo.ogr.feature_cache =
00376                     OGR_L_GetFeature(Map->fInfo.ogr.layer, FID);
00377                 if (Map->fInfo.ogr.feature_cache == NULL) {
00378                     G_fatal_error(_("Unable to get feature geometry, FID %ld"),
00379                                   FID);
00380                 }
00381                 Map->fInfo.ogr.feature_cache_id = FID;
00382             }
00383 
00384             hGeom = OGR_F_GetGeometryRef(Map->fInfo.ogr.feature_cache);
00385             if (hGeom == NULL) {
00386                 G_fatal_error(_("Unable to get feature geometry, FID %ld"),
00387                               FID);
00388             }
00389 
00390             read_line(Map, hGeom, Line->offset + 1, line_p);
00391         }
00392 
00393         /* category */
00394         if (line_c != NULL) {
00395             Vect_cat_set(line_c, 1, (int)FID);
00396         }
00397 
00398         return Line->type;
00399     }
00400 
00401     return -2;                  /* not reached */
00402 }
00403 
00404 #endif

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