[GRASS-SVN] r37304 - grass/branches/develbranch_6/lib/vector/Vlib
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed May 20 06:20:04 EDT 2009
Author: martinl
Date: 2009-05-20 06:20:04 -0400 (Wed, 20 May 2009)
New Revision: 37304
Modified:
grass/branches/develbranch_6/lib/vector/Vlib/geos.c
Log:
vlib: implement Vect__read_line_geos()
Modified: grass/branches/develbranch_6/lib/vector/Vlib/geos.c
===================================================================
--- grass/branches/develbranch_6/lib/vector/Vlib/geos.c 2009-05-20 09:28:26 UTC (rev 37303)
+++ grass/branches/develbranch_6/lib/vector/Vlib/geos.c 2009-05-20 10:20:04 UTC (rev 37304)
@@ -1,16 +1,18 @@
/*!
- * \file geos.c
- *
- * \brief Vector library - GEOS support
- *
- * Higher level functions for reading/writing/manipulating vectors.
- *
- * (C) 2009 by the GRASS Development Team
- *
- * This program is free software under the GNU General Public License
- * (>=v2). Read the file COPYING that comes with GRASS for details.
- *
- * \author Martin Landa <landa.martin gmail.com>
+ \file geos.c
+
+ \brief Vector library - GEOS support
+
+ 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
+ (>=v2). Read the file COPYING that comes with GRASS for details.
+
+ \author Martin Landa <landa.martin gmail.com>
*/
#include <grass/config.h>
@@ -23,6 +25,8 @@
static struct line_pnts *Points;
+static GEOSGeometry *Vect__read_line_geos(struct Map_info *, long);
+
/*!
\brief Read vector feature and stores it as GEOSGeometry instance
@@ -36,25 +40,27 @@
*/
GEOSGeometry *Vect_read_line_geos(const struct Map_info *Map, int line)
{
- int type, dim;
- GEOSGeometry *geom;
- GEOSCoordSequence *pseq;
-
- if (!Points)
- Points = Vect_new_line_struct();
-
- G_debug(3, "Vect_read_line_geos(line=%d)", line);
-
- if (Vect_is_3d(Map))
- dim = 3;
- else
- dim = 2;
-
- type = Vect_read_line(Map, Points, NULL, line);
- if (type < 0)
- return NULL;
-
- return Vect_line_to_geos(Map, Points, type);
+ P_LINE *Line;
+
+ G_debug(3, "Vect_read_line_geos(): line = %d", line);
+
+ if (!VECT_OPEN(Map))
+ G_fatal_error("Vect_read_line_geos(): %s", _("vector map is not opened"));
+
+ if (line < 1 || line > Map->plus.n_lines)
+ G_fatal_error(_("Vect_read_line_geos(): feature id %d is not reasonable "
+ "(max features in vector map <%s>: %d)"),
+ line, Vect_get_full_name(Map), Map->plus.n_lines);
+
+ if (Map->format != GV_FORMAT_NATIVE)
+ G_fatal_error("Vect_read_line_geos(): %s", _("only native format supported"));
+
+ Line = Map->plus.Line[line];
+ if (Line == NULL)
+ G_fatal_error("Vect_read_line_geos(): %s %d",
+ _("Attempt to read dead line"), line);
+
+ return Vect__read_line_geos(Map, Line->offset);
}
/*!
@@ -70,7 +76,7 @@
*/
GEOSGeometry *Vect_read_area_geos(const struct Map_info * Map, int area)
{
- int i, type, dim, nholes, isle;
+ int i, type, nholes, isle;
GEOSGeometry *boundary, **holes;
if (!Points)
@@ -78,11 +84,6 @@
G_debug(3, "Vect_read_area_geos(area=%d)", area);
- if (Vect_is_3d(Map))
- dim = 3;
- else
- dim = 2;
-
if (Vect_get_area_points(Map, area, Points) == -1) {
G_fatal_error(_("Vect_read_area_geos(): unable to read area id %d"),
area);
@@ -127,35 +128,35 @@
GEOSGeometry *Vect_line_to_geos(const struct Map_info * Map,
const struct line_pnts * points, int type)
{
- int i, dim;
+ int i, with_z;
GEOSGeometry *geom;
GEOSCoordSequence *pseq;
- G_debug(3, "Vect_line_to_geos(type=%d)", type);
+ G_debug(3, "Vect_line_to_geos(): type = %d", type);
+
+ with_z = Vect_is_3d(Map);
+
+ /* read only points / lines / boundaries */
+ if (!(type & (GV_POINT | GV_LINES)))
+ return NULL;
- if (Vect_is_3d(Map))
- dim = 3;
- else
- dim = 2;
-
- if (type == GV_POINT) {
+ if (type == GV_POINT) {
if (points->n_points != 1)
+ /* point is not valid */
return NULL;
}
- else { /* GV_LINES */
+ else {
if (points->n_points < 2)
+ /* line/boundary is not valid */
return NULL;
}
-
- if (!(type & (GV_POINT | GV_LINES)))
- return NULL;
-
- pseq = GEOSCoordSeq_create(points->n_points, dim);
-
+
+ pseq = GEOSCoordSeq_create(points->n_points, with_z ? 3 : 2);
+
for (i = 0; i < points->n_points; i++) {
GEOSCoordSeq_setX(pseq, i, points->x[i]);
GEOSCoordSeq_setY(pseq, i, points->y[i]);
- if (dim == 3)
+ if (with_z)
GEOSCoordSeq_setZ(pseq, i, points->z[i]);
}
@@ -163,16 +164,134 @@
geom = GEOSGeom_createPoint(pseq);
else if (type == GV_LINE)
geom = GEOSGeom_createLineString(pseq);
- else { /* GV_BOUNDARY */
+ else { /* boundary */
geom = GEOSGeom_createLineString(pseq);
if (GEOSisRing(geom)) {
/* GEOSGeom_destroy(geom); */
geom = GEOSGeom_createLinearRing(pseq);
}
}
-
+
/* GEOSCoordSeq_destroy(pseq); */
return geom;
}
+
+/*!
+ \brief Read line from coor file
+
+ \param Map pointer to Map_info
+ \param offset lineoffset
+
+ \return pointer to GEOSGeometry
+ \return NULL on error
+ \return NULL dead line
+ \return NULL end of file
+*/
+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);
+
+ Map->head.last_offset = offset;
+
+ /* reads must set in_head, but writes use default */
+ dig_set_cur_port(&(Map->head.port));
+
+ dig_fseek(&(Map->dig_fp), offset, 0);
+
+ if (0 >= dig__fread_port_C(&rhead, 1, &(Map->dig_fp)))
+ return NULL; /* end of file */
+
+ if (!(rhead & 0x01)) /* dead line */
+ return NULL;
+
+ rhead >>= 2;
+ type = dig_type_from_store((int)rhead);
+
+ /* read only points / lines / boundaries */
+ if (!(type & (GV_POINT | GV_LINES)))
+ return NULL;
+
+ if (type & GV_POINTS) {
+ n_points = 1;
+ }
+ 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);
+
+ pseq = GEOSCoordSeq_create(n_points, (Map->head.with_z) ? 3 : 2);
+
+ x = (double *) G_malloc(n_points * sizeof(double));
+ y = (double *) G_malloc(n_points * sizeof(double));
+ if (Map->head.with_z)
+ z = (double *) G_malloc(n_points * sizeof(double));
+ else
+ z = NULL;
+
+ if (0 >= dig__fread_port_D(x, n_points, &(Map->dig_fp)))
+ return NULL; /* end of file */
+
+ if (0 >= dig__fread_port_D(y, n_points, &(Map->dig_fp)))
+ return NULL; /* end of file */
+
+ if (Map->head.with_z) {
+ if (0 >= dig__fread_port_D(z, n_points, &(Map->dig_fp)))
+ return NULL; /* end of file */
+
+ }
+
+ for (i = 0; i < n_points; 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)));
+
+ 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;
+}
+
#endif /* HAVE_GEOS */
More information about the grass-commit
mailing list