[GRASS-SVN] r37316 - in grass/branches/develbranch_6: include
lib/vector/Vlib
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed May 20 17:56:34 EDT 2009
Author: martinl
Date: 2009-05-20 17:56:34 -0400 (Wed, 20 May 2009)
New Revision: 37316
Modified:
grass/branches/develbranch_6/include/Vect.h
grass/branches/develbranch_6/lib/vector/Vlib/geos.c
Log:
progress in geos support (new Vlib fns)
Modified: grass/branches/develbranch_6/include/Vect.h
===================================================================
--- grass/branches/develbranch_6/include/Vect.h 2009-05-20 21:24:51 UTC (rev 37315)
+++ grass/branches/develbranch_6/include/Vect.h 2009-05-20 21:56:34 UTC (rev 37316)
@@ -471,9 +471,11 @@
/* GEOS support */
#ifdef HAVE_GEOS
-GEOSGeometry *Vect_read_line_geos(const struct Map_info *, int);
-GEOSGeometry *Vect_read_area_geos(const struct Map_info *, int);
-GEOSGeometry *Vect_line_to_geos(const struct Map_info *, const struct line_pnts*, int);
+GEOSGeometry *Vect_read_line_geos(struct Map_info *, int);
+GEOSGeometry *Vect_line_to_geos(struct Map_info *, const struct line_pnts*, int);
+GEOSGeometry *Vect_read_area_geos(struct Map_info *, int);
+GEOSCoordSequence *Vect_get_area_points_geos(struct Map_info *, int);
+GEOSCoordSequence *Vect_get_isle_points_geos(struct Map_info *, int);
#endif
#endif /* GRASS_VECT_H */
Modified: grass/branches/develbranch_6/lib/vector/Vlib/geos.c
===================================================================
--- grass/branches/develbranch_6/lib/vector/Vlib/geos.c 2009-05-20 21:24:51 UTC (rev 37315)
+++ grass/branches/develbranch_6/lib/vector/Vlib/geos.c 2009-05-20 21:56:34 UTC (rev 37316)
@@ -5,8 +5,6 @@
Higher level functions for reading/writing/manipulating vectors.
- Note: current GEOS support is *very* slow
-
(C) 2009 by the GRASS Development Team
This program is free software under the GNU General Public License
@@ -23,22 +21,28 @@
#ifdef HAVE_GEOS
#include <geos_c.h>
-static struct line_pnts *Points;
-
static GEOSGeometry *Vect__read_line_geos(struct Map_info *, long);
+static GEOSCoordSequence *V1_read_line_geos(struct Map_info *, long);
+static GEOSCoordSequence *V2_read_line_geos(struct Map_info *, int);
/*!
\brief Read vector feature and stores it as GEOSGeometry instance
- Note: Free allocated memory by GEOSGeom_destroy().
+ Supported feature types:
+ - GV_POINT -> POINT
+ - GV_LINE -> LINESTRING
+ - GV_BOUNDARY -> LINESTRING / LINEARRING
+
+ You should free allocated memory by GEOSGeom_destroy().
\param Map pointer to Map_info structure
\param line feature id
\return pointer to GEOSGeometry instance
+ \return empty GEOSGeometry for unsupported feature type
\return NULL on error
*/
-GEOSGeometry *Vect_read_line_geos(const struct Map_info *Map, int line)
+GEOSGeometry *Vect_read_line_geos(struct Map_info *Map, int line)
{
P_LINE *Line;
@@ -64,9 +68,9 @@
}
/*!
- \brief Read vector area and stores it as GEOSGeometry instance
+ \brief Read vector area and stores it as GEOSGeometry instance (polygon)
- Note: Free allocated memory by GEOSGeom_destroy().
+ You should free allocated memory by GEOSGeom_destroy().
\param Map pointer to Map_info structure
\param area area id
@@ -74,21 +78,14 @@
\return pointer to GEOSGeometry instance
\return NULL on error
*/
-GEOSGeometry *Vect_read_area_geos(const struct Map_info * Map, int area)
+GEOSGeometry *Vect_read_area_geos(struct Map_info * Map, int area)
{
- int i, type, nholes, isle;
+ int i, nholes, isle;
GEOSGeometry *boundary, **holes;
+
+ G_debug(3, "Vect_read_area_geos(): area = %d", area);
- if (!Points)
- Points = Vect_new_line_struct();
-
- G_debug(3, "Vect_read_area_geos(area=%d)", area);
-
- if (Vect_get_area_points(Map, area, Points) == -1) {
- G_fatal_error(_("Vect_read_area_geos(): unable to read area id %d"),
- area);
- }
- boundary = Vect_line_to_geos(Map, Points, GV_BOUNDARY);
+ boundary = GEOSGeom_createLinearRing(Vect_get_area_points_geos(Map, area));
if (!boundary) {
G_fatal_error(_("Vect_read_area_geos(): unable to read area id %d"),
area);
@@ -100,12 +97,12 @@
isle = Vect_get_area_isle(Map, area, i);
if (isle < 1)
continue;
- if (Vect_get_isle_points(Map, isle, Points) < 0)
+ holes[i] = GEOSGeom_createLinearRing(Vect_get_isle_points_geos(Map, isle));
+ if (!(holes[i]))
G_fatal_error(_("Vect_read_area_geos(): unable to read isle id %d of area id %d"),
isle, area);
- holes[i] = Vect_line_to_geos(Map, Points, GV_BOUNDARY);
}
-
+
return GEOSGeom_createPolygon(boundary, holes, nholes);
}
@@ -115,9 +112,9 @@
Supported types:
- GV_POINT -> POINT
- GV_LINE -> LINESTRING
- - GV_BOUNDARY -> LINERING
+ - GV_BOUNDARY -> LINEARRING
- Note: Free allocated memory by GEOSGeom_destroy().
+ You should free allocated memory by GEOSGeom_destroy().
\param Map pointer to Map_info structure
\param type feature type (see supported types)
@@ -125,7 +122,7 @@
\return pointer to GEOSGeometry instance
\return NULL on error
*/
-GEOSGeometry *Vect_line_to_geos(const struct Map_info * Map,
+GEOSGeometry *Vect_line_to_geos(struct Map_info * Map,
const struct line_pnts * points, int type)
{
int i, with_z;
@@ -180,8 +177,10 @@
/*!
\brief Read line from coor file
+ You should free allocated memory by GEOSGeom_destroy().
+
\param Map pointer to Map_info
- \param offset lineoffset
+ \param offset line offset
\return pointer to GEOSGeometry
\return NULL on error
@@ -190,17 +189,99 @@
*/
GEOSGeometry *Vect__read_line_geos(struct Map_info *Map, long offset)
{
- int i;
- int n_points;
int type;
- char rhead;
double *x, *y, *z;
GEOSGeometry *geom;
GEOSCoordSequence *pseq;
- G_debug(3, "Vect__read_line_geos(): offset = %ld", offset);
+ pseq = V1_read_line_geos(Map, offset);
+ if (!pseq)
+ G_fatal_error(_("Unable to read line offset %ld"), offset);
+ if (type & GV_POINT) {
+ G_debug(3, " geos_type = point");
+ geom = GEOSGeom_createPoint(pseq);
+ }
+ else if (type & GV_LINE) {
+ G_debug(3, " geos_type = linestring");
+ geom = GEOSGeom_createLineString(pseq);
+ }
+ else { /* boundary */
+ geom = GEOSGeom_createLineString(pseq);
+ if (GEOSisRing(geom)) {
+ /* GEOSGeom_destroy(geom); */
+ geom = GEOSGeom_createLinearRing(pseq);
+ G_debug(3, " geos_type = linearring");
+ }
+ else {
+ G_debug(3, " geos_type = linestring");
+ }
+ }
+
+ G_free((void *) x);
+ G_free((void *) y);
+ if (z)
+ G_free((void *) z);
+
+ /* GEOSCoordSeq_destroy(pseq); */
+
+ return geom;
+}
+
+/*!
+ \brief Read line from coor file into GEOSCoordSequence
+
+ You should free allocated memory by GEOSCoordSeq_destroy().
+
+ \param Map pointer to Map_info
+ \param line line id
+
+ \return pointer to GEOSCoordSequence
+ \return empty GEOSCoordSequence for dead line or unsuppored feature type
+ \return NULL end of file
+*/
+GEOSCoordSequence *V2_read_line_geos(struct Map_info *Map, int line)
+{
+ P_LINE *Line;
+
+ G_debug(3, "V2_read_line_geos(): line = %d", line);
+
+ Line = Map->plus.Line[line];
+
+ if (Line == NULL)
+ G_fatal_error("V2_read_line_geos(): %s %d",
+ _("Attempt to read dead line"), line);
+
+ return V1_read_line_geos(Map, Line->offset);
+}
+
+
+/*!
+ \brief Read feature from coor file into GEOSCoordSequence
+
+ Note: Function reads only points, lines and boundaries, other
+ feature types are ignored (empty coord array is returned)!
+
+ You should free allocated memory by GEOSCoordSeq_destroy().
+
+ \param Map pointer to Map_info
+ \param offset line offset
+
+ \return pointer to GEOSCoordSequence
+ \return empty GEOSCoordSequence for dead line or unsuppored feature type
+ \return NULL end of file
+*/
+GEOSCoordSequence *V1_read_line_geos(struct Map_info *Map, long offset)
+{
+ int i, n_points, type;
+ char rhead;
+ double *x, *y, *z;
+
+ GEOSCoordSequence *pseq;
+
+ G_debug(3, "V1_read_line_geos(): offset = %ld", offset);
+
Map->head.last_offset = offset;
/* reads must set in_head, but writes use default */
@@ -212,14 +293,14 @@
return NULL; /* end of file */
if (!(rhead & 0x01)) /* dead line */
- return NULL;
+ return GEOSCoordSeq_create(0, (Map->head.with_z) ? 3 : 2);
rhead >>= 2;
type = dig_type_from_store((int)rhead);
/* read only points / lines / boundaries */
if (!(type & (GV_POINT | GV_LINES)))
- return NULL;
+ return GEOSCoordSeq_create(0, (Map->head.with_z) ? 3 : 2);
if (type & GV_POINTS) {
n_points = 1;
@@ -227,9 +308,6 @@
else {
if (0 >= dig__fread_port_I(&n_points, 1, &(Map->dig_fp)))
return NULL;
- if (n_points < 2)
- /* line / boundary is not valid */
- return NULL;
}
G_debug(3, " n_points = %d dim = %d", n_points, (Map->head.with_z) ? 3 : 2);
@@ -256,42 +334,182 @@
}
for (i = 0; i < n_points; i++) {
- GEOSCoordSeq_setX(pseq, i, x[i]);
+ GEOSCoordSeq_setX(pseq, i, x[i]);
GEOSCoordSeq_setY(pseq, i, y[i]);
if (Map->head.with_z)
GEOSCoordSeq_setZ(pseq, i, z[i]);
}
G_debug(3, " off = %ld", dig_ftell(&(Map->dig_fp)));
+
+ G_free((void *) x);
+ G_free((void *) y);
+ if (z)
+ G_free((void *) z);
+
+ return pseq;
+}
- if (type & GV_POINT) {
- G_debug(3, " geos_type = point");
- geom = GEOSGeom_createPoint(pseq);
+/*!
+ \brief Returns the polygon array of points, i.e. outer ring (shell)
+
+ You should free allocated memory by GEOSCoordSeq_destroy().
+
+ See also Vect_get_area_points().
+
+ \param Map pointer to Map_info
+ \param area area id
+
+ \return pointer to GEOSCoordSequence
+ \return empty GEOSCoordSequence for dead area
+ \return NULL on error
+ */
+GEOSCoordSequence *Vect_get_area_points_geos(struct Map_info *Map, int area)
+{
+ int i, j, k;
+ int line, aline;
+ unsigned int n_points, n_points_shell;
+ double x, y, z;
+
+ struct Plus_head *Plus;
+ P_AREA *Area;
+
+ GEOSCoordSequence **pseq, *pseq_shell;
+
+ G_debug(3, "Vect_get_area_points_geos(): area = %d", area);
+
+ Plus = &(Map->plus);
+ Area = Plus->Area[area];
+
+ if (Area == NULL) { /* dead area */
+ G_warning(_("Attempt to read points of nonexistent area id %d"), area);
+ return NULL; /* error , because we should not read dead areas */
}
- else if (type & GV_LINE) {
- G_debug(3, " geos_type = linestring");
- geom = GEOSGeom_createLineString(pseq);
+
+ G_debug(3, " n_lines = %d", Area->n_lines);
+ pseq = (GEOSCoordSequence **) G_malloc(Area->n_lines * sizeof(GEOSCoordSequence *));
+
+ n_points_shell = 0;
+ for (i = 0; i < Area->n_lines; i++) {
+ line = Area->lines[i];
+ aline = abs(line);
+ G_debug(3, " append line(%d) = %d", i, line);
+
+ /*
+ TODO
+ if (line > 0)
+ dir = GV_FORWARD;
+ else
+ dir = GV_BACKWARD;
+ */
+
+ pseq[i] = V2_read_line_geos(Map, aline);
+ if (!(pseq[i])) {
+ G_fatal_error(_("Unable to read feature id %d"), aline);
+ }
+
+ GEOSCoordSeq_getSize(pseq[i], &n_points);
+ G_debug(3, " line n_points = %d", n_points);
+ n_points_shell += n_points;
}
- else { /* boundary */
- geom = GEOSGeom_createLineString(pseq);
- if (GEOSisRing(geom)) {
- /* GEOSGeom_destroy(geom); */
- geom = GEOSGeom_createLinearRing(pseq);
- G_debug(3, " geos_type = linearring");
+
+ /* create shell (outer ring) */
+ pseq_shell = GEOSCoordSeq_create(n_points_shell, Map->head.with_z ? 3 : 2);
+ k = 0;
+ for (i = 0; i < Area->n_lines; i++) {
+ GEOSCoordSeq_getSize(pseq[i], &n_points);
+ for (j = 0; j < (int) n_points; j++, k++) {
+ GEOSCoordSeq_getX(pseq[i], j, &x);
+ GEOSCoordSeq_setX(pseq_shell, k, x);
+ GEOSCoordSeq_getY(pseq[i], j, &y);
+ GEOSCoordSeq_setY(pseq_shell, k, y);
+ if (Map->head.with_z) {
+ GEOSCoordSeq_getY(pseq[i], j, &z);
+ GEOSCoordSeq_setZ(pseq_shell, k, z);
+ }
}
- else {
- G_debug(3, " geos_type = linestring");
- }
}
- G_free((void *) x);
- G_free((void *) y);
- if (z)
- G_free((void *) z);
+ G_free((void *) pseq);
+
+ return pseq_shell;
+}
+
+/*!
+ \brief Returns the polygon (isle) array of points (inner ring)
+
+ You should free allocated memory by GEOSCoordSeq_destroy().
+
+ See also Vect_get_isle_points().
+
+ \param Map pointer to Map_info
+ \param isle isel id
+
+ \return pointer to GEOSGeometry
+ \return NULL on error or dead line
+ */
+GEOSCoordSequence *Vect_get_isle_points_geos(struct Map_info *Map, int isle)
+{
+ int i, j, k;
+ int line, aline;
+ unsigned n_points, n_points_shell;
+ double x, y, z;
+
+ struct Plus_head *Plus;
+ P_ISLE *Isle;
+
+ GEOSCoordSequence **pseq, *pseq_shell;
+
+ G_debug(3, "Vect_get_isle_points_geos(): isle = %d", isle);
+
+ Plus = &(Map->plus);
+ Isle = Plus->Isle[isle];
- /* GEOSCoordSeq_destroy(pseq); */
+ G_debug(3, " n_lines = %d", Isle->n_lines);
+ for (i = 0; i < Isle->n_lines; i++) {
+ line = Isle->lines[i];
+ aline = abs(line);
+ G_debug(3, " append line(%d) = %d", i, line);
+
+
+ pseq[i] = V2_read_line_geos(Map, aline);
+ if (!(pseq[i])) {
+ G_fatal_error(_("Unable to read feature id %d"), aline);
+ }
+
+ GEOSCoordSeq_getSize(pseq[i], &n_points);
+ G_debug(3, " line n_points = %d", n_points);
+ n_points_shell += n_points;
+
+ /*
+ TODO
+ if (line > 0)
+ dir = GV_FORWARD;
+ else
+ dir = GV_BACKWARD;
+ */
+ }
+
+ /* create shell (outer ring) */
+ pseq_shell = GEOSCoordSeq_create(n_points_shell, Map->head.with_z ? 3 : 2);
+ k = 0;
+ for (i = 0; i < Isle->n_lines; i++) {
+ GEOSCoordSeq_getSize(pseq[i], &n_points);
+ for (j = 0; j < (int) n_points; j++, k++) {
+ GEOSCoordSeq_getX(pseq[i], j, &x);
+ GEOSCoordSeq_setX(pseq_shell, k, x);
+ GEOSCoordSeq_getY(pseq[i], j, &y);
+ GEOSCoordSeq_setY(pseq_shell, k, y);
+ if (Map->head.with_z) {
+ GEOSCoordSeq_getY(pseq[i], j, &z);
+ GEOSCoordSeq_setZ(pseq_shell, k, z);
+ }
+ }
+ }
- return geom;
+ G_free((void *) pseq);
+
+ return pseq_shell;
}
#endif /* HAVE_GEOS */
More information about the grass-commit
mailing list