[GRASS-SVN] r49178 - grass/trunk/lib/vector/Vlib
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Nov 11 08:24:52 EST 2011
Author: martinl
Date: 2011-11-11 05:24:52 -0800 (Fri, 11 Nov 2011)
New Revision: 49178
Modified:
grass/trunk/lib/vector/Vlib/build_nat.c
grass/trunk/lib/vector/Vlib/build_ogr.c
grass/trunk/lib/vector/Vlib/read_ogr.c
grass/trunk/lib/vector/Vlib/write_ogr.c
Log:
vlib: implement V2__add_line_to_topo_ogr()
Modified: grass/trunk/lib/vector/Vlib/build_nat.c
===================================================================
--- grass/trunk/lib/vector/Vlib/build_nat.c 2011-11-11 09:08:11 UTC (rev 49177)
+++ grass/trunk/lib/vector/Vlib/build_nat.c 2011-11-11 13:24:52 UTC (rev 49178)
@@ -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 @@
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_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 @@
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 @@
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 @@
\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 @@
/*!
\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 @@
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 @@
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 @@
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 @@
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 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,15 +532,11 @@
}
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;
+ register int npoints, c;
- format = G_info_format();
-
/*
* We shall go through all primitives in coor file and
* add new node for each end point to nodes structure
Modified: grass/trunk/lib/vector/Vlib/build_ogr.c
===================================================================
--- grass/trunk/lib/vector/Vlib/build_ogr.c 2011-11-11 09:08:11 UTC (rev 49177)
+++ grass/trunk/lib/vector/Vlib/build_ogr.c 2011-11-11 13:24:52 UTC (rev 49178)
@@ -128,11 +128,13 @@
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);
Modified: grass/trunk/lib/vector/Vlib/read_ogr.c
===================================================================
--- grass/trunk/lib/vector/Vlib/read_ogr.c 2011-11-11 09:08:11 UTC (rev 49177)
+++ grass/trunk/lib/vector/Vlib/read_ogr.c 2011-11-11 13:24:52 UTC (rev 49178)
@@ -354,7 +354,7 @@
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 @@
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 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 @@
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 @@
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) {
Modified: grass/trunk/lib/vector/Vlib/write_ogr.c
===================================================================
--- grass/trunk/lib/vector/Vlib/write_ogr.c 2011-11-11 09:08:11 UTC (rev 49177)
+++ grass/trunk/lib/vector/Vlib/write_ogr.c 2011-11-11 13:24:52 UTC (rev 49178)
@@ -35,15 +35,83 @@
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 @@
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 @@
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 @@
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 @@
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 @@
struct bound_box box;
line = 0;
+ plus = &(Map->plus);
G_debug(3, "V2_write_line_ogr()");
@@ -224,17 +299,45 @@
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);
}
More information about the grass-commit
mailing list