00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
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
00024
00025
00026
00027
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
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 ) {
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
00110
00111
00112
00113
00114
00115
00116
00117
00118
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
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; }
00142
00143 hGeom = OGR_F_GetGeometryRef( hFeature );
00144 if ( hGeom == NULL ) {
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
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;
00161 }
00162
00163
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
00168
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
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;
00198 }
00199
00200
00201
00202
00203
00204
00205
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
00216
00217
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
00226
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
00268
00269
00270
00271
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
00297 if (line_p != NULL) {
00298 Vect_append_point (line_p, Node->x, Node->y, 0.0);
00299 }
00300
00301
00302 if (line_c != NULL) {
00303
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
00313 if (line_p != NULL) {
00314
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
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;
00344 }
00345
00346 #endif