GRASS Programmer's Manual  6.4.1(2011)
build_ogr.c
Go to the documentation of this file.
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
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Defines