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