vector library: Added new methods to generate Well Known Text (WTK) and Well Known Binary (WTB) out of GEOS geometries

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

include/defs/vector.h

@@ -313,7 +313,7 @@ int Vect_select_lines_by_box(struct Map_info *, const struct bound_box *,
 int Vect_select_areas_by_box(struct Map_info *, const struct bound_box *,
                              struct boxlist *);
 int Vect_select_isles_by_box(struct Map_info *, const struct bound_box *,
-			     struct boxlist *);
+                 struct boxlist *);
 int Vect_select_nodes_by_box(struct Map_info *, const struct bound_box *,
                              struct ilist *);
 int Vect_find_node(struct Map_info *, double, double, double, double, int);
@@ -595,10 +595,15 @@ int Vect_attach_centroids(struct Map_info *, const struct bound_box *);
     /* GEOS support */
 #ifdef HAVE_GEOS
 GEOSGeometry *Vect_read_line_geos(struct Map_info *, int, int*);
-GEOSGeometry *Vect_line_to_geos(struct Map_info *, const struct line_pnts*, int);
+GEOSGeometry *Vect_line_to_geos(const struct line_pnts*, int, 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);
+char *Vect_line_to_wkt(const struct line_pnts *, int, int);
+unsigned char *Vect_line_to_wkb(const struct line_pnts *,
+                                int, int, size_t *);
+char *Vect_read_area_to_wkt(struct Map_info *, int);
+unsigned char *Vect_read_area_to_wkb(struct Map_info *, int, size_t *);
 #endif
 
     /* Raster color tables */

lib/vector/Vlib/geos.c

@@ -1,15 +1,15 @@
 /*!
   \file lib/vector/Vlib/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>
  */
 
@@ -31,13 +31,13 @@ static GEOSCoordSequence *read_polygon_points(struct Map_info *, int, int*);
     - 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
    \param[out] type feature type or NULL
-   
+
    \return pointer to GEOSGeometry instance
    \return empty GEOSGeometry for unsupported feature type
    \return NULL on error
@@ -45,25 +45,25 @@ static GEOSCoordSequence *read_polygon_points(struct Map_info *, int, int*);
 GEOSGeometry *Vect_read_line_geos(struct Map_info *Map, int line, int *type)
 {
     struct 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"));
-    
+        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);
-    
+        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"));
-    
+        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);
-    
+        G_fatal_error("Vect_read_line_geos(): %s %d",
+                      _("Attempt to read dead line"), line);
+
     return Vect__read_line_geos(Map, Line->offset, type);
 }
 
@@ -73,7 +73,7 @@ GEOSGeometry *Vect_read_line_geos(struct Map_info *Map, int line, int *type)
    You should free allocated memory by GEOSGeom_destroy().
 
    \param Map pointer to Map_info structure
-   \param area area id 
+   \param area area id
 
    \return pointer to GEOSGeometry instance
    \return NULL on error
@@ -82,29 +82,29 @@ GEOSGeometry *Vect_read_area_geos(struct Map_info * Map, int area)
 {
     int i, nholes, isle;
     GEOSGeometry *boundary, *poly, **holes;
-    
+
     G_debug(3, "Vect_read_area_geos(): area = %d", area);
 
     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);
+        G_fatal_error(_("Vect_read_area_geos(): unable to read area id %d"),
+                      area);
     }
 
     nholes = Vect_get_area_num_isles(Map, area);
     holes = (GEOSGeometry **) G_malloc(nholes * sizeof(GEOSGeometry *));
     for (i = 0; i < nholes; i++) {
-	isle = Vect_get_area_isle(Map, area, i);
-	if (isle < 1) {
-	    nholes--;
-	    continue;
-	}
-	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);
+        isle = Vect_get_area_isle(Map, area, i);
+        if (isle < 1) {
+            nholes--;
+            continue;
+        }
+        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);
     }
-    
+
     poly = GEOSGeom_createPolygon(boundary, holes, nholes);
     G_free(holes);
 
@@ -116,73 +116,72 @@ GEOSGeometry *Vect_read_area_geos(struct Map_info * Map, int area)
 
    Supported types:
    - GV_POINT    -> POINT
