vlib/pg: support reading geometry collections

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@51019 15284696-431f-4ddb-bdfa-cd5b030d7da7

include/vect/dig_structs.h

@@ -409,11 +409,42 @@ struct Format_info_offset
       GRASS. This is not used for GV_CENTROID. Because one feature may
       contain more elements (geometry collection also recursively),
       offset for one line may be stored in more records. First record
-      is FID, next records are part indexes if necessary. Example:
+      is FID, next records are part indexes if necessary.
+
+      Example 1:
+      
       5. ring in 3. polygon in 7. feature (multipolygon) of geometry
       collection which has FID = 123 123 (feature 123: geometry
       colletion) 6 (7. feature in geometry collection: multiPolygon) 2
       (3. polygon) 4 (5. ring in the polygon)
+
+      Example 2: geometry collection FID '1' containing one point, one
+      linestring and one polygon
+
+      \verbatim
+      Offset:
+
+      idx  offset note
+      ----------------
+      0	   1	  FID
+      1	   0	  first part (point)
+
+      2	   1	  FID
+      3	   1	  second part (linestring)
+
+      4	   1	  FID
+      5	   2	  third part (polygon)
+      6	   0	  first ring of polygon
+
+      Topology:
+      
+      line idx
+      -----------------
+      1    0      point
+      2    2      line
+      3    4      boundary
+      4    1      centroid read from topo (idx == FID)
+      \endverbatim
     */
     int *array;
     /*!
@@ -428,6 +459,39 @@ struct Format_info_offset
 };
 
 /*!
+  \brief Lines cache for reading feature (non-native formats)
+*/
+struct Format_info_cache {
+    /*!
+      \brief Lines array
+      
+      Some features requires more allocated lines (eg. polygon
+      with more rings, multipoint, or geometrycollection)
+    */
+    struct line_pnts **lines;
+    /*!
+      \brief List of line types (GV_POINT, GV_LINE, ...)
+    */
+    int *lines_types;
+    /*!
+      \brief Number of allocated lines in cache
+    */
+    int lines_alloc;
+    /*!
+      \brief Number of lines which forms current feature
+    */
+    int lines_num;
+    /*!
+      \brief Next line to be read from cache
+    */
+    int lines_next;
+    /*!
+      \brief Feature id
+    */
+    long fid;
+};
+
+/*!
   \brief Non-native format info (OGR)
 
   \todo Structure size should not change depending on compilation I
@@ -484,34 +548,10 @@ struct Format_info_ogr
     char **layer_options;
     
     /*!
-      \brief Lines cache (per feature)
+      \brief Lines cache for reading feature
     */
-    struct {
-	/*!
-	  \brief Lines array
-	  
-	  Some features requires more allocated lines (eg. polygon
-	  with more rings, multipoint, or geometrycollection)
-	*/
-	struct line_pnts **lines;	
-	/*!
-	  \brief List of line types
-	*/
-	int *lines_types;
-	/*!
-	  \brief Number of allocated lines in cache
-	*/
-	int lines_alloc;
-	/*!
-	  \brief Number of lines which forms current feature
-	*/
-	int lines_num;
-	/*!
-	  \brief Next line to be read from cache
-	*/
-	int lines_next;
-    } cache;
-    
+    struct Format_info_cache cache;
+
     /*!
       \brief Cache to avoid repeated reading (level 2)
 
@@ -581,39 +621,11 @@ struct Format_info_pg
       \brief Next line to be read for sequential access
     */
     int next_line;
-    
+
     /*!
-      \brief Lines cache (per feature)
+      \brief Lines cache for reading feature
     */
-    struct {
-	/*!
-	  \brief Lines array
-	  
-	  Some features requires more allocated lines (eg. polygon
-	  with more rings, multipoint, or geometrycollection)
-	*/
-	struct line_pnts **lines;
-      	/*!
-	  \brief List of line types
-	*/
-	int *lines_types;
-	/*!
-	  \brief Number of allocated lines in cache
-	*/
-	int lines_alloc;
-	/*!
-	  \brief Number of lines which forms current feature
-	*/
-	int lines_num;
-	/*!
-	  \brief Next line to be read from cache
-	*/
-	int lines_next;
-	/*!
-	  \brief Feature id
-	*/
-	long fid;
-    } cache;
+    struct Format_info_cache cache;
     
     /*!
       \brief Offset list used for building pseudo-topology

lib/vector/Vlib/build_ogr.c

@@ -79,8 +79,8 @@ int Vect_build_ogr(struct Map_info *Map, int build)
 		   "PostgreSQL") == 0)
 	    G_warning(_("Feature table <%s> has no primary key defined"),
 		      ogr_info->layer_name);
-	G_warning(_("Random read is not supported by OGR for this layer, "
-		    "unable to build topology"));
+	G_warning(_("Random read is not supported by OGR for this layer. "
+		    "Unable to build topology."));
 	return 0;
     }
 
@@ -159,6 +159,8 @@ int Vect_save_fidx(struct Map_info *Map,
 				offset->array_num, &fp))
 	return 0;
     
+    G_debug(3, "Vect_save_fidx(): offset_num = %d", offset->array_num);
+    
     fclose(fp.file);
 
     return 1;

lib/vector/Vlib/build_pg.c

@@ -70,6 +70,14 @@ int Vect_build_pg(struct Map_info *Map, int build)
 	G_warning(_("No DB connection"));
 	return 0;
     }
+
+    if (!pg_info->fid_column) {
+	G_warning(_("Feature table <%s> has no primary key defined"),
+		  pg_info->table_name);
+	G_warning(_("Random read is not supported by OGR for this layer. "
+		    "Unable to build topology."));
+	return 0;
+    }
     
     G_message(_("Using external data format '%s' (feature type '%s')"),
 	      Vect_get_finfo_format_info(Map),

lib/vector/Vlib/build_sfa.c

@@ -52,8 +52,9 @@ static int add_line(struct Plus_head *, struct Format_info_offset *,
 
 static int add_geometry_pg(struct Plus_head *,
 			   struct Format_info_pg *,
-			   SF_FeatureType, int, int,
-			   struct geom_parts *);
+			   struct feat_parts *, int, 
+			   int, int,
+ struct geom_parts *);
 static void build_pg(struct Map_info *, int);
 #endif
 
@@ -176,9 +177,10 @@ int add_line(struct Plus_head *plus, struct Format_info_offset *offset,
 	dig_cidx_add_cat(plus, 0, 0, line, type);
     }
 
-    if (type != GV_CENTROID)	/* because centroids are read from topology, not from layer */
+    /* because centroids are read from topology, not from layer */
+    if (type != GV_CENTROID)
 	add_parts_to_offset(offset, parts);
