vlib: implement V2__add_line_to_topo_ogr()

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

lib/vector/Vlib/build_nat.c

@@ -34,7 +34,7 @@
  */
 int Vect_build_line_area(struct Map_info *Map, int iline, int side)
 {
-    int j, area, isle, n_lines, line, type, direction;
+    int j, area, isle, n_lines, line, direction;
     static int first = 1;
     off_t offset;
     struct Plus_head *plus;
@@ -74,7 +74,7 @@ int Vect_build_line_area(struct Map_info *Map, int iline, int side)
 	offset = BLine->offset;
 	G_debug(3, "  line[%d] = %d, offset = %lu", j, line,
 		(unsigned long)offset);
-	type = Vect_read_line(Map, Points, NULL, line);
+	Vect_read_line(Map, Points, NULL, line);
 	if (lines[j] > 0)
 	    direction = GV_FORWARD;
 	else
@@ -132,7 +132,7 @@ int Vect_build_line_area(struct Map_info *Map, int iline, int side)
  */
 int Vect_isle_find_area(struct Map_info *Map, int isle)
 {
-    int j, line, sel_area, first, area, poly;
+    int j, line, sel_area, area, poly;
     static int first_call = 1;
     const struct Plus_head *plus;
     struct P_line *Line;
@@ -181,7 +181,6 @@ int Vect_isle_find_area(struct Map_info *Map, int isle)
 
     sel_area = 0;
     cur_size = -1;
-    first = 1;
     Vect_get_isle_box(Map, isle, &box);
     for (j = 0; j < List->n_values; j++) {
 	area = List->id[j];
@@ -275,20 +274,20 @@ int Vect_attach_isle(struct Map_info *Map, int isle)
     struct P_isle *Isle;
     struct Plus_head *plus;
 
-    /* Note!: If topology is not clean and areas overlap, one island may fall to more areas
-     *  (partially or fully). Before isle is attached to area it must be check if it is not attached yet */
-    G_debug(3, "Vect_attach_isle (): isle = %d", isle);
+    /* Note!: If topology is not clean and areas overlap, one island
+       may fall to more areas (partially or fully). Before isle is
+       attached to area it must be check if it is not attached yet */
+    G_debug(3, "Vect_attach_isle(): isle = %d", isle);
 
     plus = &(Map->plus);
 
     sel_area = Vect_isle_find_area(Map, isle);
-    G_debug(3, "      isle = %d -> area outside = %d", isle, sel_area);
+    G_debug(3, "\tisle = %d -> area outside = %d", isle, sel_area);
     if (sel_area > 0) {
 	Isle = plus->Isle[isle];
 	if (Isle->area > 0) {
-	    G_debug(3,
-		    "Attempt to attach isle %d to more areas (=>topology is not clean)",
-		    isle);
+	    G_debug(3, "Attempt to attach isle %d to more areas "
+		    "(=>topology is not clean)", isle);
 	}
 	else {
 	    Isle->area = sel_area;
@@ -306,20 +305,20 @@ int Vect_attach_isle(struct Map_info *Map, int isle)
 
    \return 0
  */
-int Vect_attach_isles(struct Map_info *Map, const struct bound_box * box)
+int Vect_attach_isles(struct Map_info *Map, const struct bound_box *box)
 {
     int i, isle;
-    static int first = 1;
+    static int first = TRUE;
     static struct boxlist *List;
     struct Plus_head *plus;
 
-    G_debug(3, "Vect_attach_isles ()");
-
+    G_debug(3, "Vect_attach_isles()");
+      
     plus = &(Map->plus);
 
     if (first) {
-	List = Vect_new_boxlist(0);
-	first = 0;
+	List = Vect_new_boxlist(FALSE);
+	first = FALSE;
     }
 
     Vect_select_isles_by_box(Map, box, List);
@@ -337,6 +336,34 @@ int Vect_attach_isles(struct Map_info *Map, const struct bound_box * box)
 /*!
    \brief (Re)Attach centroids to areas in given bounding box
 
+    Warning: If map is updated on level2, it may happen that
+    previously correct island becomes incorrect. In that case,
+    centroid of area forming the island is reattached to outer area,
+    because island polygon is not excluded.
+     
+    <pre>
+      +-----------+     +-----------+
+      |   1       |     |   1       |
+      | +---+---+ |     | +---+---+ |     
+      | | 2 | 3 | |     | | 2 |     |   
+      | | x |   | |  -> | | x |     |  
+      | |   |   | |     | |   |     | 
+      | +---+---+ |     | +---+---+ |
+      |           |     |           |
+      +-----------+     +-----------+
+      centroid is       centroid is
+      attached to 2     reattached to 1
+    </pre>
+
+    Because of this, when the centroid is reattached to another area,
+    it is always necessary to check if original area exist, unregister
+    centroid from previous area.  To simplify code, this is
+    implemented so that centroid is always firs unregistered and if
+    new area is found, it is registered again.
+      
+    This problem can be avoided altogether if properly attached
+    centroids are skipped MM 2009
+
    \param Map_info vector map
    \param box bounding box
 
@@ -353,7 +380,7 @@ int Vect_attach_centroids(struct Map_info *Map, const struct bound_box * box)
     struct P_topo_c *topo;
     struct Plus_head *plus;
 
-    G_debug(3, "Vect_attach_centroids ()");
+    G_debug(3, "Vect_attach_centroids()");
 
     plus = &(Map->plus);
 
@@ -363,34 +390,8 @@ int Vect_attach_centroids(struct Map_info *Map, const struct bound_box * box)
 	first = 0;
     }
 
-    /* Warning: If map is updated on level2, it may happen that previously correct island 
-     * becomes incorrect. In that case, centroid of area forming the island is reattached 
-     * to outer area, because island polygon is not excluded. 
-     *
-     * +-----------+     +-----------+
-     * |   1       |     |   1       |
-     * | +---+---+ |     | +---+---+ |     
-     * | | 2 | 3 | |     | | 2 |     |   
-     * | | x |   | |  -> | | x |     |  
-     * | |   |   | |     | |   |     | 
-     * | +---+---+ |     | +---+---+ |
-     * |           |     |           |
-     * +-----------+     +-----------+
-     * centroid is       centroid is
-     * attached to 2     reattached to 1
-     *
-     * Because of this, when the centroid is reattached to another area, it is always necessary
-     * to check if original area exist, unregister centroid from previous area.
-     * To simplify code, this is implemented so that centroid is always firs unregistered 
-     * and if new area is found, it is registered again.
-     *
-     * This problem can be avoided altogether if properly attached centroids
-     * are skipped
-     * MM 2009
-     */
-
     Vect_select_lines_by_box(Map, box, GV_CENTROID, List);
-    G_debug(3, "  number of centroids to reattach = %d", List->n_values);
+    G_debug(3, "\tnumber of centroids to reattach = %d", List->n_values);
     for (i = 0; i < List->n_values; i++) {
 	int orig_area;
 
@@ -409,12 +410,28 @@ int Vect_attach_centroids(struct Map_info *Map, const struct bound_box * box)
 	orig_area = topo->area;
 
 	Vect_read_line(Map, Points, NULL, centr);
+	if (Points->n_points < 1) {
+	    /* try to get centroid from spatial index (OGR layers) */
+	    int found;
+	    struct boxlist list;
+	    dig_init_boxlist(&list, TRUE);
+	    Vect_select_lines_by_box(Map, box, GV_CENTROID, &list);
+
+	    found = 0;
+	    for (i = 0; i < list.n_values; i++) {
+		if (list.id[i] == centr) {
+		    found = i;
+		    break;
+		}
+	    }
+	    Vect_append_point(Points, list.box[found].E, list.box[found].N, 0.0);
+	}
 	sel_area = Vect_find_area(Map, Points->x[0], Points->y[0]);
-	G_debug(3, "  centroid %d is in area %d", centr, sel_area);
+	G_debug(3, "\tcentroid %d is in area %d", centr, sel_area);
 	if (sel_area > 0) {
 	    Area = plus->Area[sel_area];
 	    if (Area->centroid == 0) {	/* first centroid */
-		G_debug(3, "  first centroid -> attach to area");
+		G_debug(3, "\tfirst centroid -> attach to area");
 		Area->centroid = centr;
 		topo->area = sel_area;
 
@@ -424,7 +441,7 @@ int Vect_attach_centroids(struct Map_info *Map, const struct bound_box * box)
 	    else if (Area->centroid != centr) {	/* duplicate centroid */
 		/* Note: it cannot happen that Area->centroid == centr, because the centroid
 		 * was not registered or a duplicate */
-		G_debug(3, "  duplicate centroid -> do not attach to area");
+		G_debug(3, "\tduplicate centroid -> do not attach to area");
 		topo->area = -sel_area;
 
 		if (-sel_area != orig_area && plus->do_uplist)
@@ -451,14 +468,12 @@ int Vect_build_nat(struct Map_info *Map, int build)
     int i, s, type, line;
     off_t offset;
     int side, area;
-    struct line_pnts *Points, *APoints;
+    struct line_pnts *Points;
     struct line_cats *Cats;
     struct P_line *Line;
     struct P_area *Area;
     struct bound_box box;
-    struct ilist *List;
-    int print_counter = G_verbose() > G_verbose_min();
-
+    
     G_debug(3, "Vect_build_nat() build = %d", build);
 
     plus = &(Map->plus);
@@ -517,14 +532,10 @@ int Vect_build_nat(struct Map_info *Map, int build)
     }
 
     Points = Vect_new_line_struct();
-    APoints = Vect_new_line_struct();
     Cats = Vect_new_cats_struct();
-    List = Vect_new_list();
-
+    
     if (plus->built < GV_BUILD_BASE) {
-	register int npoints, format, c;
-
-	format = G_info_format();
+	register int npoints, c;
 
 	/* 
 	 *  We shall go through all primitives in coor file and 

lib/vector/Vlib/build_ogr.c

@@ -128,11 +128,13 @@ static int add_line(struct Map_info *Map, int type, struct line_pnts *Points,
     plus = &(Map->plus);
 
     if (type != GV_CENTROID) {
-	offset = Map->fInfo.ogr.offset_num;	/* beginning in the offset array */
+	/* beginning in the offset array */
+	offset = Map->fInfo.ogr.offset_num;
     }
     else {
 	/* TODO : could be used to statore category ? */
-	offset = FID;		/* because centroids are read from topology, not from layer */
+	/* because centroids are read from topology, not from layer */
+	offset = FID; 
     }
     G_debug(4, "Register line: FID = %d offset = %ld", FID, offset);
     dig_line_box(Points, &box);

lib/vector/Vlib/read_ogr.c

@@ -354,7 +354,7 @@ int read_line(const struct Map_info *Map, OGRGeometryH hGeom, long offset,
 
     eType = wkbFlatten(OGR_G_GetGeometryType(hGeom));
     G_debug(4, "OGR geometry type: %d", eType);
-
+    
     switch (eType) {
     case wkbPoint:
 	G_debug(4, "\t->Point");
@@ -382,7 +382,7 @@ int read_line(const struct Map_info *Map, OGRGeometryH hGeom, long offset,
     case wkbMultiLineString:
     case wkbMultiPolygon:
     case wkbGeometryCollection:
-	G_debug(4, " \t->more geoms -> part %d", Map->fInfo.ogr.offset[offset]);
+	G_debug(4, "\t->more geoms -> part %d", Map->fInfo.ogr.offset[offset]);
 	hGeom2 = OGR_G_GetGeometryRef(hGeom, Map->fInfo.ogr.offset[offset]);
 	line = read_line(Map, hGeom2, offset + 1, Points);
 	if (eType == wkbPolygon || wkbMultiPolygon)
@@ -483,10 +483,13 @@ int V1_read_line_ogr(struct Map_info *Map,
     int type;
     OGRGeometryH hGeom;
 
+    struct Format_info_ogr *fInfo;
+    
+    fInfo = &(Map->fInfo.ogr);
     G_debug(4, "V1_read_line_ogr() offset = %lu offset_num = %lu",
-	    (long) offset, (long) Map->fInfo.ogr.offset_num);
+	    (long) offset, (long) fInfo->offset_num);
 
-    if (offset >= Map->fInfo.ogr.offset_num)
+    if (offset >= fInfo->offset_num)
 	return -2;
     
     if (line_p != NULL)
@@ -494,27 +497,27 @@ int V1_read_line_ogr(struct Map_info *Map,
     if (line_c != NULL)
 	Vect_reset_cats(line_c);
 
-    FID = Map->fInfo.ogr.offset[offset];
+    FID = fInfo->offset[offset];
     G_debug(4, "  FID = %ld", FID);
     
     /* coordinates */
     if (line_p != NULL) {
 	/* Read feature to cache if necessary */
-	if (Map->fInfo.ogr.feature_cache_id != FID) {
+	if (fInfo->feature_cache_id != FID) {
 	    G_debug(4, "Read feature (FID = %ld) to cache", FID);
-	    if (Map->fInfo.ogr.feature_cache) {
-		OGR_F_Destroy(Map->fInfo.ogr.feature_cache);
+	    if (fInfo->feature_cache) {
+		OGR_F_Destroy(fInfo->feature_cache);
 	    }
-	    Map->fInfo.ogr.feature_cache =
-		OGR_L_GetFeature(Map->fInfo.ogr.layer, FID);
-	    if (Map->fInfo.ogr.feature_cache == NULL) {
+	    fInfo->feature_cache =
+		OGR_L_GetFeature(fInfo->layer, FID);
+	    if (fInfo->feature_cache == NULL) {
 		G_fatal_error(_("Unable to get feature geometry, FID %ld"),
 			      FID);
 	    }
-	    Map->fInfo.ogr.feature_cache_id = FID;
+	    fInfo->feature_cache_id = FID;
 	}
 	
-	hGeom = OGR_F_GetGeometryRef(Map->fInfo.ogr.feature_cache);
+	hGeom = OGR_F_GetGeometryRef(fInfo->feature_cache);
 	if (hGeom == NULL) {
 	    G_fatal_error(_("Unable to get feature geometry, FID %ld"),
 			  FID);
@@ -560,30 +563,32 @@ int V2_read_line_ogr(struct Map_info *Map, struct line_pnts *line_p,
 	G_fatal_error(_("Attempt to read dead feature %d"), line);
 
     if (Line->type == GV_CENTROID) {
-	G_debug(4, "Centroid");
-	
 	if (line_p != NULL) {
 	    int i, found;
 	    struct bound_box box;
 	    struct boxlist list;
 	    struct P_topo_c *topo = (struct P_topo_c *)Line->topo;
+
+	    G_debug(4, "Centroid: area = %d", topo->area);
+	    Vect_reset_line(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);
-	    Vect_select_lines_by_box(Map, &box, Line->type, &list);
-	    
-	    found = 0;
-	    for (i = 0; i < list.n_values; i++) {
-		if (list.id[i] == line) {
-		    found = i;
-		    break;
+	    if (topo->area > 0 && topo->area <= Map->plus.n_areas) {
+		/* 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);
+		Vect_select_lines_by_box(Map, &box, Line->type, &list);
+		
+		found = 0;
+		for (i = 0; i < list.n_values; i++) {
+		    if (list.id[i] == line) {
+			found = i;
+			break;
+		    }
 		}
+		
+		Vect_append_point(line_p, list.box[found].E, list.box[found].N, 0.0);
 	    }
-
-	    Vect_reset_line(line_p);
-	    Vect_append_point(line_p, list.box[found].E, list.box[found].N, 0.0);
 	}
 
 	if (line_c != NULL) {

lib/vector/Vlib/write_ogr.c

@@ -35,15 +35,83 @@ static void V2__add_line_to_topo_ogr(struct Map_info *Map, int line,
 				     const struct line_pnts *points,
 				     const struct line_cats *cats)
 {
-    /* recycle code from build_ogr */
-    G_warning("feature not yet implemented, coming soon...");
+    int first, s, i;
+    int type, area, side, new_area[2];
 
+    struct Plus_head *plus;
+    struct P_line *Line;
+    
+    struct bound_box box, abox;
+    
+    G_debug(3, "V2__add_line_to_topo_ogr(): line = %d", line);
+
+    plus = &(Map->plus);
+    Line = plus->Line[line];
+    type = Line->type;
+
+    if (plus->built >= GV_BUILD_AREAS &&
+	type == GV_BOUNDARY) {	
+	struct P_topo_b *topo = (struct P_topo_b *)Line->topo;
+	
+	if (topo->N1 != topo->N2) {
+	    G_warning(_("Boundary is not closed. Skipping."));
+	    return;
+	}
+	
+	/* Build new areas/isles */
+	for (s = 0; s < 2; s++) {
+	    side = (s == 0 ? GV_LEFT : GV_RIGHT);
+	    area = Vect_build_line_area(Map, line, side);
+	    if (area > 0) {	/* area */
+		Vect_get_area_box(Map, area, &box);
+		if (first) {
+		    Vect_box_copy(&abox, &box);
+		    first = FALSE;
+		}
+		else
+		    Vect_box_extend(&abox, &box);
+	    }
+	    else if (area < 0) {
+		/* isle -> must be attached -> add to abox */
+		Vect_get_isle_box(Map, -area, &box);
+		if (first) {
+		    Vect_box_copy(&abox, &box);
+		    first = FALSE;
+		}
+		else
+		    Vect_box_extend(&abox, &box);
+	    }
+	    new_area[s] = area;
+	    G_debug(4, "Vect_build_line_area(): -> area = %d", area);
+	}
+
+	/* Attach centroid/isle to the new area */
+	if (plus->built >= GV_BUILD_ATTACH_ISLES)
+	    Vect_attach_isles(Map, &abox);
+	if (plus->built >= GV_BUILD_CENTROIDS)
+	    Vect_attach_centroids(Map, &abox);
+    }
+    
+    /* Add category index */
+    for (i = 0; i < cats->n_cats; i++) {
+	dig_cidx_add_cat_sorted(plus, cats->field[i], cats->cat[i], line,
+				type);
+    }
+    
     return;
 }
 
 /*!
   \brief Writes feature on level 1 (OGR interface)
 
+  Note:
+   - centroids are not supported in OGR, pseudotopo holds virtual
+     centroids
+   - boundaries are not supported in OGR, pseudotopo treats polygons
+     as boundaries
+     
+  \todo How to deal with OGRNullFID ?
+  
   \param Map pointer to Map_info structure
   \param type feature type
   \param points pointer to line_pnts structure (feature geometry) 
@@ -61,12 +129,15 @@ off_t V1_write_line_ogr(struct Map_info *Map, int type,
     struct field_info *Fi;
     struct Format_info_ogr *fInfo;
     
+    off_t offset;
+    
     OGRGeometryH       Ogr_geometry;
     OGRFeatureH        Ogr_feature;
     OGRFeatureDefnH    Ogr_featuredefn;
     OGRwkbGeometryType Ogr_geom_type;
 
     if (!Map->fInfo.ogr.layer) {
+	/* create OGR layer if doesn't exist */
 	if (V2_open_new_ogr(Map, type) < 0)
 	    return -1;
     }
@@ -91,11 +162,6 @@ off_t V1_write_line_ogr(struct Map_info *Map, int type,
     Ogr_geom_type = OGR_FD_GetGeomType(Ogr_featuredefn);
     
     /* determine matching OGR feature geometry type */
-    /* NOTE: centroids are not supported in OGR,
-     *       pseudotopo holds virtual centroids */
-    /* NOTE: boundaries are not supported in OGR,
-     *       pseudotopo treats polygons as boundaries */
-    
     if (type & (GV_POINT | GV_KERNEL)) {
 	if (Ogr_geom_type != wkbPoint &&
 	    Ogr_geom_type != wkbPoint25D) {
@@ -183,9 +249,15 @@ off_t V1_write_line_ogr(struct Map_info *Map, int type,
 					  fInfo->offset_alloc *
 					  sizeof(int));
     }
-    /* how to deal with OGRNullFID ? */
-    fInfo->offset[fInfo->offset_num] = (int)OGR_F_GetFID(Ogr_feature);
+
+    offset = fInfo->offset_num;
     
+    fInfo->offset[fInfo->offset_num++] = (int) OGR_F_GetFID(Ogr_feature);
+    if (Ogr_geom_type == wkbPolygon || Ogr_geom_type == wkbPolygon25D) {
+	/* register exterior ring in offset array */
+	fInfo->offset[fInfo->offset_num++] = 0; 
+    }
+      
     /* destroy */
     OGR_G_DestroyGeometry(Ogr_geometry);
     OGR_F_Destroy(Ogr_feature);
@@ -193,7 +265,9 @@ off_t V1_write_line_ogr(struct Map_info *Map, int type,
     if (ret != OGRERR_NONE)
 	return -1;
     
-    return fInfo->offset_num++;
+    G_debug(3, "V1_write_line_ogr(): -> offset = %d", offset);
+
+    return offset;
 }
 
 /*!
@@ -216,6 +290,7 @@ off_t V2_write_line_ogr(struct Map_info *Map, int type,
     struct bound_box box;
 
     line = 0;
+    plus = &(Map->plus);
     
     G_debug(3, "V2_write_line_ogr()");
     
@@ -224,17 +299,45 @@ off_t V2_write_line_ogr(struct Map_info *Map, int type,
 	return -1;
     
     /* Update topology */
-    plus = &(Map->plus);
-    /* Add line */
     if (plus->built >= GV_BUILD_BASE) {
 	dig_line_box(points, &box);
 	line = dig_add_line(plus, type, points, &box, offset);
-	G_debug(3, "  line added to topo with id = %d", line);
+	G_debug(3, "\tline added to topo with line = %d", line);
 	if (line == 1)
 	    Vect_box_copy(&(plus->box), &box);
 	else
 	    Vect_box_extend(&(plus->box), &box);
 
+	if (type == GV_BOUNDARY) {
+	    int ret, cline, lines[1];
+	    long FID;
+	    double x, y, area_size;
+	    
+	    struct bound_box box;
+	    struct line_pnts *CPoints;
+
+	    /* add virtual centroid to pseudo-topology */
+	    ret = Vect_get_point_in_poly(points, &x, &y);
+	    if (ret == 0) {
+		CPoints = Vect_new_line_struct();
+		Vect_append_point(CPoints, x, y, 0.0);
+		
+		FID = Map->fInfo.ogr.offset[offset];
+
+		dig_line_box(CPoints, &box);
+		cline = dig_add_line(plus, GV_CENTROID,
+				     CPoints, &box, FID);
+		G_debug(4, "\tCentroid: x = %f, y = %f, cat = %d, line = %d",
+			x, y, FID, cline);	  
+		dig_cidx_add_cat(plus, 1, (int) FID,
+				 cline, GV_CENTROID);
+		
+		Vect_destroy_line_struct(CPoints);
+	    }
+	    else {
+		G_warning(_("Unable to calculate centroid for area"));
+	    }
+	}
 	V2__add_line_to_topo_ogr(Map, line, points, cats);
     }