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