+   - GV_CENTROID -> POINT
    - GV_LINE     -> LINESTRING
    - GV_BOUNDARY -> LINEARRING
 
    You should free allocated memory by GEOSGeom_destroy().
 
-   \param Map pointer to Map_info structure
    \param points pointer to line_pnts structure
    \param type feature type (see supported types)
+   \param with_z Set to 1 if the feature is 3d, 0 otherwise
 
    \return pointer to GEOSGeometry instance
    \return NULL on error
  */
-GEOSGeometry *Vect_line_to_geos(struct Map_info *Map,
-				const struct line_pnts *points, int type)
+GEOSGeometry *Vect_line_to_geos(const struct line_pnts *points,
+                                int type, int with_z)
 {
-    int i, with_z;
+    int i;
     GEOSGeometry *geom;
     GEOSCoordSequence *pseq;
 
     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 | GV_CENTROID | GV_LINES)))
+        return NULL;
 
-    if (type == GV_POINT) { 
-	if (points->n_points != 1)
-	    /* point is not valid */
-	    return NULL;
+    if (type == GV_POINT || type == GV_CENTROID) {
+        if (points->n_points != 1)
+            /* point is not valid */
+            return NULL;
     }
-    else {			
-	if (points->n_points < 2)
-	    /* line/boundary is not valid */
-	    return NULL;
+    else {
+        if (points->n_points < 2)
+            /* line/boundary is not valid */
+            return NULL;
     }
-    
+
     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 (with_z)
-	    GEOSCoordSeq_setZ(pseq, i, points->z[i]);
+        GEOSCoordSeq_setX(pseq, i, points->x[i]);
+        GEOSCoordSeq_setY(pseq, i, points->y[i]);
+        if (with_z)
+            GEOSCoordSeq_setZ(pseq, i, points->z[i]);
     }
 
-    if (type == GV_POINT)
-	geom = GEOSGeom_createPoint(pseq);
+    if (type == GV_POINT || type == GV_CENTROID)
+        geom = GEOSGeom_createPoint(pseq);
     else if (type == GV_LINE)
-	geom = GEOSGeom_createLineString(pseq);
+        geom = GEOSGeom_createLineString(pseq);
     else { /* boundary */
-	geom = GEOSGeom_createLineString(pseq);
-	if (GEOSisRing(geom)) {
-	    /* GEOSGeom_destroy(geom); */
-	    geom = GEOSGeom_createLinearRing(pseq);
-	}
+        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
@@ -197,50 +196,50 @@ GEOSGeometry *Vect_line_to_geos(struct Map_info *Map,
 GEOSGeometry *Vect__read_line_geos(struct Map_info *Map, long offset, int *type)
 {
     int ftype;
-    
+
     GEOSGeometry *geom;
     GEOSCoordSequence *pseq;
-    
+
     pseq = V1_read_line_geos(Map, offset, &ftype);
     if (!pseq)
-	G_fatal_error(_("Unable to read line offset %ld"), offset);
-    
+        G_fatal_error(_("Unable to read line offset %ld"), offset);
+
     if (ftype & GV_POINT) {
-	G_debug(3, "    geos_type = point");
-	geom = GEOSGeom_createPoint(pseq);
+        G_debug(3, "    geos_type = point");
+        geom = GEOSGeom_createPoint(pseq);
     }
     else if (ftype & GV_LINE) {
-	G_debug(3, "    geos_type = linestring");
-	geom = GEOSGeom_createLineString(pseq);
+        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");
-	}
+        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");
+        }
     }
-        
+
     /* GEOSCoordSeq_destroy(pseq); */
-    
+
     if (type)
       *type = ftype;
-    
+
     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
@@ -249,31 +248,31 @@ GEOSCoordSequence *V2_read_line_geos(struct Map_info *Map, int line)
 {
     int ftype;
     struct 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);
-    
+        G_fatal_error("V2_read_line_geos(): %s %d",
+                      _("Attempt to read dead line"), line);
+
     return V1_read_line_geos(Map, Line->offset, &ftype);
 }
 
 
-/*!  
+/*!
   \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
   \param[out] type feature type
-  
+
   \return pointer to GEOSCoordSequence
   \return empty GEOSCoordSequence for dead line or unsuppored feature type
   \return NULL end of file
@@ -285,103 +284,103 @@ GEOSCoordSequence *V1_read_line_geos(struct Map_info *Map, long offset, int *typ
     char rhead, nc;
     long size;
     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);
-
-    if (rhead & 0x02)		/* categories exists */
-	do_cats = 1;		/* do not return here let file offset moves forward to next */
-    else			/* line */
-	do_cats = 0;
-    
+        return NULL;            /* end of file */
+
+    if (!(rhead & 0x01))        /* dead line */
+        return GEOSCoordSeq_create(0, (Map->head.with_z) ? 3 : 2);
+
+    if (rhead & 0x02)           /* categories exists */
+        do_cats = 1;            /* do not return here let file offset moves forward to next */
+    else                        /* line */
+        do_cats = 0;
+
     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);
