GRASS Programmer's Manual
6.4.1(2011)
|
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