00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include <string.h>
00017 #include <stdlib.h>
00018 #include <stdio.h>
00019 #include <grass/gis.h>
00020 #include <grass/Vect.h>
00021 #include <grass/glocale.h>
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #ifdef HAVE_OGR
00034 #include <ogr_api.h>
00035
00036 extern FILE *Msgout;
00037 extern int prnmsg (char *msg, ...);
00038
00039
00040
00041
00042
00043 typedef struct {
00044 int *part;
00045 int a_parts;
00046 int n_parts;
00047 } GEOM_PARTS;
00048
00049
00050 static void init_parts ( GEOM_PARTS *parts )
00051 {
00052 parts->part = NULL;
00053 parts->a_parts = parts->n_parts = 0;
00054 }
00055
00056
00057 static void reset_parts ( GEOM_PARTS *parts )
00058 {
00059 parts->n_parts = 0;
00060 }
00061
00062
00063 static void free_parts ( GEOM_PARTS *parts )
00064 {
00065 free ( parts->part );
00066 parts->a_parts = parts->n_parts = 0;
00067 }
00068
00069
00070 static void add_part ( GEOM_PARTS *parts, int part )
00071 {
00072 if ( parts->a_parts == parts->n_parts ) {
00073 parts->a_parts += 10;
00074 parts->part = (int *) G_realloc ( (void*) parts->part, parts->a_parts * sizeof(int) );
00075 }
00076 parts->part[parts->n_parts] = part;
00077 parts->n_parts++;
00078 }
00079
00080
00081 static void del_part ( GEOM_PARTS *parts )
00082 {
00083 parts->n_parts--;
00084 }
00085
00086
00087 static void add_parts_to_offset ( struct Map_info *Map, GEOM_PARTS *parts, int line )
00088 {
00089 int i, j;
00090
00091 if ( Map->fInfo.ogr.offset_num + parts->n_parts >= Map->fInfo.ogr.offset_alloc ) {
00092 Map->fInfo.ogr.offset_alloc += parts->n_parts + 1000;
00093 Map->fInfo.ogr.offset = (int *) G_realloc ( Map->fInfo.ogr.offset,
00094 Map->fInfo.ogr.offset_alloc * sizeof(int) );
00095 }
00096 j = Map->fInfo.ogr.offset_num;
00097 for ( i = 0; i < parts->n_parts; i++ ) {
00098 G_debug (4, "add offset %d", parts->part[i] );
00099 Map->fInfo.ogr.offset[j] = parts->part[i];
00100 j++;
00101 }
00102 Map->fInfo.ogr.offset_num += parts->n_parts;
00103 }
00104
00105
00106 static int add_line ( struct Map_info *Map, int type, struct line_pnts *Points,
00107 int FID, GEOM_PARTS *parts )
00108 {
00109 int line;
00110 struct Plus_head *plus;
00111 long offset;
00112 BOUND_BOX box;
00113
00114 plus = &(Map->plus);
00115
00116 if ( type != GV_CENTROID ) {
00117 offset = Map->fInfo.ogr.offset_num;
00118 } else {
00119
00120 offset = FID;
00121 }
00122 G_debug (4, "Register line: FID = %d offset = %ld", FID, offset );
00123 line = dig_add_line (plus, type, Points, offset );
00124 G_debug (4, "Line registered with line = %d", line);
00125
00126
00127 dig_line_box (Points, &box);
00128 if (line == 1) Vect_box_copy (&(plus->box), &box);
00129 else Vect_box_extend (&(plus->box), &box);
00130
00131 if ( type != GV_BOUNDARY ) {
00132 dig_cidx_add_cat (plus, 1, (int)FID, line, type);
00133 } else {
00134 dig_cidx_add_cat (plus, 0, 0, line, type);
00135 }
00136
00137 if ( type != GV_CENTROID )
00138 add_parts_to_offset ( Map, parts, line );
00139
00140 return line;
00141 }
00142
00143
00144
00145
00146
00147
00148 static int add_geometry ( struct Map_info *Map, OGRGeometryH hGeom, int FID, GEOM_PARTS *parts )
00149 {
00150 struct Plus_head *plus;
00151 int i, ret;
00152 int line;
00153 int area, isle, outer_area=0;
00154 int lines[1];
00155 static struct line_pnts **Points = NULL;
00156 static int alloc_points = 0;
00157 BOUND_BOX box;
00158 P_LINE *Line;
00159 double area_size, x, y;
00160 int eType, nRings, iPart, nParts, nPoints;
00161 OGRGeometryH hGeom2, hRing;
00162
00163 G_debug (4, "add_geometry() FID = %d", FID );
00164 plus = &(Map->plus);
00165
00166 if ( !Points ) {
00167 alloc_points = 1;
00168 Points = (struct line_pnts **) G_malloc ( sizeof (struct line_pnts *) );
00169 Points[0] = Vect_new_line_struct ();
00170 }
00171 Vect_reset_line ( Points[0] );
00172
00173 eType = wkbFlatten (OGR_G_GetGeometryType (hGeom));
00174 G_debug (4, "OGR type = %d", eType);
00175
00176 switch ( eType ) {
00177 case wkbPoint:
00178 G_debug (4, "Point" );
00179 Vect_append_point ( Points[0], OGR_G_GetX(hGeom,0), OGR_G_GetY(hGeom,0), OGR_G_GetZ(hGeom,0) );
00180 add_line ( Map, GV_POINT, Points[0], FID, parts );
00181 break;
00182
00183 case wkbLineString:
00184 G_debug (4, "LineString" );
00185 nPoints = OGR_G_GetPointCount(hGeom);
00186 for( i = 0; i < nPoints; i++ ) {
00187 Vect_append_point ( Points[0],
00188 OGR_G_GetX(hGeom,i), OGR_G_GetY(hGeom,i), OGR_G_GetZ(hGeom,i) );
00189 }
00190 add_line ( Map, GV_LINE, Points[0], FID, parts );
00191 break;
00192
00193 case wkbPolygon:
00194 G_debug (4, "Polygon");
00195
00196 nRings = OGR_G_GetGeometryCount (hGeom);
00197 G_debug (4, "Number of rings: %d", nRings);
00198
00199
00200 if ( nRings >= alloc_points ) {
00201 Points = (struct line_pnts **) G_realloc ( (void*)Points,
00202 nRings * sizeof (struct line_pnts *));
00203 for ( i = alloc_points; i < nRings; i++) {
00204 Points[i] = Vect_new_line_struct ();
00205 }
00206 }
00207
00208 for (iPart = 0; iPart < nRings; iPart++) {
00209 hRing = OGR_G_GetGeometryRef (hGeom, iPart);
00210 nPoints = OGR_G_GetPointCount (hRing);
00211 G_debug (4, " ring %d : nPoints = %d", iPart, nPoints );
00212
00213
00214 Vect_reset_line ( Points[iPart] );
00215 for ( i = 0; i < nPoints; i++ ) {
00216 Vect_append_point ( Points[iPart],
00217 OGR_G_GetX (hRing, i), OGR_G_GetY (hRing, i), OGR_G_GetZ (hRing, i));
00218 }
00219
00220
00221 add_part ( parts, iPart );
00222 line = add_line ( Map, GV_BOUNDARY, Points[iPart], FID, parts );
00223 del_part ( parts );
00224
00225
00226 dig_line_box ( Points[iPart], &box);
00227 dig_find_area_poly (Points[iPart], &area_size);
00228
00229 if ( area_size > 0 )
00230 lines[0] = line;
00231 else
00232 lines[0] = -line;
00233
00234 area = dig_add_area (plus, 1, lines);
00235 dig_area_set_box ( plus, area, &box );
00236
00237
00238 lines[0] = -lines[0];
00239
00240 isle = dig_add_isle (plus, 1, lines);
00241 dig_isle_set_box (plus, isle, &box);
00242
00243 if ( iPart == 0 ) {
00244 outer_area = area;
00245 } else {
00246 P_ISLE *Isle;
00247
00248 Isle = plus->Isle[isle];
00249 Isle->area = outer_area;
00250
00251 dig_area_add_isle ( plus, outer_area, isle );
00252 }
00253 }
00254
00255
00256 ret = Vect_get_point_in_poly_isl ( Points[0], Points+1, nRings-1, &x, &y );
00257 if (ret < -1) {
00258 G_warning (_("Cannot calculate centroid for area %d"), outer_area );
00259 } else {
00260 P_AREA *Area;
00261
00262 G_debug (4, " Centroid: %f, %f", x, y );
00263 Vect_reset_line ( Points[0] );
00264 Vect_append_point ( Points[0], x, y, 0.0);
00265 line = add_line ( Map, GV_CENTROID, Points[0], FID, parts );
00266 dig_line_box (Points[0], &box);
00267 dig_line_set_box (plus, line, &box);
00268
00269 Line = plus->Line[line];
00270 Line->left = outer_area;
00271
00272
00273 Area = plus->Area[outer_area];
00274 Area->centroid = line;
00275 }
00276 break;
00277
00278 case wkbMultiPoint:
00279 case wkbMultiLineString:
00280 case wkbMultiPolygon:
00281 case wkbGeometryCollection:
00282 nParts = OGR_G_GetGeometryCount( hGeom );
00283 G_debug (4, "%d geoms -> next level", nParts );
00284 for ( i = 0 ; i < nParts; i++ ) {
00285 add_part ( parts, i );
00286 hGeom2 = OGR_G_GetGeometryRef( hGeom, i );
00287 add_geometry ( Map, hGeom2, FID, parts );
00288 del_part ( parts );
00289 }
00290 break;
00291
00292 default:
00293 G_warning (_("OGR feature type %d not supported)"), eType);
00294 break;
00295 }
00296
00297 return 0;
00298 }
00299
00306 int
00307 Vect_build_ogr (struct Map_info *Map, int build, FILE * msgout)
00308 {
00309 int iFeature, count, FID;
00310 GEOM_PARTS parts;
00311 OGRFeatureH hFeature;
00312 OGRGeometryH hGeom;
00313
00314 if ( build != GV_BUILD_ALL ) G_fatal_error ("Partial build for OGR is not supported.");
00315
00316 Msgout = msgout;
00317
00318
00319 Map->fInfo.ogr.offset = NULL;
00320 Map->fInfo.ogr.offset_num = 0;
00321 Map->fInfo.ogr.offset_alloc = 0;
00322
00323
00324 if ( !OGR_L_TestCapability ( Map->fInfo.ogr.layer, OLCRandomRead ) ) {
00325 G_warning ("Random read is not supported by OGR for this layer, cannot build support." );
00326 return 0;
00327 }
00328
00329 init_parts (&parts);
00330
00331
00332 prnmsg ("Feature: ");
00333
00334 OGR_L_ResetReading ( Map->fInfo.ogr.layer );
00335 count = iFeature = 0;
00336 while ((hFeature = OGR_L_GetNextFeature ( Map->fInfo.ogr.layer )) != NULL) {
00337 iFeature++;
00338 count++;
00339
00340 G_debug (4, "---- Feature %d ----", iFeature);
00341
00342
00343 if ( count == 1000 ) {
00344 prnmsg ("%7d\b\b\b\b\b\b\b", iFeature);
00345 count = 0;
00346 }
00347
00348 hGeom = OGR_F_GetGeometryRef (hFeature);
00349 if (hGeom == NULL) {
00350 G_warning (_("Feature %d without geometry ignored"), iFeature);
00351 OGR_F_Destroy( hFeature );
00352 continue;
00353 }
00354
00355 FID = (int) OGR_F_GetFID ( hFeature );
00356 if ( FID == OGRNullFID ) {
00357 G_warning (_("OGR feature without ID ignored."));
00358 OGR_F_Destroy( hFeature );
00359 continue;
00360 }
00361 G_debug (3, "FID = %d", FID);
00362
00363 reset_parts (&parts);
00364 add_part ( &parts, FID );
00365 add_geometry ( Map, hGeom, FID, &parts );
00366
00367 OGR_F_Destroy( hFeature );
00368 }
00369 free_parts (&parts);
00370
00371 prnmsg ("%4d\n", iFeature);
00372
00373 Map->plus.built = GV_BUILD_ALL;
00374 return 1;
00375 }
00376 #endif