GRASS Programmer's Manual 6.4.1(2011)
|
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 * Line offset is 00029 * - centroids : FID 00030 * - other types : index of the first record (which is FID) in offset array. 00031 * 00032 * Category: FID, not all layer have FID, OGRNullFID is defined (5/2004) as -1, so FID should be only >= 0 00033 * 00034 */ 00035 00036 #ifdef HAVE_OGR 00037 #include <ogr_api.h> 00038 00039 /* 00040 * This structure keeps info about geometry parts above current geometry, path to curent geometry in the 00041 * feature. First 'part' number however is feature Id 00042 */ 00043 typedef struct 00044 { 00045 int *part; 00046 int a_parts; 00047 int n_parts; 00048 } GEOM_PARTS; 00049 00050 /* Init parts */ 00051 static void init_parts(GEOM_PARTS * parts) 00052 { 00053 parts->part = NULL; 00054 parts->a_parts = parts->n_parts = 0; 00055 } 00056 00057 /* Reset parts */ 00058 static void reset_parts(GEOM_PARTS * parts) 00059 { 00060 parts->n_parts = 0; 00061 } 00062 00063 /* Free parts */ 00064 static void free_parts(GEOM_PARTS * parts) 00065 { 00066 free(parts->part); 00067 parts->a_parts = parts->n_parts = 0; 00068 } 00069 00070 /* Add new part number to parts */ 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 /* Remove last part */ 00084 static void del_part(GEOM_PARTS * parts) 00085 { 00086 parts->n_parts--; 00087 } 00088 00089 /* Add parts to offset */ 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 /* add line to support structures */ 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; /* beginning in the offset array */ 00123 } 00124 else { 00125 /* TODO : could be used to statore category ? */ 00126 offset = FID; /* because centroids are read from topology, not from layer */ 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 /* Set box */ 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) /* because centroids are read from topology, not from layer */ 00147 add_parts_to_offset(Map, parts); 00148 00149 return line; 00150 } 00151 00152 /* Recursively add geometry to topology */ 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 /* Alloc space for islands */ 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 /* register line */ 00232 add_part(parts, iPart); 00233 line = add_line(Map, GV_BOUNDARY, Points[iPart], FID, parts); 00234 del_part(parts); 00235 00236 /* add area (each inner ring is also area) */ 00237 dig_line_box(Points[iPart], &box); 00238 dig_find_area_poly(Points[iPart], &area_size); 00239 00240 if (area_size > 0) /* clockwise */ 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 /* Each area is also isle */ 00249 lines[0] = -lines[0]; /* island is counter clockwise */ 00250 00251 isle = dig_add_isle(plus, 1, lines); 00252 dig_isle_set_box(plus, isle, &box); 00253 00254 if (iPart == 0) { /* outer ring */ 00255 outer_area = area; 00256 } 00257 else { /* inner ring */ 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 /* create virtual centroid */ 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 /* register centroid to area */ 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 /* TODO move this init to better place (Vect_open_ ?), because in theory build may be reused on level2 */ 00336 Map->fInfo.ogr.offset = NULL; 00337 Map->fInfo.ogr.offset_num = 0; 00338 Map->fInfo.ogr.offset_alloc = 0; 00339 00340 /* test layer capabilities */ 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 /* Note: Do not use OGR_L_GetFeatureCount (it may scan all features)!!! */ 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 } /* while */ 00380 free_parts(&parts); 00381 00382 Map->plus.built = GV_BUILD_ALL; 00383 return 1; 00384 } 00385 #endif