-
+    
     return line;
 }
 
@@ -188,15 +190,19 @@ int add_line(struct Plus_head *plus, struct Format_info_offset *offset,
 */
 int add_geometry_pg(struct Plus_head *plus,
 		    struct Format_info_pg *pg_info,
-		    SF_FeatureType ftype, int FID, int build,
-		    struct geom_parts *parts)
+		    struct feat_parts *fparts, int ipart,
+		    int FID, int build, struct geom_parts *parts)
 {
-    int line, iPart, area, isle, outer_area, ret;
+    int line, i, idx, area, isle, outer_area, ret;
     int lines[1];
     double area_size, x, y;
+    SF_FeatureType ftype;
     struct bound_box box;
     struct Format_info_offset *offset;
-
+    struct line_pnts *line_i;
+    
+    ftype = fparts->ftype[ipart];
+    
     G_debug(4, "add_geometry_pg() FID = %d ftype = %d", FID, ftype);
     
     offset = &(pg_info->offset);
@@ -206,28 +212,34 @@ int add_geometry_pg(struct Plus_head *plus,
     switch(ftype) {	
     case SF_POINT:
 	G_debug(4, "Point");
-	add_line(plus, offset, GV_POINT, pg_info->cache.lines[0], FID, parts);
+	line_i = pg_info->cache.lines[fparts->idx[ipart]];
+	add_line(plus, offset, GV_POINT, line_i,
+		 FID, parts);
 	break;
     case SF_LINESTRING:
 	G_debug(4, "LineString");
-	add_line(plus, offset, GV_LINE, pg_info->cache.lines[0], FID, parts);
+	line_i = pg_info->cache.lines[fparts->idx[ipart]];
+	add_line(plus, offset, GV_LINE, line_i,
+		 FID, parts);
 	break;
     case SF_POLYGON:
 	G_debug(4, "Polygon");
 
 	/* register boundaries */
-	for (iPart = 0; iPart < pg_info->cache.lines_num; iPart++) {
-	    add_part(parts, iPart);
-	    line = add_line(plus, offset, GV_BOUNDARY, pg_info->cache.lines[iPart],
-			    FID, parts);
+	idx = fparts->idx[ipart];
+	for (i = 0; i < fparts->nlines[ipart]; i++) {
+	    line_i = pg_info->cache.lines[idx++];
+	    add_part(parts, i);
+	    line = add_line(plus, offset, GV_BOUNDARY,
+			    line_i, FID, parts);
 	    del_part(parts);
 
 	    if (build < GV_BUILD_AREAS)
 		continue;
 	    
 	    /* add area (each inner ring is also area) */
-	    dig_line_box(pg_info->cache.lines[iPart], &box);
-	    dig_find_area_poly(pg_info->cache.lines[iPart], &area_size);
+	    dig_line_box(line_i, &box);
+	    dig_find_area_poly(line_i, &area_size);
 
 	    if (area_size > 0)	        /* area clockwise */
 		lines[0] = line;
@@ -244,7 +256,7 @@ int add_geometry_pg(struct Plus_head *plus,
 	    if (build < GV_BUILD_ATTACH_ISLES)
 		continue;
 	    
-	    if (iPart == 0) {	/* outer ring */
+	    if (i == 0) {	/* outer ring */
 		outer_area = area;
 	    }
 	    else {		/* inner ring */
@@ -259,9 +271,9 @@ int add_geometry_pg(struct Plus_head *plus,
 
 	if (build >= GV_BUILD_CENTROIDS) {
 	    /* create virtual centroid */
-	    ret = Vect_get_point_in_poly_isl((const struct line_pnts *) pg_info->cache.lines[0],
-					     (const struct line_pnts **) pg_info->cache.lines + 1,
-					     pg_info->cache.lines_num - 1, &x, &y);
+	    ret = Vect_get_point_in_poly_isl((const struct line_pnts *) pg_info->cache.lines[fparts->idx[ipart]],
+					     (const struct line_pnts **) pg_info->cache.lines[fparts->idx[ipart]] + 1,
+					     fparts->nlines[ipart] - 1, &x, &y);
 	    if (ret < -1) {
 		G_warning(_("Unable to calculate centroid for area %d"),
 			  outer_area);
@@ -270,11 +282,12 @@ int add_geometry_pg(struct Plus_head *plus,
 		struct P_area *Area;
 		struct P_topo_c *topo;
 		struct P_line *Line;
-		
+		struct line_pnts *line_c;
+
 		G_debug(4, "  Centroid: %f, %f", x, y);
-		Vect_reset_line(pg_info->cache.lines[0]);
-		Vect_append_point(pg_info->cache.lines[0], x, y, 0.0);
-		line = add_line(plus, offset, GV_CENTROID, pg_info->cache.lines[0], FID, parts);
+		line_c = Vect_new_line_struct();
+		Vect_append_point(line_c, x, y, 0.0);
+		line = add_line(plus, offset, GV_CENTROID, line_c, FID, parts);
 		
 		Line = plus->Line[line];
 		topo = (struct P_topo_c *)Line->topo;
@@ -283,6 +296,7 @@ int add_geometry_pg(struct Plus_head *plus,
 		/* register centroid to area */
 		Area = plus->Area[outer_area];
 		Area->centroid = line;
+		Vect_destroy_line_struct(line_c);
 	    }
 	}
 	break;
@@ -299,19 +313,21 @@ int add_geometry_pg(struct Plus_head *plus,
 */
 void build_pg(struct Map_info *Map, int build)
 {
-    int iFeature, FID, nrecords;
+    int iFeature, ipart, fid, nrecords;
     char stmt[DB_SQL_MAX];
-    SF_FeatureType ftype;
+    char *wkb_data;
     
     struct Format_info_pg *pg_info;
     
+    struct feat_parts fparts;
     struct geom_parts parts;
 
     pg_info = &(Map->fInfo.pg);
 
     /* initialize data structures */
     init_parts(&parts);
-
+    G_zero(&fparts, sizeof(struct feat_parts));
+    
     /* get records */
     sprintf(stmt, "SELECT %s,%s FROM %s",
 	    pg_info->fid_column,
@@ -330,23 +346,29 @@ void build_pg(struct Map_info *Map, int build)
     G_debug(4, "build_pg(): nrecords = %d", nrecords);
     for (iFeature = 0; iFeature < nrecords; iFeature++) {
 	/* get feature id */
-	FID  = atoi(PQgetvalue(pg_info->res, iFeature, 0));
-	/* cache feature (lines) */
-    	ftype = cache_feature(PQgetvalue(pg_info->res, iFeature, 1),
-			      FALSE, pg_info);
-	
-	if (pg_info->cache.lines_num < 1) {
-	    G_warning(_("Feature %d without geometry ignored"), iFeature);
-	    continue;
-	}
+	fid  = atoi(PQgetvalue(pg_info->res, iFeature, 0));
+	wkb_data = PQgetvalue(pg_info->res, iFeature, 1);
 	
-	G_debug(3, "Feature %d: cat=%d", iFeature, FID);
+	/* cache feature (lines) */
+	cache_feature(wkb_data,
+		      FALSE, &(pg_info->cache), &fparts);
 
 	/* register topo */
 	reset_parts(&parts);
-	add_part(&parts, FID);
-	add_geometry_pg(&(Map->plus), pg_info,
-			ftype, FID, build, &parts);
+	add_part(&parts, fid);
+	for (ipart = 0; ipart < fparts.n_parts; ipart++) {
+	    if (fparts.nlines[ipart] < 1) {
+		G_warning(_("Feature %d without geometry skipped"), fid);
+		continue;
+	    }
+	    
+	    G_debug(4, "Feature: fid = %d part = %d", fid, ipart);
+	    
+	    add_part(&parts, ipart);
+	    add_geometry_pg(&(Map->plus), pg_info, &fparts, ipart,
+			    fid, build, &parts);
+	    del_part(&parts);
+	}
     }
 
     Map->plus.built = GV_BUILD_BASE;
@@ -594,14 +616,14 @@ void build_ogr(struct Map_info *Map, int build)
 	
 	hGeom = OGR_F_GetGeometryRef(hFeature);
 	if (hGeom == NULL) {
-	    G_warning(_("Feature %d without geometry ignored"), iFeature);
+	    G_warning(_("Feature %d without geometry skipped"), iFeature);
 	    OGR_F_Destroy(hFeature);
 	    continue;
 	}
 	
 	FID = (int) OGR_F_GetFID(hFeature);
 	if (FID == OGRNullFID) {
-	    G_warning(_("OGR feature %d without ID ignored"), iFeature);
+	    G_warning(_("OGR feature %d without ID skipped"), iFeature);
 	    OGR_F_Destroy(hFeature);
 	    continue;
 	}

lib/vector/Vlib/field.c

@@ -817,7 +817,8 @@ static int read_dblinks_pg(struct Map_info *Map)
     pg_info = &(Map->fInfo.pg);
     
     if (!pg_info->fid_column) {
-	G_warning(_("No FID column defined"));
+	G_warning(_("Feature table <%s> has no primary key defined. "
+		    "Unable to define DB links."), pg_info->table_name);
 	return -1;
     }
     G_debug(3, "Using FID column <%s>", pg_info->fid_column);

lib/vector/Vlib/pg_local_proto.h

@@ -4,10 +4,22 @@
 #ifdef HAVE_POSTGRES
 #include <libpq-fe.h>
 
+/* used for building pseudo-topology (requires some extra information
+ * about lines in cache) */
+struct feat_parts
+{
+    int             a_parts; /* number of allocated items */
+    int             n_parts; /* number of parts which forms given feature */
+    SF_FeatureType *ftype;   /* simple feature type */
+    int            *nlines;  /* number of lines used in cache */
+    int            *idx;     /* index in cache where to start */
+};
+
 /* functions used in *_pg.c files */
 int execute(PGconn *, const char *);
-int cache_feature(const char *, int,
-		  struct Format_info_pg *);
+SF_FeatureType cache_feature(const char *, int,
+			     struct Format_info_cache *,
+			     struct feat_parts *);
 
 #endif /* HAVE_POSTGRES */
 

lib/vector/Vlib/poly.c

@@ -398,22 +398,21 @@ Vect_find_poly_centroid(const struct line_pnts *points,
    \brief Get point inside polygon but outside the islands specifiled in IPoints.
 
    Take a line and intersect it with the polygon and any islands.
-   sort the list of X values from these intersections.  This will
-   be a list of segments alternating  IN/OUT/IN/OUt of the polygon.
-   Pick the largest IN segment and take the midpoint. 
+   sort the list of X values from these intersections. This will be a
+   list of segments alternating IN/OUT/IN/OUT of the polygon. Pick the
+   largest IN segment and take the midpoint.
 
-   \param Points polygon
-   \param IPoints isles
+   \param Points polygon (boundary)
+   \param IPoints isles (list of isle boundaries)
    \param n_isles number of isles
    \param[out] att_x,att_y point coordinates
 
    \return 0 on success
    \return -1 on error
  */
-int
-Vect_get_point_in_poly_isl(const struct line_pnts *Points,
-			   const struct line_pnts **IPoints, int n_isles,
-			   double *att_x, double *att_y)
+int Vect_get_point_in_poly_isl(const struct line_pnts *Points,
+			       const struct line_pnts **IPoints, int n_isles,
+			       double *att_x, double *att_y)
 {
     static struct line_pnts *Intersects;
     static int first_time = 1;
@@ -771,7 +770,6 @@ int Vect_point_in_island(double X, double Y, const struct Map_info *Map,
     static struct line_pnts *Points;
     struct bound_box lbox;
     const struct Plus_head *Plus;
-    struct P_line *Line;
     struct P_isle *Isle;
 
     G_debug(3, "Vect_point_in_island(): x = %f y = %f isle = %d", X, Y, isle);
@@ -791,8 +789,6 @@ int Vect_point_in_island(double X, double Y, const struct Map_info *Map,
     for (i = 0; i < Isle->n_lines; i++) {
 	line = abs(Isle->lines[i]);
 
-	Line = Plus->Line[line];
-
 	/* this is slow, but the fastest of all alternatives */
 	Vect_get_line_box(Map, line, &lbox);
 

lib/vector/Vlib/read.c

@@ -126,8 +126,8 @@ int Vect_read_next_line(const struct Map_info *Map,
     ret = (*Read_next_line_array[Map->format][Map->level]) (Map, line_p,
 							    line_c);
     if (ret == -1)
-	G_fatal_error(_("Unable to read vector map <%s>"),
-		      Vect_get_full_name(Map));
+	G_fatal_error(_("Unable to read feature %d vector map <%s>"),
+		      Map->next_line, Vect_get_full_name(Map));
     
     return ret;
 }
@@ -168,7 +168,7 @@ int Vect_read_line(const struct Map_info *Map,
 
     ret = (*Read_line_array[Map->format]) (Map, line_p, line_c, line);
     if (ret == -1)
-	G_fatal_error(_("Unable to read feature %d vector map <%s>"),
+	G_fatal_error(_("Unable to read feature %d from vector map <%s>"),
 		      line, Vect_get_full_name(Map));
     
     return ret;

lib/vector/Vlib/read_ogr.c

@@ -180,61 +180,61 @@ int V1_read_line_ogr(struct Map_info *Map,
 		     struct line_pnts *line_p, struct line_cats *line_c, off_t offset)
 {
 #ifdef HAVE_OGR
-    long FID;
+    long fid;
     int type;
     OGRGeometryH hGeom;
 
     struct Format_info_ogr *ogr_info;
     
     ogr_info = &(Map->fInfo.ogr);
-    G_debug(4, "V1_read_line_ogr(): offset = %lu offset_num = %lu",
+    G_debug(3, "V1_read_line_ogr(): offset = %lu offset_num = %lu",
 	    (long) offset, (long) ogr_info->offset.array_num);
 
     if (offset >= ogr_info->offset.array_num)
-	return -2;
+	return -2; /* nothing to read */
     
     if (line_p != NULL)
 	Vect_reset_line(line_p);
     if (line_c != NULL)
 	Vect_reset_cats(line_c);
 
-    FID = ogr_info->offset.array[offset];
-    G_debug(4, "  FID = %ld", FID);
+    fid = ogr_info->offset.array[offset];
+    G_debug(4, "  fid = %ld", fid);
     
     /* coordinates */
     if (line_p != NULL) {
-	/* Read feature to cache if necessary */
-	if (ogr_info->feature_cache_id != FID) {
-	    G_debug(4, "Read feature (FID = %ld) to cache", FID);
+	/* read feature to cache if necessary */
+	if (ogr_info->feature_cache_id != fid) {
+	    G_debug(4, "Read feature (fid = %ld) to cache", fid);
 	    if (ogr_info->feature_cache) {
 		OGR_F_Destroy(ogr_info->feature_cache);
 	    }
 	    ogr_info->feature_cache =
-		OGR_L_GetFeature(ogr_info->layer, FID);
+		OGR_L_GetFeature(ogr_info->layer, fid);
 	    if (ogr_info->feature_cache == NULL) {
-		G_warning(_("Unable to get feature geometry, FID %ld"),
-			  FID);
+		G_warning(_("Unable to get feature geometry, fid %ld"),
+			  fid);
 		return -1;
 	    }
-	    ogr_info->feature_cache_id = FID;
+	    ogr_info->feature_cache_id = fid;
 	}
 	
 	hGeom = OGR_F_GetGeometryRef(ogr_info->feature_cache);
 	if (hGeom == NULL) {
-	    G_warning(_("Unable to get feature geometry, FID %ld"),
-		      FID);
+	    G_warning(_("Unable to get feature geometry, fid %ld"),
+		      fid);
 	    return -1;
 	}
 	
 	type = read_line(Map, hGeom, offset + 1, line_p);
     }
     else {
-	type = get_line_type(Map, FID);
+	type = get_line_type(Map, fid);
     }
 
     /* category */
     if (line_c != NULL) {
-	Vect_cat_set(line_c, 1, (int) FID);
+	Vect_cat_set(line_c, 1, (int) fid);
     }
 
     return type;
@@ -275,6 +275,7 @@ int V2_read_line_ogr(struct Map_info *Map, struct line_pnts *line_p,
     }
     
     if (Line->type == GV_CENTROID) {
+	/* read centroid for topo */
 	if (line_p != NULL) {
 	    int i, found;
 	    struct bound_box box;
@@ -288,7 +289,7 @@ int V2_read_line_ogr(struct Map_info *Map, struct line_pnts *line_p,
 		/* get area bbox */
 		Vect_get_area_box(Map, topo->area, &box);
 		/* search in spatial index for centroid with area bbox */
-		dig_init_boxlist(&list, 1);
+		dig_init_boxlist(&list, TRUE);
 		Vect_select_lines_by_box(Map, &box, Line->type, &list);
 		
 		found = 0;
@@ -304,7 +305,7 @@ int V2_read_line_ogr(struct Map_info *Map, struct line_pnts *line_p,
 	}
 
 	if (line_c != NULL) {
-	  /* cat = FID and offset = FID for centroid */
+	  /* cat = fid and offset = fid for centroid */
 	  Vect_reset_cats(line_c);
 	  Vect_cat_set(line_c, 1, (int) Line->offset);
 	}
@@ -605,7 +606,7 @@ int read_line(const struct Map_info *Map, OGRGeometryH hGeom, long offset,
   \return feature type
   \return -1 on error
 */
-int get_line_type(const struct Map_info *Map, long FID)
+int get_line_type(const struct Map_info *Map, long fid)
 {
     int eType;
 
@@ -614,11 +615,11 @@ int get_line_type(const struct Map_info *Map, long FID)
     OGRFeatureH hFeat;
     OGRGeometryH hGeom;
 
-    G_debug(4, "get_line_type() fid = %ld", FID);
+    G_debug(4, "get_line_type() fid = %ld", fid);
 
     ogr_info = &(Map->fInfo.ogr);
     
-    hFeat = OGR_L_GetFeature(ogr_info->layer, FID);
+    hFeat = OGR_L_GetFeature(ogr_info->layer, fid);
     if (hFeat == NULL)
 	return -1;
 

lib/vector/Vlib/read_pg.c

@@ -60,12 +60,17 @@ SF_FeatureType get_feature(struct Format_info_pg *, int);
 static unsigned char *hex_to_wkb(const char *, int *);
 static int point_from_wkb(const unsigned char *, int, int, int,
 			  struct line_pnts *);
-static int line_from_wkb(const unsigned char *, int, int, int,
-			 struct line_pnts *, int);
+static int linestring_from_wkb(const unsigned char *, int, int, int,
+			       struct line_pnts *, int);
 static int polygon_from_wkb(const unsigned char *, int, int, int,
-			    struct Format_info_pg *);
+			    struct Format_info_cache *);
+static int geometry_collection_from_wkb(const unsigned char *, int, int, int,
+					struct Format_info_cache *,
+					struct feat_parts *);
 static int error_corrupted_data(const char *);
 static int set_initial_query();
+static void reallocate_cache(struct Format_info_cache *, int);
+static void add_fpart(struct feat_parts *, SF_FeatureType, int, int);
 #endif
 
 /*!
@@ -170,7 +175,7 @@ int V2_read_next_line_pg(struct Map_info *Map, struct line_pnts *line_p,
 		dig_init_boxlist(&list, TRUE);
 		Vect_select_lines_by_box(Map, &box, Line->type, &list);
 		
-		found = 0;
+		found = -1;
 		for (i = 0; i < list.n_values; i++) {
 		    if (list.id[i] == line) {
 			found = i;
@@ -178,8 +183,10 @@ int V2_read_next_line_pg(struct Map_info *Map, struct line_pnts *line_p,
 		    }
 		}
 		
-		Vect_reset_line(line_p);
-		Vect_append_point(line_p, list.box[found].E, list.box[found].N, 0.0);
+		if (found > -1) {
+		    Vect_reset_line(line_p);
+		    Vect_append_point(line_p, list.box[found].E, list.box[found].N, 0.0);
+		}
 	    }
 	    if (line_c != NULL) {
 		/* cat = FID and offset = FID for centroid */
@@ -235,7 +242,7 @@ int V1_read_line_pg(struct Map_info *Map,
 {
 #ifdef HAVE_POSTGRES
     long fid;
-    int i, type;
+    int i, ipart, type;
     SF_FeatureType sf_type;
     
     struct line_pnts      *line_i;
@@ -246,7 +253,7 @@ int V1_read_line_pg(struct Map_info *Map,
 	    (long) offset, (long) pg_info->offset.array_num);
     
     if (offset >= pg_info->offset.array_num)
-	return -2;
+	return -2; /* nothing to read */
     
     if (line_p != NULL)
 	Vect_reset_line(line_p);
@@ -267,9 +274,12 @@ int V1_read_line_pg(struct Map_info *Map,
 		return (int) sf_type;
 	}
 	
+	ipart = pg_info->offset.array[offset + 1];
+	G_debug(4, "read feature part: %d", ipart);
+	
 	/* get data from cache */
-	type  = pg_info->cache.lines_types[pg_info->cache.lines_next];
-	line_i = pg_info->cache.lines[pg_info->cache.lines_next];
+	type  = pg_info->cache.lines_types[ipart];
+	line_i = pg_info->cache.lines[ipart];
 	for (i = 0; i < line_i->n_points; i++) {
 	    Vect_append_point(line_p,
 			      line_i->x[i], line_i->y[i], line_i->z[i]);
@@ -318,6 +328,7 @@ int V2_read_line_pg(struct Map_info *Map, struct line_pnts *line_p,
     }
 
     if (Line->type == GV_CENTROID) {
+	/* read centroid from topo */
 	if (line_p != NULL) {
 	    int i, found;
 	    struct bound_box box;
@@ -331,10 +342,10 @@ int V2_read_line_pg(struct Map_info *Map, struct line_pnts *line_p,
 		/* get area bbox */
 		Vect_get_area_box(Map, topo->area, &box);
 		/* search in spatial index for centroid with area bbox */
-		dig_init_boxlist(&list, 1);
+		dig_init_boxlist(&list, TRUE);
 		Vect_select_lines_by_box(Map, &box, Line->type, &list);
 		
-		found = 0;
+		found = -1;
 		for (i = 0; i < list.n_values; i++) {
 		    if (list.id[i] == line) {
 			found = i;
@@ -342,7 +353,16 @@ int V2_read_line_pg(struct Map_info *Map, struct line_pnts *line_p,
 		    }
 		}
 		
-		Vect_append_point(line_p, list.box[found].E, list.box[found].N, 0.0);
+		if (found > -1) {
+		    Vect_append_point(line_p, list.box[found].E, list.box[found].N, 0.0);
+		}
+		else {
+		    G_warning(_("Unable to construct centroid for area %d. Skipped."),
+			      topo->area);
+		}
+	    }
+	    else {
+		G_warning(_("Centroid %d: invalid area %d"), line, topo->area);
 	    }
 	}
 
@@ -411,8 +431,9 @@ int read_next_line_pg(struct Map_info *Map,
 	    if ((int) sf_type < 0) /* -1 || - 2 */
 		return (int) sf_type;
 
-	    if (sf_type & (SF_UNKNOWN | SF_NONE)) {
+	    if (sf_type == SF_UNKNOWN || sf_type == SF_NONE) {
 		G_warning(_("Feature without geometry. Skipped."));
+		pg_info->cache.lines_next = pg_info->cache.lines_num = 0;
 		continue;
 	    }
 	    
@@ -547,7 +568,7 @@ SF_FeatureType get_feature(struct Format_info_pg *pg_info, int fid)
     }
     data = (char *)PQgetvalue(pg_info->res, pg_info->next_line, 0);
     
-    ftype = cache_feature(data, FALSE, pg_info);
+    ftype = cache_feature(data, FALSE, &(pg_info->cache), NULL);
     if (fid < 0) {
 	pg_info->cache.fid = atoi(PQgetvalue(pg_info->res, pg_info->next_line, 1));
 	pg_info->next_line++;
@@ -632,14 +653,15 @@ static unsigned char *hex_to_wkb(const char *hex_data, int *nbytes)
 
   \param data HEX data
   \param skip_polygon skip polygons (level 1)
-  \param[in,out] pg_info pointer to Format_info_pg struct
-  (geometry is stored in lines cache)
+  \param[out] cache lines cache
+  \param[out] fparts used for building pseudo-topology (or NULL)
   
-  \return 0 on success
-  \return -1 on error
+  \return simple feature type
+  \return SF_UNKNOWN on error
 */
-int cache_feature(const char *data, int skip_polygon,
-		  struct Format_info_pg *pg_info)
+SF_FeatureType cache_feature(const char *data, int skip_polygon,
+			     struct Format_info_cache *cache,
+			     struct feat_parts *fparts)
 {
     int ret, byte_order, nbytes, is3D;
     unsigned char *wkb_data;
@@ -652,7 +674,7 @@ int cache_feature(const char *data, int skip_polygon,
     if (nbytes < 5) {
 	G_warning(_("Invalid WKB content: %d bytes"), nbytes);
 	G_free(wkb_data);
-	return -1;
+	return SF_UNKNOWN;
     }
     
     /* parsing M coordinate not supported */
@@ -665,7 +687,7 @@ int cache_feature(const char *data, int skip_polygon,
         G_warning(_("Reading EWKB with 4-dimensional coordinates (XYZM) "
 		    "is not supported"));
 	G_free(wkb_data);
-	return -1;
+	return SF_UNKNOWN;
     }
 
     /* PostGIS EWKB format includes an  SRID, but this won't be       
@@ -685,7 +707,7 @@ int cache_feature(const char *data, int skip_polygon,
     
     if (nbytes < 9 && nbytes != -1) {
 	G_free(wkb_data);
-	return -1;
+	return SF_UNKNOWN;
     }
     
    /* Get the geometry feature type. For now we assume that geometry
@@ -699,45 +721,52 @@ int cache_feature(const char *data, int skip_polygon,
         ftype = (SF_FeatureType) wkb_data[4];
 	is3D = wkb_data[1] & 0x80 || wkb_data[3] & 0x80;
     }
-    
+    G_debug(5, "cache_feature(): sf_type = %d", ftype);
+   
     /* allocate space in lines cache - be minimalistic
        
        more lines require eg. polygon with more rings, multi-features
        or geometry collections
     */
-    if (!pg_info->cache.lines) {
-	pg_info->cache.lines_alloc = 1;
-	pg_info->cache.lines = (struct line_pnts **) G_malloc(sizeof(struct line_pnts *));
-	
-	pg_info->cache.lines_types = (int *) G_malloc(sizeof(int));
-	pg_info->cache.lines[0] = Vect_new_line_struct();
-	pg_info->cache.lines_types[0] = -1;
+    if (!cache->lines) {
+	reallocate_cache(cache, 1);
     }
-    pg_info->cache.lines_num = 0;
+    cache->lines_num = 0;
     
     ret = -1;
     if (ftype == SF_POINT) {
-	pg_info->cache.lines_num = 1;
-	pg_info->cache.lines_types[0] = GV_POINT;
+	cache->lines_num = 1;
+	cache->lines_types[0] = GV_POINT;
 	ret = point_from_wkb(wkb_data, nbytes, byte_order,
-			     is3D, pg_info->cache.lines[0]);
+			     is3D, cache->lines[0]);
+	add_fpart(fparts, ftype, 0, 1);
     }
     else if (ftype == SF_LINESTRING) {
-	pg_info->cache.lines_num = 1;
-	pg_info->cache.lines_types[0] = GV_LINE;
-	ret = line_from_wkb(wkb_data, nbytes, byte_order,
-			    is3D, pg_info->cache.lines[0], FALSE);
+	cache->lines_num = 1;
+	cache->lines_types[0] = GV_LINE;
+	ret = linestring_from_wkb(wkb_data, nbytes, byte_order,
+				  is3D, cache->lines[0], FALSE);
+	add_fpart(fparts, ftype, 0, 1);
     }
-    else if (ftype == SF_POLYGON && !skip_polygon)
+    else if (ftype == SF_POLYGON && !skip_polygon) {
 	ret = polygon_from_wkb(wkb_data, nbytes, byte_order,
-			       is3D, pg_info) > 0 ? 0 : -1;
+			       is3D, cache);
+	add_fpart(fparts, ftype, 0, 1);
+    }
+    else if (ftype == SF_MULTIPOINT ||
+	     ftype == SF_MULTILINESTRING ||
+	     ftype == SF_MULTIPOLYGON ||
+	     ftype == SF_GEOMETRYCOLLECTION) {
+	ret = geometry_collection_from_wkb(wkb_data, nbytes, byte_order,
+					   is3D, cache, fparts);
+    }
     else  {
-	G_warning(_("Unsupported geometry type %d"), ftype);
+	G_warning(_("Unsupported feature type %d"), ftype);
     }
     
     G_free(wkb_data);
     
-    return ret;
+    return ret > 0 ? ftype : SF_UNKNOWN;
 }
 
 /*!
@@ -751,7 +780,7 @@ int cache_feature(const char *data, int skip_polygon,
   \param with_z WITH_Z for 3D data
   \param[out] line_p point geometry (pointer to line_pnts struct)
 
-  \return 0 on success
+  \return wkb size
   \return -1 on error
 */
 int point_from_wkb(const unsigned char *wkb_data, int nbytes, int byte_order,
@@ -788,7 +817,7 @@ int point_from_wkb(const unsigned char *wkb_data, int nbytes, int byte_order,
 	Vect_append_point(line_p, x, y, z);
     }
     
-    return 0;
+    return 5 + 8 * (with_z == WITH_Z ? 3 : 2);
 }
 
 /*!
@@ -802,11 +831,11 @@ int point_from_wkb(const unsigned char *wkb_data, int nbytes, int byte_order,
   \param with_z WITH_Z for 3D data
   \param[out] line_p line geometry (pointer to line_pnts struct)
 
-  \return 0 on success
+  \return wkb size
   \return -1 on error
 */
-int line_from_wkb(const unsigned char *wkb_data, int nbytes, int byte_order,
-		  int with_z, struct line_pnts *line_p, int is_ring)
+int linestring_from_wkb(const unsigned char *wkb_data, int nbytes, int byte_order,
+			int with_z, struct line_pnts *line_p, int is_ring)
 {
     int npoints, point_size, buff_min_size, offset;
     int i;
@@ -862,7 +891,7 @@ int line_from_wkb(const unsigned char *wkb_data, int nbytes, int byte_order,
 	    Vect_append_point(line_p, x, y, z);
     }
     
-    return 0;
+    return (9 - offset) + (with_z == WITH_Z ? 3 : 2) * 8 * line_p->n_points;
 }
 
 /*!
@@ -876,13 +905,13 @@ int line_from_wkb(const unsigned char *wkb_data, int nbytes, int byte_order,
   \param with_z WITH_Z for 3D data
   \param[out] line_p array of rings (pointer to line_pnts struct)
 
-  \return number of rings
+  \return wkb size
   \return -1 on error
 */
 int polygon_from_wkb(const unsigned char *wkb_data, int nbytes, int byte_order,
-		     int with_z, struct Format_info_pg *pg_info)
+		     int with_z, struct Format_info_cache *cache)
 {
-    int nrings, data_offset, i, nsize;
+    int nrings, data_offset, i, nsize, isize;
     struct line_pnts *line_i;
     
     if (nbytes < 9 && nbytes != -1)
@@ -898,21 +927,8 @@ int polygon_from_wkb(const unsigned char *wkb_data, int nbytes, int byte_order,
     }
     
     /* reallocate space for islands if needed */
-    if (nrings > pg_info->cache.lines_alloc) {
-	pg_info->cache.lines_alloc += 20;
-	pg_info->cache.lines = (struct line_pnts **) G_realloc(pg_info->cache.lines,
-							       pg_info->cache.lines_alloc *
-							       sizeof(struct line_pnts *));
-	pg_info->cache.lines_types = (int *) G_realloc(pg_info->cache.lines_types,
-						       pg_info->cache.lines_alloc *
-						       sizeof(int));
-	
-	for (i = pg_info->cache.lines_alloc - 20; i < pg_info->cache.lines_alloc; i++) {
-	    pg_info->cache.lines[i] = Vect_new_line_struct();
-	    pg_info->cache.lines_types[i] = -1;
-	}
-    }
-    pg_info->cache.lines_num = nrings;
+    reallocate_cache(cache, nrings);
+    cache->lines_num += nrings;
     
     /* each ring has a minimum of 4 bytes (point count) */
     if (nbytes != -1 && nbytes - 9 < nrings * 4) {
@@ -924,27 +940,131 @@ int polygon_from_wkb(const unsigned char *wkb_data, int nbytes, int byte_order,
         nbytes -= data_offset;
     
     /* get the rings */
-    nsize = 0;
+    nsize = 9;
     for (i = 0; i < nrings; i++ ) {
-	line_i = pg_info->cache.lines[i];
-	pg_info->cache.lines_types[i] = GV_BOUNDARY;
+	line_i = cache->lines[cache->lines_next];
+	cache->lines_types[cache->lines_next++] = GV_BOUNDARY;
 	
-	line_from_wkb(wkb_data + data_offset, nbytes, byte_order,
-		      with_z, line_i, TRUE);
+	linestring_from_wkb(wkb_data + data_offset, nbytes, byte_order,
+			    with_z, line_i, TRUE);
 	
         if (nbytes != -1) {
-	    if (with_z)
-		nsize = 4 + 24 * line_i->n_points;
-	    else
-		nsize = 4 + 16 * line_i->n_points;
-	    
+	    isize = 4 + 8 * (with_z == WITH_Z ? 3 : 2) * line_i->n_points;
+	    nbytes -= isize;
+	}
+	
+	nsize += isize;
+        data_offset += isize;
+    }
+    
+    return nsize;
+}
+
+/*!
+  \brief Read geometry collection for WKB data
+
+  See OGRGeometryCollection::importFromWkbInternal() from GDAL/OGR library
+
+  \param wkb_data WKB data
+  \param nbytes number of bytes (WKB data buffer)
+  \param byte_order byte order (ENDIAN_LITTLE, ENDIAN_BIG)
+  \param with_z WITH_Z for 3D data
+  \param ipart part to cache (starts at 0)
+  \param[out] cache lines cache
+  \param[in,out] fparts feature parts (required for building pseudo-topology)
+
+  \return number of parts
+  \return -1 on error
+*/
+int geometry_collection_from_wkb(const unsigned char *wkb_data, int nbytes, int byte_order,
+				 int with_z, struct Format_info_cache *cache,
+				 struct feat_parts *fparts)
+{
+    int ipart, nparts, data_offset, nsize;
+    unsigned char *wkb_subdata;
+    SF_FeatureType ftype;
+    
+    if (nbytes < 9 && nbytes != -1)
+       return error_corrupted_data(NULL);
+    
+    /* get the geometry count */
+    memcpy(&nparts, wkb_data + 5, 4 );
+    if (byte_order == ENDIAN_BIG) {
+        nparts = SWAP32(nparts);
+    }
+    if (nparts < 0 || nparts > INT_MAX / 9) {
+        return error_corrupted_data(NULL);
+    }
+    G_debug(5, "\t(geometry collections) parts: %d", nparts);
+    
+    /* each geometry has a minimum of 9 bytes */
+    if (nbytes != -1 && nbytes - 9 < nparts * 9) {
+        return error_corrupted_data(_("Length of input WKB is too small"));
+    }
+
+    data_offset = 9;
+    if (nbytes != -1)
+        nbytes -= data_offset;
+
+    /* reallocate space for parts if needed */
+    reallocate_cache(cache, nparts);
+    
+    /* get parts */
+    cache->lines_next = cache->lines_num = 0;
+    for (ipart = 0; ipart < nparts; ipart++) {
+	wkb_subdata = (unsigned char *)wkb_data + data_offset;
+	if (nbytes < 9 && nbytes != -1)
+	    return error_corrupted_data(NULL);
+	
+	if (byte_order == ENDIAN_LITTLE) {
+	    ftype = (SF_FeatureType) wkb_subdata[1];
+	}
+	else {
+	    ftype = (SF_FeatureType) wkb_subdata[4];
+	}
+	
+	if (ftype == SF_POINT) {
+	    cache->lines_types[cache->lines_next] = GV_POINT;
+	    nsize = point_from_wkb(wkb_subdata, nbytes, byte_order, with_z,
+				   cache->lines[cache->lines_next]);
+	    cache->lines_num++;
+	    add_fpart(fparts, ftype, cache->lines_next, 1);
+	    cache->lines_next++;
+	}
+	else if (ftype == SF_LINESTRING) {
+	    cache->lines_types[cache->lines_next] = GV_LINE;
+	    nsize = linestring_from_wkb(wkb_subdata, nbytes, byte_order, with_z,
+					cache->lines[cache->lines_next],
+					FALSE);
+	    cache->lines_num++;
+	    add_fpart(fparts, ftype, cache->lines_next, 1);
+	    cache->lines_next++;
+	}
+	else if (ftype == SF_POLYGON) {
+	    int idx = cache->lines_next;
+	    nsize = polygon_from_wkb(wkb_subdata, nbytes, byte_order,
+				     with_z, cache);
+	    add_fpart(fparts, ftype, idx, cache->lines_num - idx);
+	}
+	else if (ftype == SF_GEOMETRYCOLLECTION ||
+		 ftype == SF_MULTIPOLYGON ||
+		 ftype == SF_MULTILINESTRING ||
+		 ftype == SF_MULTIPOLYGON) {
+	    // geometry_collection_from_wkb();
+	}
+	else  {
+	    G_warning(_("Unsupported feature type %d"), ftype);
+	}
+
+        if (nbytes != -1) {
 	    nbytes -= nsize;
 	}
 	
         data_offset += nsize;
     }
+    cache->lines_next = 0;
     
-    return nrings;
+    return nparts;
 }
 
 /*!
@@ -1025,4 +1145,67 @@ int execute(PGconn *conn, const char *stmt)
     return 0;
 }
 
+/*!
+  \brief Reallocate lines cache
+*/
+void reallocate_cache(struct Format_info_cache *cache, int num)
+{
+    int i;
+    
+    if (cache->lines_alloc >= num)
+	return;
+
+    if (!cache->lines) {
+	/* most of features requires only one line cache */
+	cache->lines_alloc = 1; 
+    }
+    else {
+	cache->lines_alloc += 20;
+    }
+
+    cache->lines = (struct line_pnts **) G_realloc(cache->lines,
+						   cache->lines_alloc *
+						   sizeof(struct line_pnts *));
+    cache->lines_types = (int *) G_realloc(cache->lines_types,
+					   cache->lines_alloc *
+					   sizeof(int));
+    
+    if (cache->lines_alloc > 1) {
+	for (i = cache->lines_alloc - 20; i < cache->lines_alloc; i++) {
+	    cache->lines[i] = Vect_new_line_struct();
+	    cache->lines_types[i] = -1;
+	}
+    }
+    else {
+	cache->lines[0] = Vect_new_line_struct();
+	cache->lines_types[0] = -1;
+    }
+}
+
+void add_fpart(struct feat_parts *fparts, SF_FeatureType ftype,
+	       int idx, int nlines)
+{
+    if (!fparts)
+	return;
+    
+    if (fparts->a_parts == 0 || fparts->n_parts >= fparts->a_parts) {
+	if (fparts->a_parts == 0)
+	    fparts->a_parts = 1;
+	else
+	    fparts->a_parts += 20;
+	
+	fparts->ftype  = (SF_FeatureType *) G_realloc(fparts->ftype,
+						      fparts->a_parts * sizeof(SF_FeatureType));
+	fparts->nlines = (int *) G_realloc(fparts->nlines,
+					   fparts->a_parts * sizeof(int));
+	fparts->idx    = (int *) G_realloc(fparts->idx,
+					   fparts->a_parts * sizeof(int));
+    }
+    
+    fparts->ftype[fparts->n_parts]  = ftype;
+    fparts->idx[fparts->n_parts]    = idx;
+    fparts->nlines[fparts->n_parts] = nlines;
+    
+    fparts->n_parts++;
+}
 #endif