- 
+        return GEOSCoordSeq_create(0, (Map->head.with_z) ? 3 : 2);
+
     /* skip categories */
     if (do_cats) {
-	if (Map->head.coor_version.minor == 1) {	/* coor format 5.1 */
-	    if (0 >= dig__fread_port_I(&n_cats, 1, &(Map->dig_fp)))
-		return NULL;
-	}
-	else {			                /* coor format 5.0 */
-	    if (0 >= dig__fread_port_C(&nc, 1, &(Map->dig_fp)))
-		return NULL;
-	    n_cats = (int) nc;
-	}
-	G_debug(3, "    n_cats = %d", n_cats);
-
-	if (Map->head.coor_version.minor == 1) {	/* coor format 5.1 */
-	    size = (2 * PORT_INT) * n_cats;
-	}
-	else {		                /* coor format 5.0 */
-	    size = (PORT_SHORT + PORT_INT) * n_cats;
-	}
-	dig_fseek(&(Map->dig_fp), size, SEEK_CUR);
+        if (Map->head.coor_version.minor == 1) {        /* coor format 5.1 */
+            if (0 >= dig__fread_port_I(&n_cats, 1, &(Map->dig_fp)))
+                return NULL;
+        }
+        else {                                  /* coor format 5.0 */
+            if (0 >= dig__fread_port_C(&nc, 1, &(Map->dig_fp)))
+                return NULL;
+            n_cats = (int) nc;
+        }
+        G_debug(3, "    n_cats = %d", n_cats);
+
+        if (Map->head.coor_version.minor == 1) {        /* coor format 5.1 */
+            size = (2 * PORT_INT) * n_cats;
+        }
+        else {                          /* coor format 5.0 */
+            size = (PORT_SHORT + PORT_INT) * n_cats;
+        }
+        dig_fseek(&(Map->dig_fp), size, SEEK_CUR);
     }
 
     if (*type & GV_POINTS) {
-	    n_points = 1;
+            n_points = 1;
     }
     else {
-	if (0 >= dig__fread_port_I(&n_points, 1, &(Map->dig_fp)))
-	    return NULL;
+        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));
+        z = (double *) G_malloc(n_points * sizeof(double));
     else
-	z = NULL;
-    
+        z = NULL;
+
     if (0 >= dig__fread_port_D(x, n_points, &(Map->dig_fp)))
-	return NULL; /* end of file */
+        return NULL; /* end of file */
 
     if (0 >= dig__fread_port_D(y, n_points, &(Map->dig_fp)))
-	return NULL; /* end of file */
+        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 */
+        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]);
+        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", (long) dig_ftell(&(Map->dig_fp)));
-    
+
     G_free((void *) x);
     G_free((void *) y);
     if (z)
-	G_free((void *) z);
-    
+        G_free((void *) z);
+
     return pseq;
 }
 
@@ -403,17 +402,17 @@ GEOSCoordSequence *Vect_get_area_points_geos(struct Map_info *Map, int area)
 {
     struct Plus_head *Plus;
     struct P_area *Area;
-    
+
     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 */
+    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 */
     }
