[GRASS-SVN] r37339 - in grass/trunk: lib/vector/Vlib vector/v.select
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu May 21 14:48:39 EDT 2009
Author: martinl
Date: 2009-05-21 14:48:39 -0400 (Thu, 21 May 2009)
New Revision: 37339
Added:
grass/trunk/vector/v.select/args.c
Modified:
grass/trunk/lib/vector/Vlib/geos.c
grass/trunk/vector/v.select/main.c
grass/trunk/vector/v.select/proto.h
grass/trunk/vector/v.select/v.select.html
Log:
v.select: code reorgranization
vlib: some bugfixes in GEOS support
(merge from devbr6, r37338)
Modified: grass/trunk/lib/vector/Vlib/geos.c
===================================================================
--- grass/trunk/lib/vector/Vlib/geos.c 2009-05-21 18:43:04 UTC (rev 37338)
+++ grass/trunk/lib/vector/Vlib/geos.c 2009-05-21 18:48:39 UTC (rev 37339)
@@ -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 @@
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 *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 @@
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_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 */
Copied: grass/trunk/vector/v.select/args.c (from rev 37338, grass/branches/develbranch_6/vector/v.select/args.c)
===================================================================
--- grass/trunk/vector/v.select/args.c (rev 0)
+++ grass/trunk/vector/v.select/args.c 2009-05-21 18:48:39 UTC (rev 37339)
@@ -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
+}
Modified: grass/trunk/vector/v.select/main.c
===================================================================
--- grass/trunk/vector/v.select/main.c 2009-05-21 18:43:04 UTC (rev 37338)
+++ grass/trunk/vector/v.select/main.c 2009-05-21 18:48:39 UTC (rev 37339)
@@ -34,13 +34,8 @@
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 @@
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);
Modified: grass/trunk/vector/v.select/proto.h
===================================================================
--- grass/trunk/vector/v.select/proto.h 2009-05-21 18:43:04 UTC (rev 37338)
+++ grass/trunk/vector/v.select/proto.h 2009-05-21 18:48:39 UTC (rev 37339)
@@ -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 @@
void add_aarea(struct Map_info *, int, int *);
int line_overlap_area(struct Map_info *, int, struct Map_info *, int);
+#endif /* PROTO_H */
Modified: grass/trunk/vector/v.select/v.select.html
===================================================================
--- grass/trunk/vector/v.select/v.select.html 2009-05-21 18:43:04 UTC (rev 37338)
+++ grass/trunk/vector/v.select/v.select.html 2009-05-21 18:48:39 UTC (rev 37339)
@@ -53,6 +53,10 @@
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>
More information about the grass-commit
mailing list