v.select: code reorgranization

vlib: some bugfixes in GEOS support
      (merge from devbr6, https://trac.osgeo.org/grass/changeset/37338)


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

lib/vector/Vlib/geos.c

@@ -24,6 +24,7 @@
 static GEOSGeometry *Vect__read_line_geos(struct Map_info *, long, int *);
 static GEOSCoordSequence *V1_read_line_geos(struct Map_info *, long, int *);
 static GEOSCoordSequence *V2_read_line_geos(struct Map_info *, int);
+static GEOSCoordSequence *read_polygon_points(struct Map_info *, int, int*);
 
 /*!
    \brief Read vector feature and stores it as GEOSGeometry instance
@@ -96,8 +97,10 @@ GEOSGeometry *Vect_read_area_geos(struct Map_info * Map, int 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)
+	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"),
@@ -396,16 +399,9 @@ GEOSCoordSequence *V1_read_line_geos(struct Map_info *Map, long offset, int *typ
  */
 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);
@@ -416,53 +412,7 @@ GEOSCoordSequence *Vect_get_area_points_geos(struct Map_info *Map, int 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;
+    return read_polygon_points(Map, Area->n_lines, Area->lines);
 }
 
 /*!
@@ -480,66 +430,90 @@ GEOSCoordSequence *Vect_get_area_points_geos(struct Map_info *Map, int area)
  */
 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];