-    
+
     return read_polygon_points(Map, Area->n_lines, Area->lines);
 }
 
@@ -421,7 +420,7 @@ GEOSCoordSequence *Vect_get_area_points_geos(struct Map_info *Map, int area)
    \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
@@ -434,7 +433,7 @@ GEOSCoordSequence *Vect_get_isle_points_geos(struct Map_info *Map, int isle)
 {
     struct Plus_head *Plus;
     struct P_isle *Isle;
-    
+
     G_debug(3, "Vect_get_isle_points_geos(): isle = %d", isle);
 
     Plus = &(Map->plus);
@@ -450,7 +449,7 @@ GEOSCoordSequence *read_polygon_points(struct Map_info *Map, int n_lines, int *l
     unsigned int n_points, n_points_shell;
     double x, y, z;
     int *dir;
-    
+
     GEOSCoordSequence **pseq, *pseq_shell;
 
     G_debug(3, "  n_lines = %d", n_lines);
@@ -459,64 +458,64 @@ GEOSCoordSequence *read_polygon_points(struct Map_info *Map, int n_lines, int *l
 
     n_points_shell = 0;
     for (i = 0; i < n_lines; i++) {
-	line = lines[i];
-	aline = abs(line);
-	G_debug(3, "  append line(%d) = %d", i, line);
-
-	if (line > 0)
-	    dir[i] = GV_FORWARD;
-	else
-	    dir[i] = 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;
+        line = lines[i];
+        aline = abs(line);
+        G_debug(3, "  append line(%d) = %d", i, line);
+
+        if (line > 0)
+            dir[i] = GV_FORWARD;
+        else
+            dir[i] = 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 < n_lines; i++) {
-	GEOSCoordSeq_getSize(pseq[i], &n_points);
-	if (dir[i] == GV_FORWARD) {
-	    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 { /* GV_BACKWARD */
-	    for (j = (int) n_points - 1; j > -1; 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);
-		}
-	    }
-	}
-	GEOSCoordSeq_destroy(pseq[i]);
+        GEOSCoordSeq_getSize(pseq[i], &n_points);
+        if (dir[i] == GV_FORWARD) {
+            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 { /* GV_BACKWARD */
+            for (j = (int) n_points - 1; j > -1; 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);
+                }
+            }
+        }
+        GEOSCoordSeq_destroy(pseq[i]);
     }
-    
+
     G_free((void *) pseq);
     G_free((void *) dir);
-    
+
     return pseq_shell;
 }
 #endif /* HAVE_GEOS */

lib/vector/Vlib/geos_to_wktb.c

