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
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) {
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
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 }
00169
00170 hGeom = OGR_F_GetGeometryRef(hFeature);
00171 if (hGeom == NULL) {
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
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;
00188 }
00189
00190
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
00195
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
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;
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
00270
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
00351 if (line_p != NULL) {
00352 Vect_append_point(line_p, Node->x, Node->y, 0.0);
00353 }
00354
00355
00356 if (line_c != NULL) {
00357
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
00368 if (line_p != NULL) {
00369
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
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;
00402 }
00403
00404 #endif