build_ogr.c

Go to the documentation of this file.
00001 /****************************************************************************
00002 *
00003 * MODULE:       Vector library 
00004 *               
00005 * AUTHOR(S):    Radim Blazek, Piero Cavalieri
00006 *
00007 * PURPOSE:      Higher level functions for reading/writing/manipulating vectors.
00008 *
00009 * COPYRIGHT:    (C) 2001 by the GRASS Development Team
00010 *
00011 *               This program is free software under the GNU General Public
00012 *               License (>=v2). Read the file COPYING that comes with GRASS
00013 *               for details.
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 *   Line offset is
00024 *      - centroids   : FID
00025 *      - other types : index of the first record (which is FID) in offset array.
00026 *
00027 *   Category: FID, not all layer have FID, OGRNullFID is defined (5/2004) as -1, so FID should be only >= 0
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 *  This structure keeps info about geometry parts above current geometry, path to curent geometry in the 
00040 *  feature. First 'part' number however is feature Id 
00041 */
00042 typedef struct {
00043     int *part;
00044     int a_parts;
00045     int n_parts;
00046 } GEOM_PARTS;
00047 
00048 /* Init parts */
00049 static void init_parts ( GEOM_PARTS *parts ) 
00050 {
00051     parts->part = NULL;
00052     parts->a_parts = parts->n_parts = 0;
00053 }
00054 
00055 /* Reset parts */
00056 static void reset_parts ( GEOM_PARTS *parts ) 
00057 {
00058     parts->n_parts = 0;
00059 }
00060 
00061 /* Free parts */
00062 static void free_parts ( GEOM_PARTS *parts ) 
00063 {
00064     free ( parts->part );
00065     parts->a_parts = parts->n_parts = 0;
00066 }
00067 
00068 /* Add new part number to parts */
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 /* Remove last part */
00080 static void del_part ( GEOM_PARTS *parts ) 
00081 {
00082     parts->n_parts--;
00083 }
00084 
00085 /* Add parts to offset */
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 /* add line to support structures */
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; /* beginning in the offset array */
00117     } else { 
00118         /* TODO : could be used to statore category ? */
00119         offset = FID; /* because centroids are read from topology, not from layer */
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     /* Set box */
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 ) /* because centroids are read from topology, not from layer */ 
00137         add_parts_to_offset ( Map, parts, line );
00138 
00139     return line;
00140 }
00141 
00142 /* Recursively add geometry to topology 
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             /* Alloc space for islands */
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                 /* register line */
00220                 add_part ( parts, iPart );
00221                 line = add_line ( Map, GV_BOUNDARY, Points[iPart], FID, parts );
00222                 del_part ( parts );
00223 
00224                 /* add area (each inner ring is also area) */
00225                 dig_line_box ( Points[iPart], &box);
00226                 dig_find_area_poly (Points[iPart], &area_size);
00227 
00228                 if ( area_size > 0 ) /* clockwise */
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                 /* Each area is also isle */
00237                 lines[0] = -lines[0]; /* island is counter clockwise */ 
00238                 
00239                 isle = dig_add_isle (plus, 1, lines);
00240                 dig_isle_set_box (plus, isle, &box);
00241                 
00242                 if ( iPart == 0 ) { /* outer ring */
00243                     outer_area = area;
00244                 } else { /* inner ring */
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             /* create virtual centroid */
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                 /* register centroid to area */
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     /* TODO move this init to better place (Vect_open_ ?), because in theory build may be reused on level2 */
00318     Map->fInfo.ogr.offset = NULL;
00319     Map->fInfo.ogr.offset_num = 0;
00320     Map->fInfo.ogr.offset_alloc = 0;
00321 
00322     /* test layer capabilities */
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     /* Note: Do not use OGR_L_GetFeatureCount (it may scan all features)!!! */
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         /* print progress */
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     }                           /* while */
00369     free_parts (&parts);
00370 
00371     prnmsg ("%4d\n", iFeature);
00372 
00373     Map->plus.built = GV_BUILD_ALL;
00374     return 1;
00375 }
00376 #endif

Generated on Sat Jul 22 22:05:57 2006 for GRASS by  doxygen 1.4.7