@@ -0,0 +1,195 @@
+/*!
+  \file lib/vector/Vlib/geos_to_wktb.c
+
+  \brief Vector library - GEOS powered WKT and WKB export
+
+  Higher level functions for reading/writing/manipulating vectors.
+
+  (C) 2015 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 Soeren Gebbert <soerengebbert googlemail.com>
+ */
+
+#include <stdlib.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+
+#ifdef HAVE_GEOS
+
+/*!
+   \brief Read vector area and stores it as WKB unsigned char array
+
+   \param Map pointer to Map_info structure
+   \param area area idmetry instance
+   \param size The size of the returned byte array
+   \return NULL on error
+
+   \return pointer to char array
+   \return NULL on error
+ */
+unsigned char *Vect_read_area_to_wkb(struct Map_info * Map, int area, size_t *size)
+{
+    static int init = 0;
+    /* The writer is static for performance reasons */
+    static GEOSWKBWriter *writer = NULL;
+    unsigned char *wkb = NULL;
+
+    if(init == 0) {
+        initGEOS(NULL, NULL);
+        writer = GEOSWKBWriter_create();
+        init += 1;
+    }
+
+    GEOSWKBWriter_setOutputDimension(writer, 2);
+
+    GEOSGeometry *geom = Vect_read_area_geos(Map, area);
+
+    if(!geom) {
+        return(NULL);
+    }
+
+    wkb = GEOSWKBWriter_write(writer, geom, size);
+
+    GEOSGeom_destroy(geom);
+
+    return(wkb);
+}
+
+/*!
+   \brief Read vector area and stores it as WKT char array
+
+   \param Map pointer to Map_info structure
+   \param area area idmetry instance
+   \return NULL on error
+
+   \return pointer to char array
+   \return NULL on error
+ */
+char *Vect_read_area_to_wkt(struct Map_info * Map, int area)
+{
+    static int init = 0;
+    /* The writer is static for performance reasons */
+    static GEOSWKTWriter *writer = NULL;
+    char *wkt = NULL;
+
+    if(init == 0) {
+        initGEOS(NULL, NULL);
+        writer = GEOSWKTWriter_create();
+        init += 1;
+    }
+
+    GEOSWKTWriter_setOutputDimension(writer, 2);
+
+    GEOSGeometry *geom = Vect_read_area_geos(Map, area);
+
+    if(!geom) {
+        return(NULL);
+    }
+
+    wkt = GEOSWKTWriter_write(writer, geom);
+
+    GEOSGeom_destroy(geom);
+
+    return(wkt);
+}
+
+/*!
+   \brief Create a WKB representation of given type from feature points.
+
+   This function is not thread safe, it uses static variables for speedup.
+
+   Supported types:
+   - GV_POINT    -> POINT
+   - GV_CENTROID -> POINT
+   - GV_LINE     -> LINESTRING
+   - GV_BOUNDARY -> LINEARRING
+
+   \param points pointer to line_pnts structure
+   \param type feature type (see supported types)
+   \param with_z Set to 1 if the feature is 3d, 0 otherwise
+   \param size The size of the returned byte array
+
+   \return pointer to char array
+   \return NULL on error
+ */
+unsigned char *Vect_line_to_wkb(const struct line_pnts *points,
+                       int type, int with_z, size_t *size)
+{
+    static int init = 0;
+    /* The writer is static for performance reasons */
+    static GEOSWKBWriter *writer = NULL;
+    unsigned char *wkb = NULL;
+
+    if(init == 0) {
+        initGEOS(NULL, NULL);
+        writer = GEOSWKBWriter_create();
+        init += 1;
+    }
+
+    GEOSWKBWriter_setOutputDimension(writer, with_z ? 3 : 2);
+
+    GEOSGeometry *geom = Vect_line_to_geos(points, type, with_z);
+
+    if(!geom) {
+        return(NULL);
+    }
+
+    wkb = GEOSWKBWriter_write(writer, geom, size);
+
+    GEOSGeom_destroy(geom);
+
+    return(wkb);
+}
+
+/*!
+   \brief Create a WKT representation of given type from feature points.
+
+   This function is not thread safe, it uses static variables for speedup.
+
+   Supported types:
+   - GV_POINT    -> POINT
+   - GV_CENTROID -> POINT
+   - GV_LINE     -> LINESTRING
+   - GV_BOUNDARY -> LINEARRING
+
+   \param points pointer to line_pnts structure
+   \param type feature type (see supported types)
+   \param with_z Set to 1 if the feature is 3d, 0 otherwise
+   \param with_z Set to 1 if the feature is 3d, 0 otherwise
+
+   \return pointer to char array
+   \return NULL on error
+ */
+char *Vect_line_to_wkt(const struct line_pnts *points,
+                       int type, int with_z)
+{
+    static int init = 0;
+    /* The writer is static for performance reasons */
+    static GEOSWKTWriter *writer = NULL;
+    char *wkt = NULL;
+
+    if(init == 0) {
+        initGEOS(NULL, NULL);
+        writer = GEOSWKTWriter_create();
+        init += 1;
+    }
+
+    GEOSWKTWriter_setOutputDimension(writer, with_z ? 3 : 2);
+
+    GEOSGeometry *geom = Vect_line_to_geos(points, type, with_z);
+
+    if(!geom) {
+        return(NULL);
+    }
+
+    wkt = GEOSWKTWriter_write(writer, geom);
+
+    GEOSGeom_destroy(geom);
+
+    return(wkt);
+}
+
+#endif /* HAVE_GEOS */