+
+    return read_polygon_points(Map, Isle->n_lines, Isle->lines);
+}
+
+GEOSCoordSequence *read_polygon_points(struct Map_info *Map, int n_lines, int *lines)
+{
+    int i, j, k;
+    int line, aline;
+    unsigned int n_points, n_points_shell;
+    double x, y, z;
+    int *dir;
     
-    G_debug(3, "  n_lines = %d", Isle->n_lines);
-    for (i = 0; i < Isle->n_lines; i++) {
-	line = Isle->lines[i];
+    GEOSCoordSequence **pseq, *pseq_shell;
+
+    G_debug(3, "  n_lines = %d", n_lines);
+    pseq = (GEOSCoordSequence **) G_malloc(n_lines * sizeof(GEOSCoordSequence *));
+    dir  = (int*) G_malloc(n_lines * sizeof(int));
+
+    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;
-
-	/*
-	  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++) {
+    for (i = 0; i < 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);
+	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);
+		}
 	    }
 	}
     }
     
     G_free((void *) pseq);
-
+    G_free((void *) dir);
+    
     return pseq_shell;
 }
-
 #endif /* HAVE_GEOS */

vector/v.select/args.c

@@ -0,0 +1,87 @@
+#include <grass/gis.h>
+#include <grass/Vect.h>
+#include <grass/glocale.h>
+
+#include "proto.h"
+
+void parse_options(struct GParm *parm, struct GFlag *flag)
+{
+    parm->type[0] = G_define_standard_option(G_OPT_V_TYPE);
+    parm->type[0]->label = _("Feature type (vector map A)");
+    parm->type[0]->key = "atype";
+    parm->type[0]->guisection = _("Selection");
+
+    parm->field[0] = G_define_standard_option(G_OPT_V_FIELD);
+    parm->field[0]->label = _("Layer number (vector map A)");
+    parm->field[0]->key = "alayer";
+    parm->field[0]->guisection = _("Selection");
+
+    parm->input[1] = G_define_standard_option(G_OPT_V_INPUT);
+    parm->input[1]->description = _("Name of input vector map (B)");
+    parm->input[1]->key = "binput";
+
+    parm->type[1] = G_define_standard_option(G_OPT_V_TYPE);
+    parm->type[1]->label = _("Feature type (vector map B)");
+    parm->type[1]->key = "btype";
+    parm->type[1]->guisection = _("Selection");
+    
+    parm->field[1] = G_define_standard_option(G_OPT_V_FIELD);
+    parm->field[1]->label = _("Layer number (vector map B)");
+    parm->field[1]->key = "blayer";
+    parm->field[1]->guisection = _("Selection");
+    
+    parm->output = G_define_standard_option(G_OPT_V_OUTPUT);
+
+    parm->operator = G_define_option();
+    parm->operator->key = "operator";
+    parm->operator->type = TYPE_STRING;
+    parm->operator->required = YES;
+    parm->operator->multiple = NO;
+    parm->operator->label =
+	_("Operator defines required relation between features");
+    parm->operator->description =
+	_("A feature is written to output if the result of operation 'ainput operator binput' is true. "
+	  "An input feature is considered to be true, if category of given layer is defined.");
+#ifndef HAVE_GEOS
+    parm->operator->options = "overlaps";
+    parm->operator->answer = "overlaps";
+    parm->operator->descriptions = _("overlaps;features partially or completely overlap");
+
+    parm->relate = NULL;
+#else
+    parm->operator->options = "equals,disjoint,intersects,touches,crosses,within,contains,overlaps,relate";
+    parm->operator->descriptions = _("equals;features are spatially equals (requires flag 'g');"
+				    "disjoint;features do not spatially intersect (requires flag 'g');"
+				    "intersects;features spatially intersect (requires flag 'g');"
+				    "touches;features spatially touches (requires flag 'g');"
+				    "crosses;features spatially crosses (requires flag 'g');"
+				    "within;feature A is completely inside feature B (requires flag 'g');"
+				    "contains;feature B is completely inside feature A (requires flag 'g');"
+				    "overlaps;features spatilly overlap;"
+				    "relate;feature A is spatially related to feature B "
+				    "(requires 'relate' option and flag 'g');");
+    
+    parm->relate = G_define_option();
+    parm->relate->key = "relate";
+    parm->relate->type = TYPE_STRING;
+    parm->relate->required = NO;
+    parm->relate->multiple = NO;
+    parm->relate->description = _("Intersection Matrix Pattern used for 'relate' operator");
+#endif
+
+    flag->table = G_define_flag();
+    flag->table->key = 't';
+    flag->table->description = _("Do not create attribute table");
+
+    flag->reverse = G_define_flag();
+    flag->reverse->key = 'r';
+    flag->reverse->description = _("Reverse selection");
+    flag->reverse->guisection = _("Selection");
+#ifdef HAVE_GEOS
+    flag->geos = G_define_flag();
+    flag->geos->key = 'g';
+    flag->geos->description = _("Use GEOS operators");
+#else
+    flag->geos = NULL;
+#endif
+}

vector/v.select/main.c

@@ -34,13 +34,8 @@ int main(int argc, char *argv[])
     int **cats, *ncats, nfields, *fields;
     char *pre[2];
     struct GModule *module;
-    struct {
-	struct Option *input[2], *output, *type[2], *field[2],
-	    *operator, *relate;
-    } parm;
-    struct {
-	struct Flag *table, *reverse, *geos;
-    } flag;
+    struct GParm parm;
+    struct GFlag flag;
     struct Map_info In[2], Out;
     struct field_info *IFi, *OFi;
     struct line_pnts *APoints, *BPoints;
@@ -58,89 +53,8 @@ int main(int argc, char *argv[])
     module->description =
 	_("Selects features from vector map (A) by features from other vector map (B).");
 
-    parm.input[0] = G_define_standard_option(G_OPT_V_INPUT);
-    parm.input[0]->description = _("Name of input vector map (A)");
-    parm.input[0]->key = "ainput";
+    parse_options(&parm, &flag);
     
-    parm.type[0] = G_define_standard_option(G_OPT_V_TYPE);
-    parm.type[0]->label = _("Feature type (vector map A)");
-    parm.type[0]->key = "atype";
-    parm.type[0]->guisection = _("Selection");
-
-    parm.field[0] = G_define_standard_option(G_OPT_V_FIELD);
-    parm.field[0]->label = _("Layer number (vector map A)");
-    parm.field[0]->key = "alayer";
-    parm.field[0]->guisection = _("Selection");
-
-    parm.input[1] = G_define_standard_option(G_OPT_V_INPUT);
-    parm.input[1]->description = _("Name of input vector map (B)");
-    parm.input[1]->key = "binput";
-
-    parm.type[1] = G_define_standard_option(G_OPT_V_TYPE);
-    parm.type[1]->label = _("Feature type (vector map B)");
-    parm.type[1]->key = "btype";
-    parm.type[1]->guisection = _("Selection");
-    
-    parm.field[1] = G_define_standard_option(G_OPT_V_FIELD);
-    parm.field[1]->label = _("Layer number (vector map B)");
-    parm.field[1]->key = "blayer";
-    parm.field[1]->guisection = _("Selection");
-    
-    parm.output = G_define_standard_option(G_OPT_V_OUTPUT);
-
-    parm.operator = G_define_option();
-    parm.operator->key = "operator";
-    parm.operator->type = TYPE_STRING;
-    parm.operator->required = YES;
-    parm.operator->multiple = NO;
-    parm.operator->label =
-	_("Operator defines required relation between features");
-    parm.operator->description =
-	_("A feature is written to output if the result of operation 'ainput operator binput' is true. "
-	  "An input feature is considered to be true, if category of given layer is defined.");
-#ifndef HAVE_GEOS
-    parm.operator->options = "overlaps";
-    parm.operator->answer = "overlaps";
-    parm.operator->descriptions = _("overlaps;features partially or completely overlap");
-
-    parm.relate = NULL;
-#else
-    parm.operator->options = "equals,disjoint,intersects,touches,crosses,within,contains,overlaps,relate";
-    parm.operator->descriptions = _("equals;features are spatially equals (requires flag 'g');"
-				    "disjoint;features do not spatially intersect (requires flag 'g');"
-				    "intersects;features spatially intersect (requires flag 'g');"
-				    "touches;features spatially touches (requires flag 'g');"
-				    "crosses;features spatially crosses (requires flag 'g');"
-				    "within;feature A is completely inside feature B (requires flag 'g');"
-				    "contains;feature B is completely inside feature A (requires flag 'g');"
-				    "overlaps;features spatilly overlap;"
-				    "relate;feature A is spatially related to feature B "
-				    "(requires 'relate' option and flag 'g');");
-    
-    parm.relate = G_define_option();
-    parm.relate->key = "relate";
-    parm.relate->type = TYPE_STRING;
-    parm.relate->required = NO;
-    parm.relate->multiple = NO;
-    parm.relate->description = _("Intersection Matrix Pattern used for 'relate' operator");
-#endif
-
-    flag.table = G_define_flag();
-    flag.table->key = 't';
-    flag.table->description = _("Do not create attribute table");
-
-    flag.reverse = G_define_flag();
-    flag.reverse->key = 'r';
-    flag.reverse->description = _("Reverse selection");
-    flag.reverse->guisection = _("Selection");
-#ifdef HAVE_GEOS
-    flag.geos = G_define_flag();
-    flag.geos->key = 'g';
-    flag.geos->description = _("Use GEOS operators");
-#else
-    flag.geos = NULL;
-#endif
-
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
     

vector/v.select/proto.h

@@ -1,3 +1,6 @@
+#ifndef PROTO_H
+#define PROTO_H
+
 #define OP_EQUALS     0
 #define OP_DISJOINT   1
 #define OP_INTERSECTS 2
@@ -8,6 +11,17 @@
 #define OP_OVERLAPS   7
 #define OP_RELATE     8
 
+struct GParm {
+    struct Option *input[2], *output, *type[2], *field[2],
+	*operator, *relate;
+};
+struct GFlag {
+    struct Flag *table, *reverse, *geos;
+};
+
+/* args.c */
+void parse_options(struct GParm *, struct GFlag *);
+
 #ifdef HAVE_GEOS
 /* geos.c */
 int line_relate_geos(struct Map_info *, const GEOSGeometry *,
@@ -20,3 +34,4 @@ int area_relate_geos(struct Map_info *, const GEOSGeometry *,
 void add_aarea(struct Map_info *, int, int *);
 int line_overlap_area(struct Map_info *, int, struct Map_info *, int);
 
+#endif /* PROTO_H */

vector/v.select/v.select.html

@@ -53,6 +53,10 @@ v.select ainput=watrcrsl_eurnasia_wgs84 binput=italy_area \
          output=watrcrsl_italy operator=overlaps
 </pre></div>
 
+<h2>TODO</h2>
+ 
+Processing areas with GEOS is currently incredibly slow. Significant
+speed-up is required.
 
 <h2>SEE ALSO</h2>