progress in geos support (new Vlib fns)

	 (merge from devbr6, https://trac.osgeo.org/grass/changeset/37316)


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

include/Vect.h

@@ -470,9 +470,11 @@ int Vect_attach_centroids(struct Map_info *, const BOUND_BOX *);
 
     /* 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 */

lib/vector/Vlib/geos.c

@@ -1,16 +1,16 @@
 /*!
- * \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.
+  
+  (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>
@@ -21,44 +21,56 @@
 #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)
 {
-    int type, dim;
+    P_LINE *Line;
     
-    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 < 1) /* ignore dead lines */
-	return NULL;
-
-    return Vect_line_to_geos(Map, Points, type);
+    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);
 }
 
 /*!
-   \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 
@@ -66,26 +78,14 @@ GEOSGeometry *Vect_read_line_geos(const struct Map_info *Map, int line)
    \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, dim, 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_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);
-    }
-    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);
@@ -97,12 +97,12 @@ GEOSGeometry *Vect_read_area_geos(const struct Map_info * Map, int area)
 	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);
 }
 
@@ -112,9 +112,9 @@ GEOSGeometry *Vect_read_area_geos(const struct Map_info * Map, int area)
    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)
@@ -122,38 +122,38 @@ GEOSGeometry *Vect_read_area_geos(const struct Map_info * Map, int area)
    \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, dim;
+    int i, with_z;
     GEOSGeometry *geom;
     GEOSCoordSequence *pseq;
 
-    G_debug(3, "Vect_line_to_geos(type=%d)", type);
-
-    if (Vect_is_3d(Map))
-	dim = 3;
-    else
-	dim = 2;
+    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 (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]);
     }
 
@@ -161,16 +161,355 @@ GEOSGeometry *Vect_line_to_geos(const struct Map_info * Map,
 	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
+  
+  You should free allocated memory by GEOSGeom_destroy().
+
+  \param Map pointer to Map_info
+  \param offset line offset
+ 
+  \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 type;
+    double *x, *y, *z;
+    
+    GEOSGeometry *geom;
+    GEOSCoordSequence *pseq;
+    
+    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 */
+    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 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 GEOSCoordSeq_create(0, (Map->head.with_z) ? 3 : 2);
+    
+    if (type & GV_POINTS) {
+	n_points = 1;
+    }
+    else {
+	if (0 >= dig__fread_port_I(&n_points, 1, &(Map->dig_fp)))
+	    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)));
+    
+    G_free((void *) x);
+    G_free((void *) y);
+    if (z)
+	G_free((void *) z);
+    
+    return 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 */
+    }
+    
+    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;
+    }
+
+    /* 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);
+	    }
+	}
+    }
+    
+    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];
+    
+    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);
+	    }
+	}
+    }
+    
+    G_free((void *) pseq);
+
+    return pseq_shell;
+}
+
 #endif /* HAVE_GEOS */