[GRASS-SVN] r38385 - grass/trunk/lib/vector/Vlib
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Jul 13 08:27:16 EDT 2009
Author: mmetz
Date: 2009-07-13 08:27:16 -0400 (Mon, 13 Jul 2009)
New Revision: 38385
Modified:
grass/trunk/lib/vector/Vlib/build.c
grass/trunk/lib/vector/Vlib/build_nat.c
grass/trunk/lib/vector/Vlib/close.c
grass/trunk/lib/vector/Vlib/intersect.c
grass/trunk/lib/vector/Vlib/open.c
grass/trunk/lib/vector/Vlib/select.c
grass/trunk/lib/vector/Vlib/sindex.c
Log:
Vlib new spatial index
Modified: grass/trunk/lib/vector/Vlib/build.c
===================================================================
--- grass/trunk/lib/vector/Vlib/build.c 2009-07-13 12:26:17 UTC (rev 38384)
+++ grass/trunk/lib/vector/Vlib/build.c 2009-07-13 12:27:16 UTC (rev 38385)
@@ -19,6 +19,7 @@
#include <stdlib.h>
#include <stdio.h>
#include <stdarg.h>
+#include <unistd.h>
#include <grass/glocale.h>
#include <grass/gis.h>
#include <grass/vector.h>
@@ -118,7 +119,8 @@
Map->level = 1; /* may be not needed, because V1_read is used directly by Vect_build_ */
Map->support_updated = 1;
- Map->plus.Spidx_built = 1;
+ if (Map->plus.Spidx_built == 0)
+ Vect_open_sidx(Map, 2);
plus = &(Map->plus);
if (build > GV_BUILD_NONE) {
@@ -161,11 +163,11 @@
if (plus->n_flines > 0)
G_message(_("Number of faces: %d"), plus->n_flines);
-
+
if (plus->n_klines > 0)
G_message(_("Number of kernels: %d"), plus->n_klines);
}
-
+
if (plus->built >= GV_BUILD_AREAS) {
int line, nlines, area, nareas, err_boundaries, err_centr_out,
err_centr_dupl, err_nocentr;
@@ -209,19 +211,17 @@
if (err_boundaries)
G_message(_("Number of incorrect boundaries: %d"),
- err_boundaries);
+ err_boundaries);
if (err_centr_out)
G_message(_("Number of centroids outside area: %d"),
- err_centr_out);
+ err_centr_out);
if (err_centr_dupl)
- G_message(_("Number of duplicate centroids: %d"),
- err_centr_dupl);
+ G_message(_("Number of duplicate centroids: %d"), err_centr_dupl);
if (err_nocentr)
- G_message(_("Number of areas without centroid: %d"),
- err_nocentr);
+ G_message(_("Number of areas without centroid: %d"), err_nocentr);
}
else if (build > GV_BUILD_NONE) {
@@ -327,8 +327,8 @@
Line = plus->Line[i];
fprintf(out, "line = %d, type = %d, offset = %lu n1 = %d, n2 = %d, "
"left/area = %d, right = %d\n",
- i, Line->type, (unsigned long) Line->offset, Line->N1, Line->N2,
- Line->left, Line->right);
+ i, Line->type, (unsigned long)Line->offset, Line->N1,
+ Line->N2, Line->left, Line->right);
fprintf(out, "N,S,E,W,T,B: %f, %f, %f, %f, %f, %f\n", Line->N,
Line->S, Line->E, Line->W, Line->T, Line->B);
}
@@ -383,44 +383,180 @@
}
/*!
- \brief Save spatial index file
+ \brief Create spatial index if necessary.
+ To be used in modules.
+ Map must be opened on level 2.
+
+ \param[in,out] Map pointer to vector map
+
+ \return 0 OK
+ \return 1 error
+ */
+int Vect_build_sidx(struct Map_info *Map)
+{
+ if (Map->level < 2) {
+ G_fatal_error(_("Unable to build spatial index from topology, "
+ "vector map is not opened at topology level 2"));
+ }
+ if (!(Map->plus.Spidx_built)) {
+ return (Vect_build_sidx_from_topo(Map));
+ }
+ return 0;
+}
+
+/*!
+ \brief Create spatial index from topology if necessary
+
+ \param Map pointer to vector map
+
+ \return 0 OK
+ \return 1 error
+ */
+int Vect_build_sidx_from_topo(struct Map_info *Map)
+{
+ int i, total, done;
+ struct Plus_head *plus;
+ BOUND_BOX box;
+ P_LINE *Line;
+ P_NODE *Node;
+ P_AREA *Area;
+ P_ISLE *Isle;
+
+ G_debug(3, "Vect_build_sidx_from_topo()");
+
+ plus = &(Map->plus);
+
+ Vect_open_sidx(Map, 2);
+
+ total = plus->n_nodes + plus->n_lines + plus->n_areas + plus->n_isles;
+
+ /* Nodes */
+ for (i = 1; i <= plus->n_nodes; i++) {
+ G_percent(i, total, 3);
+
+ Node = plus->Node[i];
+ if (!Node)
+ G_fatal_error(_("BUG (Vect_build_sidx_from_topo): node does not exist"));
+
+ dig_spidx_add_node(plus, i, Node->x, Node->y, Node->z);
+ }
+
+ /* Lines */
+ done = plus->n_nodes;
+ for (i = 1; i <= plus->n_lines; i++) {
+ G_percent(done + i, total, 3);
+
+ Line = plus->Line[i];
+ if (!Line)
+ G_fatal_error(_("BUG (Vect_build_sidx_from_topo): line does not exist"));
+
+ box.N = Line->N;
+ box.S = Line->S;
+ box.E = Line->E;
+ box.W = Line->W;
+ box.T = Line->T;
+ box.B = Line->B;
+
+ dig_spidx_add_line(plus, i, &box);
+ }
+
+ /* Areas */
+ done += plus->n_lines;
+ for (i = 1; i <= plus->n_areas; i++) {
+ G_percent(done + i, total, 3);
+
+ Area = plus->Area[i];
+ if (!Area)
+ G_fatal_error(_("BUG (Vect_build_sidx_from_topo): area does not exist"));
+
+ box.N = Area->N;
+ box.S = Area->S;
+ box.E = Area->E;
+ box.W = Area->W;
+ box.T = Area->T;
+ box.B = Area->B;
+
+ dig_spidx_add_area(plus, i, &box);
+ }
+
+ /* Isles */
+ done += plus->n_areas;
+ for (i = 1; i <= plus->n_isles; i++) {
+ G_percent(done + i, total, 3);
+
+ Isle = plus->Isle[i];
+ if (!Isle)
+ G_fatal_error(_("BUG (Vect_build_sidx_from_topo): isle does not exist"));
+
+ box.N = Isle->N;
+ box.S = Isle->S;
+ box.E = Isle->E;
+ box.W = Isle->W;
+ box.T = Isle->T;
+ box.B = Isle->B;
+
+ dig_spidx_add_isle(plus, i, &box);
+ }
+
+ Map->plus.Spidx_built = 1;
+
+ G_debug(3, "Spatial index was built");
+
+ return 0;
+}
+
+/*!
+ \brief Save spatial index file for vector map
+
\param Map vector map
\return 1 on success
\return 0 on error
*/
-int Vect_save_spatial_index(struct Map_info *Map)
+int Vect_save_sidx(struct Map_info *Map)
{
struct Plus_head *plus;
char fname[GPATH_MAX], buf[GPATH_MAX];
- GVFILE fp;
G_debug(1, "Vect_save_spatial_index()");
plus = &(Map->plus);
- /* write out rtrees to the sidx file */
- sprintf(buf, "%s/%s", GRASS_VECT_DIRECTORY, Map->name);
- G__file_name(fname, buf, GV_SIDX_ELEMENT, Map->mapset);
- G_debug(1, "Open sidx: %s", fname);
- dig_file_init(&fp);
- fp.file = fopen(fname, "w");
- if (fp.file == NULL) {
- G_warning(_("Unable open spatial index file for write <%s>"), fname);
+ if (Map->plus.Spidx_built == 0) {
+ G_warning("Spatial index not available, can not be saved");
return 0;
}
- /* set portable info */
- dig_init_portable(&(plus->spidx_port), dig__byte_order_out());
+ /* new or update mode ? */
+ if (Map->plus.Spidx_new == 1) {
- if (0 > dig_write_spidx(&fp, plus)) {
- G_warning(_("Error writing out spatial index file"));
- return 0;
+ /* write out rtrees to sidx file */
+ sprintf(buf, "%s/%s", GRASS_VECT_DIRECTORY, Map->name);
+ G__file_name(fname, buf, GV_SIDX_ELEMENT, Map->mapset);
+ G_debug(1, "Open sidx: %s", fname);
+ dig_file_init(&(Map->plus.spidx_fp));
+ Map->plus.spidx_fp.file = fopen(fname, "w+");
+ if (Map->plus.spidx_fp.file == NULL) {
+ G_warning(_("Unable open spatial index file for write <%s>"),
+ fname);
+ return 0;
+ }
+
+ /* set portable info */
+ dig_init_portable(&(plus->spidx_port), dig__byte_order_out());
+
+ if (0 > dig_Wr_spidx(&(Map->plus.spidx_fp), plus)) {
+ G_warning(_("Error writing out spatial index file"));
+ return 0;
+ }
+ Map->plus.Spidx_new = 0;
}
- fclose(fp.file);
+ fclose(Map->plus.spidx_fp.file);
+ Map->plus.Spidx_built = 0;
+
return 1;
}
@@ -433,7 +569,7 @@
\return 1 on success
\return 0 on error
*/
-int Vect_spatial_index_dump(struct Map_info *Map, FILE * out)
+int Vect_sidx_dump(struct Map_info *Map, FILE * out)
{
if (!(Map->plus.Spidx_built)) {
Vect_build_sidx_from_topo(Map);
Modified: grass/trunk/lib/vector/Vlib/build_nat.c
===================================================================
--- grass/trunk/lib/vector/Vlib/build_nat.c 2009-07-13 12:26:17 UTC (rev 38384)
+++ grass/trunk/lib/vector/Vlib/build_nat.c 2009-07-13 12:27:16 UTC (rev 38385)
@@ -73,16 +73,17 @@
line = abs(lines[j]);
BLine = plus->Line[line];
offset = BLine->offset;
- G_debug(3, " line[%d] = %d, offset = %lu", j, line, (unsigned long) offset);
+ G_debug(3, " line[%d] = %d, offset = %lu", j, line,
+ (unsigned long)offset);
type = Vect_read_line(Map, Points, NULL, line);
if (lines[j] > 0)
direction = GV_FORWARD;
else
direction = GV_BACKWARD;
Vect_append_points(APoints, Points, direction);
- APoints->n_points--; /* skip last point, avoids duplicates */
+ APoints->n_points--; /* skip last point, avoids duplicates */
}
- APoints->n_points++; /* close polygon */
+ APoints->n_points++; /* close polygon */
dig_find_area_poly(APoints, &area_size);
@@ -220,9 +221,9 @@
/* This is slow, but should not be called often */
Vect_get_area_points(Map, sel_area, APoints);
/* G_begin_polygon_area_calculations();
- cur_size =
- G_area_of_polygon(APoints->x, APoints->y,
- APoints->n_points); */
+ cur_size =
+ G_area_of_polygon(APoints->x, APoints->y,
+ APoints->n_points); */
/* this is faster, but there may be latlon problems: the poles */
dig_find_area_poly(APoints, &cur_size);
G_debug(3, " first area size = %f (n points = %d)",
@@ -232,8 +233,8 @@
Vect_get_area_points(Map, area, APoints);
/* size =
- G_area_of_polygon(APoints->x, APoints->y,
- APoints->n_points); */
+ G_area_of_polygon(APoints->x, APoints->y,
+ APoints->n_points); */
/* this is faster, but there may be latlon problems: the poles */
dig_find_area_poly(APoints, &size);
G_debug(3, " area size = %f (n points = %d)", size,
@@ -394,7 +395,7 @@
* 3) it's faster
*/
if (Line->left > 0)
- continue;
+ continue;
orig_area = Line->left;
@@ -539,7 +540,7 @@
offset = Map->head.last_offset;
- G_debug(3, "Register line: offset = %lu", (unsigned long) offset);
+ G_debug(3, "Register line: offset = %lu", (unsigned long)offset);
lineid = dig_add_line(plus, type, Points, offset);
dig_line_box(Points, &box);
if (lineid == 1)
@@ -565,11 +566,11 @@
else
fprintf(stderr, "%11d\b\b\b\b\b\b\b\b\b\b\b", i);
}
-
+
i++;
}
-
- if ( (G_verbose() > G_verbose_min() ) && format != G_INFO_FORMAT_PLAIN )
+
+ if ((G_verbose() > G_verbose_min()) && format != G_INFO_FORMAT_PLAIN)
fprintf(stderr, "\r");
G_message(_("%d primitives registered"), plus->n_lines);
Modified: grass/trunk/lib/vector/Vlib/close.c
===================================================================
--- grass/trunk/lib/vector/Vlib/close.c 2009-07-13 12:26:17 UTC (rev 38384)
+++ grass/trunk/lib/vector/Vlib/close.c 2009-07-13 12:27:16 UTC (rev 38385)
@@ -12,7 +12,7 @@
\author Original author CERL, probably Dave Gerdes or Mike Higgins.
\author Update to GRASS 5.7 Radim Blazek and David D. Gray.
-*/
+ */
#include <grass/config.h>
#include <stdlib.h>
@@ -94,8 +94,7 @@
Vect_save_topo(Map);
- /* Spatial index is not saved */
- /* Vect_save_spatial_index ( Map ); */
+ Vect_save_sidx(Map);
Vect_cidx_save(Map);
@@ -104,15 +103,18 @@
V2_close_ogr(Map);
#endif
}
+ else {
+ /* spatial index must also be closed when opened with topo but not modified */
+ if (Map->plus.Spidx_built == 1 && Map->plus.built == GV_BUILD_ALL)
+ Vect_save_sidx(Map);
+ }
if (Map->level == 2 && Map->plus.release_support) {
G_debug(1, "free topology");
dig_free_plus(&(Map->plus));
- if (!Map->head_only) {
- G_debug(1, "free spatial index");
- dig_spidx_free(&(Map->plus));
- }
+ G_debug(1, "free spatial index");
+ dig_spidx_free(&(Map->plus));
G_debug(1, "free category index");
dig_cidx_free(&(Map->plus));
Modified: grass/trunk/lib/vector/Vlib/intersect.c
===================================================================
--- grass/trunk/lib/vector/Vlib/intersect.c 2009-07-13 12:26:17 UTC (rev 38384)
+++ grass/trunk/lib/vector/Vlib/intersect.c 2009-07-13 12:27:16 UTC (rev 38385)
@@ -6,49 +6,49 @@
Higher level functions for reading/writing/manipulating vectors.
Some parts of code taken from grass50 v.spag/linecros.c
-
+
Based on the following:
-
+
<code>
(ax2-ax1)r1 - (bx2-bx1)r2 = ax2 - ax1
(ay2-ay1)r1 - (by2-by1)r2 = ay2 - ay1
</code>
-
+
Solving for r1 and r2, if r1 and r2 are between 0 and 1, then line
segments (ax1,ay1)(ax2,ay2) and (bx1,by1)(bx2,by2) intersect.
Intersect 2 line segments.
-
+
Returns: 0 - do not intersect
- 1 - intersect at one point
+ 1 - intersect at one point
<pre>
- \ / \ / \ /
- \/ \/ \/
- /\ \
- / \ \
- 2 - partial overlap ( \/ )
- ------ a ( distance < threshold )
- ------ b ( )
- 3 - a contains b ( /\ )
- ---------- a ----------- a
- ---- b ----- b
- 4 - b contains a
- ---- a ----- a
- ---------- b ----------- b
- 5 - identical
- ---------- a
- ---------- b
- </pre>
- Intersection points:
- <pre>
- return point1 breakes: point2 breaks: distance1 on: distance2 on:
- 0 - - - -
- 1 a,b - a b
- 2 a b a b
- 3 a a a a
- 4 b b b b
- 5 - - - -
- </pre>
+ \ / \ / \ /
+ \/ \/ \/
+ /\ \
+ / \ \
+ 2 - partial overlap ( \/ )
+ ------ a ( distance < threshold )
+ ------ b ( )
+ 3 - a contains b ( /\ )
+ ---------- a ----------- a
+ ---- b ----- b
+ 4 - b contains a
+ ---- a ----- a
+ ---------- b ----------- b
+ 5 - identical
+ ---------- a
+ ---------- b
+ </pre>
+ Intersection points:
+ <pre>
+ return point1 breakes: point2 breaks: distance1 on: distance2 on:
+ 0 - - - -
+ 1 a,b - a b
+ 2 a b a b
+ 3 a a a a
+ 4 b b b b
+ 5 - - - -
+ </pre>
Sometimes (often) is important to get the same coordinates for a x
b and b x a. To reach this, the segments a,b are 'sorted' at the
beginning, so that for the same switched segments, results are
@@ -67,6 +67,8 @@
#include <grass/config.h>
#include <stdlib.h>
+#include <stdio.h>
+#include <unistd.h>
#include <math.h>
#include <grass/gis.h>
#include <grass/vector.h>
@@ -608,7 +610,7 @@
double x, y, rethresh;
struct Rect rect;
struct line_pnts **XLines, *Points;
- struct Node *RTree;
+ struct RTree *MyRTree;
struct line_pnts *Points1, *Points2; /* first, second points */
int seg1, seg2, vert1, vert2;
@@ -675,7 +677,7 @@
* in bound box */
/* Create rtree for B line */
- RTree = RTreeNewIndex();
+ MyRTree = RTreeNewIndex(2);
for (i = 0; i < BPoints->n_points - 1; i++) {
if (BPoints->x[i] <= BPoints->x[i + 1]) {
rect.boundary[0] = BPoints->x[i];
@@ -704,7 +706,7 @@
rect.boundary[5] = BPoints->z[i];
}
- RTreeInsertRect(&rect, i + 1, &RTree, 0); /* B line segment numbers in rtree start from 1 */
+ RTreeInsertRect(&rect, i + 1, MyRTree); /* B line segment numbers in rtree start from 1 */
}
/* Break segments in A by segments in B */
@@ -735,11 +737,11 @@
rect.boundary[5] = APoints->z[i];
}
- j = RTreeSearch(RTree, &rect, (void *)cross_seg, &i); /* A segment number from 0 */
+ j = RTreeSearch(MyRTree, &rect, (void *)cross_seg, &i); /* A segment number from 0 */
}
/* Free RTree */
- RTreeDestroyNode(RTree);
+ RTreeFreeIndex(MyRTree);
G_debug(2, "n_cross = %d", n_cross);
/* Lines do not cross each other */
@@ -1087,7 +1089,7 @@
if (!IPnts)
IPnts = Vect_new_line_struct();
-
+
switch (ret) {
case 0:
case 5:
@@ -1132,7 +1134,7 @@
int i;
double dist, rethresh;
struct Rect rect;
- struct Node *RTree;
+ struct RTree *MyRTree;
rethresh = 0.000001; /* TODO */
APnts = APoints;
@@ -1212,7 +1214,7 @@
* in bound box */
/* Create rtree for B line */
- RTree = RTreeNewIndex();
+ MyRTree = RTreeNewIndex(2);
for (i = 0; i < BPoints->n_points - 1; i++) {
if (BPoints->x[i] <= BPoints->x[i + 1]) {
rect.boundary[0] = BPoints->x[i];
@@ -1241,7 +1243,7 @@
rect.boundary[5] = BPoints->z[i];
}
- RTreeInsertRect(&rect, i + 1, &RTree, 0); /* B line segment numbers in rtree start from 1 */
+ RTreeInsertRect(&rect, i + 1, MyRTree); /* B line segment numbers in rtree start from 1 */
}
/* Find intersection */
@@ -1274,7 +1276,7 @@
rect.boundary[5] = APoints->z[i];
}
- RTreeSearch(RTree, &rect, (void *)find_cross, &i); /* A segment number from 0 */
+ RTreeSearch(MyRTree, &rect, (void *)find_cross, &i); /* A segment number from 0 */
if (cross_found) {
break;
@@ -1282,7 +1284,7 @@
}
/* Free RTree */
- RTreeDestroyNode(RTree);
+ RTreeFreeIndex(MyRTree);
return cross_found;
}
Modified: grass/trunk/lib/vector/Vlib/open.c
===================================================================
--- grass/trunk/lib/vector/Vlib/open.c 2009-07-13 12:26:17 UTC (rev 38384)
+++ grass/trunk/lib/vector/Vlib/open.c 2009-07-13 12:27:16 UTC (rev 38385)
@@ -212,7 +212,8 @@
_("Unable to open vector map <%s> on level %d. "
"Try to rebuild vector topology by v.build."),
Vect_get_full_name(Map), level_request);
- G_warning(_("Unable to read head file of vector <%s>"), Vect_get_full_name(Map));
+ G_warning(_("Unable to read head file of vector <%s>"),
+ Vect_get_full_name(Map));
}
G_debug(1, "Level request = %d", level_request);
@@ -229,7 +230,7 @@
/* open topo */
ret = Vect_open_topo(Map, head_only);
if (ret == 1) { /* topo file is not available */
- G_debug(1, "Topo file for vector '%s' not available.",
+ G_debug(1, "topo file for vector '%s' not available.",
Vect_get_full_name(Map));
level = 1;
}
@@ -237,23 +238,25 @@
G_fatal_error(_("Unable to open topology file for vector map <%s>"),
Vect_get_full_name(Map));
}
- /* open spatial index, not needed for head_only */
- /* spatial index is not loaded anymore */
- /*
- if ( level == 2 && !head_only ) {
- if ( Vect_open_spatial_index(Map) == -1 ) {
- G_debug( 1, "Cannot open spatial index file for vector '%s'.", Vect_get_full_name (Map) );
- dig_free_plus ( &(Map->plus) );
- level = 1;
- }
- }
- */
+ /* open spatial index */
+ if (level == 2) {
+ ret = Vect_open_sidx(Map, (update != 0));
+ if (ret == 1) { /* sidx file is not available */
+ G_debug(1, "sidx file for vector '%s' not available.",
+ Vect_get_full_name(Map));
+ level = 1;
+ }
+ else if (ret == -1) {
+ G_fatal_error(_("Unable to open spatial index file for vector map <%s>"),
+ Vect_get_full_name(Map));
+ }
+ }
/* open category index */
if (level == 2) {
ret = Vect_cidx_open(Map, head_only);
if (ret == 1) { /* category index is not available */
G_debug(1,
- "Category index file for vector '%s' not available.",
+ "cidx file for vector '%s' not available.",
Vect_get_full_name(Map));
dig_free_plus(&(Map->plus)); /* free topology */
dig_spidx_free(&(Map->plus)); /* free spatial index */
@@ -356,7 +359,7 @@
fatal_error(ferror, errmsg);
return (-1);
}
- fseek(Map->hist_fp, (off_t)0, SEEK_END);
+ fseek(Map->hist_fp, (off_t) 0, SEEK_END);
Vect_hist_write(Map,
"---------------------------------------------------------------------------------\n");
@@ -445,8 +448,10 @@
Map->plus.n_upnodes = 0;
Map->plus.alloc_upnodes = 0;
+ /* read spatial index */
+
/* Build spatial index from topo */
- Vect_build_sidx_from_topo(Map);
+ /* Vect_build_sidx_from_topo(Map); */
}
return ret;
@@ -518,7 +523,7 @@
int Vect_open_new(struct Map_info *Map, const char *name, int with_z)
{
int ret, ferror;
- char errmsg[2000], buf[200];
+ char errmsg[2000], buf[500];
char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
G_debug(2, "Vect_open_new(): name = %s", name);
@@ -582,8 +587,12 @@
Open_level = 0;
+ /* initialize topo */
dig_init_plus(&(Map->plus));
+ /* initialize spatial index */
+ Vect_open_sidx(Map, 2);
+
Map->open = VECT_OPEN_CODE;
Map->level = 1;
Map->head_only = 0;
@@ -625,7 +634,7 @@
Info->mtime = -1L;
}
else {
- Info->size = (off_t)stat_buf.st_size; /* file size */
+ Info->size = (off_t) stat_buf.st_size; /* file size */
Info->mtime = (long)stat_buf.st_mtime; /* last modified time */
}
/* stat does not give correct size on MINGW
@@ -643,8 +652,8 @@
Info->mtime = 0L;
break;
}
- G_debug(1, "Info->size = %lu, Info->mtime = %ld", (unsigned long) Info->size,
- Info->mtime);
+ G_debug(1, "Info->size = %lu, Info->mtime = %ld",
+ (unsigned long)Info->size, Info->mtime);
return 1;
}
@@ -727,7 +736,7 @@
return -1;
G_debug(1, "Topo head: coor size = %lu, coor mtime = %ld",
- (unsigned long) Plus->coor_size, Plus->coor_mtime);
+ (unsigned long)Plus->coor_size, Plus->coor_mtime);
/* do checks */
err = 0;
@@ -767,73 +776,101 @@
* \brief Open spatial index file
*
* \param[in,out] Map vector map
+ * \param mode 0 old, 1 update, 2 new
*
* \return 0 on success
* \return -1 on error
*/
-int Vect_open_spatial_index(struct Map_info *Map)
+int Vect_open_sidx(struct Map_info *Map, int mode)
{
- char buf[500];
- GVFILE fp;
-
- /* struct Coor_info CInfo; */
+ char buf[500], spidxbuf[500], file_path[2000];
+ int err;
+ struct Coor_info CInfo;
struct Plus_head *Plus;
+ struct stat info;
- G_debug(1, "Vect_open_spatial_index(): name = %s mapset= %s", Map->name,
- Map->mapset);
+ G_debug(1, "Vect_open_sidx(): name = %s mapset= %s mode = %s", Map->name,
+ Map->mapset, mode == 0 ? "old" : (mode == 1 ? "update" : "new"));
+ if (Map->plus.Spidx_built == 1) {
+ G_warning("spatial index already opened");
+ return 0;
+ }
+
Plus = &(Map->plus);
- sprintf(buf, "%s/%s", GRASS_VECT_DIRECTORY, Map->name);
- dig_file_init(&fp);
- fp.file = G_fopen_old(buf, GV_SIDX_ELEMENT, Map->mapset);
+ dig_file_init(&(Map->plus.spidx_fp));
- if (fp.file == NULL) { /* spatial index file is not available */
- G_debug(1, "Cannot open spatial index file for vector '%s@%s'.",
- Map->name, Map->mapset);
- return -1;
+ if (mode < 2) {
+ sprintf(buf, "%s/%s", GRASS_VECT_DIRECTORY, Map->name);
+ G__file_name(file_path, buf, GV_SIDX_ELEMENT, Map->mapset);
+
+ if (stat(file_path, &info) != 0) /* does not exist */
+ return 1;
+
+ Map->plus.spidx_fp.file =
+ G_fopen_old(buf, GV_SIDX_ELEMENT, Map->mapset);
+
+ if (Map->plus.spidx_fp.file == NULL) { /* sidx file is not available */
+ G_debug(1, "Cannot open spatial index file for vector '%s@%s'.",
+ Map->name, Map->mapset);
+ return -1;
+ }
+
+ /* get coor info */
+ Vect_coor_info(Map, &CInfo);
+
+ /* initialize spatial index */
+ Map->plus.Spidx_new = 0;
+
+ dig_spidx_init(Plus);
+
+ /* load head */
+ if (dig_Rd_spidx_head(&(Map->plus.spidx_fp), Plus) == -1) {
+ fclose(Map->plus.spidx_fp.file);
+ return -1;
+ }
+
+ G_debug(1, "Sidx head: coor size = %lu, coor mtime = %ld",
+ (unsigned long)Plus->coor_size, Plus->coor_mtime);
+
+ /* do checks */
+ err = 0;
+ if (CInfo.size != Plus->coor_size) {
+ G_warning(_("Size of 'coor' file differs from value saved in sidx file"));
+ err = 1;
+ }
+ /* Do not check mtime because mtime is changed by copy */
+ /*
+ if ( CInfo.mtime != Plus->coor_mtime ) {
+ G_warning ( "Time of last modification for 'coor' file differs from value saved in topo file.\n");
+ err = 1;
+ }
+ */
+ if (err) {
+ G_warning(_("Please rebuild topology for vector map <%s@%s>"),
+ Map->name, Map->mapset);
+ fclose(Map->plus.spidx_fp.file);
+ return -1;
+ }
}
- /* TODO: checks */
- /* load head */
- /*
- dig_Rd_spindx_head (fp, Plus);
- G_debug ( 1, "Spindx head: coor size = %ld, coor mtime = %ld",
- Plus->coor_size, Plus->coor_mtime);
+ if (mode) {
+ /* open new spatial index */
+ Map->plus.Spidx_new = 1;
- */
- /* do checks */
- /*
- err = 0;
- if ( CInfo.size != Plus->coor_size ) {
- G_warning ( "Size of 'coor' file differs from value saved in topo file.\n");
- err = 1;
- }
- */
- /* Do not check mtime because mtime is changed by copy */
- /*
- if ( CInfo.mtime != Plus->coor_mtime ) {
- G_warning ( "Time of last modification for 'coor' file differs from value saved in topo file.\n");
- err = 1;
- }
- */
- /*
- if ( err ) {
- G_warning ( "Please rebuild topology for vector '%s@%s'\n", Map->name,
- Map->mapset );
- return -1;
- }
- */
+ dig_spidx_init(Plus);
- /* load file to the memory */
- /* dig_file_load ( &fp); */
+ if (mode == 1) {
+ /* load spatial index for update */
+ if (dig_Rd_spidx(&(Map->plus.spidx_fp), Plus) == -1) {
+ fclose(Map->plus.spidx_fp.file);
+ return -1;
+ }
+ }
+ }
- /* load topo to memory */
- dig_spidx_init(Plus);
- dig_read_spidx(&fp, Plus);
+ Map->plus.Spidx_built = 1;
- fclose(fp.file);
- /* dig_file_free ( &fp); */
-
return 0;
}
Modified: grass/trunk/lib/vector/Vlib/select.c
===================================================================
--- grass/trunk/lib/vector/Vlib/select.c 2009-07-13 12:26:17 UTC (rev 38384)
+++ grass/trunk/lib/vector/Vlib/select.c 2009-07-13 12:27:16 UTC (rev 38385)
@@ -1,9 +1,9 @@
/*!
- \file select.c
+ \file Vlib/select.c
- \brief Vector library - select vector features
+ \brief Vector library - spatial index
- Higher level functions for reading/writing/manipulating vectors.
+ Higher level functions for a custom spatial index.
(C) 2001-2009 by the GRASS Development Team
@@ -15,345 +15,134 @@
#include <grass/config.h>
#include <stdlib.h>
+#include <unistd.h>
+#include <sys/stat.h>
+#include <string.h>
+#include <grass/glocale.h>
#include <grass/gis.h>
#include <grass/vector.h>
+#include <grass/glocale.h>
/*!
- \brief Select lines by box.
+ \brief Initialize spatial index structure
- Select lines whose boxes overlap specified box!!! It means that
- selected line may or may not overlap the box.
+ \param si pointer to spatial index structure
- \param Map vector map
- \param Box bounding box
- \param type line type
- \param[out] list output list, must be initialized
-
- \return number of lines
+ \return void
*/
-int
-Vect_select_lines_by_box(struct Map_info *Map, const BOUND_BOX * Box,
- int type, struct ilist *list)
+void Vect_spatial_index_init(SPATIAL_INDEX * si, int with_z)
{
- int i, line, nlines;
- struct Plus_head *plus;
- P_LINE *Line;
- static struct ilist *LocList = NULL;
+ G_debug(1, "Vect_spatial_index_init()");
- G_debug(3, "Vect_select_lines_by_box()");
- G_debug(3, " Box(N,S,E,W,T,B): %e, %e, %e, %e, %e, %e", Box->N, Box->S,
- Box->E, Box->W, Box->T, Box->B);
- plus = &(Map->plus);
-
- if (!(plus->Spidx_built)) {
- G_debug(3, "Building spatial index.");
- Vect_build_sidx_from_topo(Map);
- }
-
- list->n_values = 0;
- if (!LocList)
- LocList = Vect_new_list();
-
- nlines = dig_select_lines(plus, Box, LocList);
- G_debug(3, " %d lines selected (all types)", nlines);
-
- /* Remove lines of not requested types */
- for (i = 0; i < nlines; i++) {
- line = LocList->value[i];
- if (plus->Line[line] == NULL)
- continue; /* Should not happen */
- Line = plus->Line[line];
- if (!(Line->type & type))
- continue;
- dig_list_add(list, line);
- }
-
- G_debug(3, " %d lines of requested type", list->n_values);
-
- return list->n_values;
+ si->si_tree = RTreeNewIndex(2 + with_z);
}
/*!
- \brief Select areas by box.
+ \brief Destroy existing spatial index
- Select areas whose boxes overlap specified box!!!
- It means that selected area may or may not overlap the box.
+ Vect_spatial_index_init() must be call before new use.
- \param Map vector map
- \param Box bounding box
- \param[out] output list, must be initialized
+ \param si pointer to spatial index structure
- \return number of areas
+ \return void
*/
-int
-Vect_select_areas_by_box(struct Map_info *Map, const BOUND_BOX * Box,
- struct ilist *list)
+void Vect_spatial_index_destroy(SPATIAL_INDEX * si)
{
- int i;
- const char *dstr;
- int debug_level;
+ G_debug(1, "Vect_spatial_index_destroy()");
- dstr = G__getenv("DEBUG");
-
- if (dstr != NULL)
- debug_level = atoi(dstr);
- else
- debug_level = 0;
-
- G_debug(3, "Vect_select_areas_by_box()");
- G_debug(3, "Box(N,S,E,W,T,B): %e, %e, %e, %e, %e, %e", Box->N, Box->S,
- Box->E, Box->W, Box->T, Box->B);
-
- if (!(Map->plus.Spidx_built)) {
- G_debug(3, "Building spatial index.");
- Vect_build_sidx_from_topo(Map);
- }
-
- dig_select_areas(&(Map->plus), Box, list);
- G_debug(3, " %d areas selected", list->n_values);
- /* avoid loop when not debugging */
- if (debug_level > 2) {
- for (i = 0; i < list->n_values; i++) {
- G_debug(3, " area = %d pointer to area structure = %lx",
- list->value[i],
- (unsigned long)Map->plus.Area[list->value[i]]);
- }
- }
-
- return list->n_values;
+ RTreeFreeIndex(si->si_tree);
}
-
/*!
- \brief Select isles by box.
+ \brief Add a new item to spatial index structure
- Select isles whose boxes overlap specified box!!!
- It means that selected isle may or may not overlap the box.
+ \param[in,out] si pointer to spatial index structure
+ \param id item identifier
+ \param box pointer to item bounding box
- \param Map vector map
- \param Box bounding box
- \param[out] list output list, must be initialized
-
- \return number of isles
+ \return void
*/
-int
-Vect_select_isles_by_box(struct Map_info *Map, const BOUND_BOX * Box,
- struct ilist *list)
+void Vect_spatial_index_add_item(SPATIAL_INDEX * si, int id,
+ const BOUND_BOX * box)
{
- G_debug(3, "Vect_select_isles_by_box()");
- G_debug(3, "Box(N,S,E,W,T,B): %e, %e, %e, %e, %e, %e", Box->N, Box->S,
- Box->E, Box->W, Box->T, Box->B);
+ struct Rect rect;
- if (!(Map->plus.Spidx_built)) {
- G_debug(3, "Building spatial index.");
- Vect_build_sidx_from_topo(Map);
- }
+ G_debug(3, "Vect_spatial_index_add_item(): id = %d", id);
- dig_select_isles(&(Map->plus), Box, list);
- G_debug(3, " %d isles selected", list->n_values);
-
- return list->n_values;
+ rect.boundary[0] = box->W;
+ rect.boundary[1] = box->S;
+ rect.boundary[2] = box->B;
+ rect.boundary[3] = box->E;
+ rect.boundary[4] = box->N;
+ rect.boundary[5] = box->T;
+ RTreeInsertRect(&rect, id, si->si_tree);
}
/*!
- \brief Select nodes by box.
+ \brief Delete item from spatial index structure
- \param Map vector map
- \param Box bounding box
- \param[out] list output list, must be initialized
+ \param[in,out] si pointer to spatial index structure
+ \param id item identifier
- \return number of nodes
+ \return void
*/
-int
-Vect_select_nodes_by_box(struct Map_info *Map, const BOUND_BOX * Box,
- struct ilist *list)
+void Vect_spatial_index_del_item(SPATIAL_INDEX * si, int id,
+ const BOUND_BOX * box)
{
- struct Plus_head *plus;
+ int ret;
+ struct Rect rect;
- G_debug(3, "Vect_select_nodes_by_box()");
- G_debug(3, "Box(N,S,E,W,T,B): %e, %e, %e, %e, %e, %e", Box->N, Box->S,
- Box->E, Box->W, Box->T, Box->B);
+ G_debug(3, "Vect_spatial_index_del_item(): id = %d", id);
- plus = &(Map->plus);
+ rect.boundary[0] = box->W;
+ rect.boundary[1] = box->S;
+ rect.boundary[2] = box->B;
+ rect.boundary[3] = box->E;
+ rect.boundary[4] = box->N;
+ rect.boundary[5] = box->T;
- if (!(plus->Spidx_built)) {
- G_debug(3, "Building spatial index.");
- Vect_build_sidx_from_topo(Map);
- }
+ ret = RTreeDeleteRect(&rect, id, si->si_tree);
- list->n_values = 0;
-
- dig_select_nodes(plus, Box, list);
- G_debug(3, " %d nodes selected", list->n_values);
-
- return list->n_values;
+ if (ret)
+ G_fatal_error(_("Unable to delete item %d from spatial index"), id);
}
-/*!
- \brief Select lines by Polygon with optional isles.
-
- Polygons should be closed, i.e. first and last points must be identical.
-
- \param Map vector map
- \param Polygon outer ring
- \param nisles number of islands or 0
- \param Isles array of islands or NULL
- \param type line type
- \param[out] list output list, must be initialised
-
- \return number of lines
- */
-int
-Vect_select_lines_by_polygon(struct Map_info *Map, struct line_pnts *Polygon,
- int nisles, struct line_pnts **Isles, int type,
- struct ilist *List)
+/************************* SELECT BY BOX *********************************/
+/* This function is called by RTreeSearch() to add selected item to the list */
+static int _add_item(int id, struct ilist *list)
{
- int i;
- BOUND_BOX box;
- static struct line_pnts *LPoints = NULL;
- static struct ilist *LocList = NULL;
-
- /* TODO: this function was not tested with isles */
- G_debug(3, "Vect_select_lines_by_polygon() nisles = %d", nisles);
-
- List->n_values = 0;
- if (!LPoints)
- LPoints = Vect_new_line_struct();
- if (!LocList)
- LocList = Vect_new_list();
-
- /* Select first all lines by box */
- dig_line_box(Polygon, &box);
- Vect_select_lines_by_box(Map, &box, type, LocList);
- G_debug(3, " %d lines selected by box", LocList->n_values);
-
- /* Check all lines if intersect the polygon */
- for (i = 0; i < LocList->n_values; i++) {
- int j, line, intersect = 0;
-
- line = LocList->value[i];
- /* Read line points */
- Vect_read_line(Map, LPoints, NULL, line);
-
- /* Check if any of line vertices is within polygon */
- for (j = 0; j < LPoints->n_points; j++) {
- if (Vect_point_in_poly(LPoints->x[j], LPoints->y[j], Polygon) >= 1) { /* inside polygon */
- int k, inisle = 0;
-
- for (k = 0; k < nisles; k++) {
- if (Vect_point_in_poly(LPoints->x[j], LPoints->y[j], Isles[k]) >= 1) { /* in isle */
- inisle = 1;
- break;
- }
- }
-
- if (!inisle) { /* inside polygon, outside isles -> select */
- intersect = 1;
- break;
- }
- }
- }
- if (intersect) {
- dig_list_add(List, line);
- continue;
- }
-
- /* Check intersections of the line with area/isles boundary */
- /* Outer boundary */
- if (Vect_line_check_intersection(LPoints, Polygon, 0)) {
- dig_list_add(List, line);
- continue;
- }
-
- /* Islands */
- for (j = 0; j < nisles; j++) {
- if (Vect_line_check_intersection(LPoints, Isles[j], 0)) {
- intersect = 1;
- break;
- }
- }
- if (intersect) {
- dig_list_add(List, line);
- }
- }
-
- G_debug(4, " %d lines selected by polygon", List->n_values);
-
- return List->n_values;
+ dig_list_add(list, id);
+ return 1;
}
-
/*!
- \brief Select areas by Polygon with optional isles.
+ \brief Select items by bounding box to list
- Polygons should be closed, i.e. first and last points must be identical.
+ \param si pointer to spatial index structure
+ \param box bounding box
+ \param[out] list pointer to list where selected items are stored
- Warning : values in list may be duplicate!
-
- \param Map vector map
- \param Polygon outer ring
- \param nisles number of islands or 0
- \param Isles array of islands or NULL
- \param[out] list output list, must be initialised
-
- \return number of areas
+ \return number of selected items
*/
int
-Vect_select_areas_by_polygon(struct Map_info *Map, struct line_pnts *Polygon,
- int nisles, struct line_pnts **Isles,
- struct ilist *List)
+Vect_spatial_index_select(const SPATIAL_INDEX * si, const BOUND_BOX * box,
+ struct ilist *list)
{
- int i, area;
- static struct ilist *BoundList = NULL;
+ struct Rect rect;
- /* TODO: this function was not tested with isles */
- G_debug(3, "Vect_select_areas_by_polygon() nisles = %d", nisles);
+ G_debug(3, "Vect_spatial_index_select()");
- List->n_values = 0;
- if (!BoundList)
- BoundList = Vect_new_list();
+ Vect_reset_list(list);
- /* Select boundaries by polygon */
- Vect_select_lines_by_polygon(Map, Polygon, nisles, Isles, GV_BOUNDARY,
- BoundList);
+ rect.boundary[0] = box->W;
+ rect.boundary[1] = box->S;
+ rect.boundary[2] = box->B;
+ rect.boundary[3] = box->E;
+ rect.boundary[4] = box->N;
+ rect.boundary[5] = box->T;
+ RTreeSearch(si->si_tree, &rect, (void *)_add_item, list);
- /* Add areas on left/right side of selected boundaries */
- for (i = 0; i < BoundList->n_values; i++) {
- int line, left, right;
+ G_debug(3, " %d items selected", list->n_values);
- line = BoundList->value[i];
-
- Vect_get_line_areas(Map, line, &left, &right);
- G_debug(4, "boundary = %d left = %d right = %d", line, left, right);
-
- if (left > 0) {
- dig_list_add(List, left);
- }
- else if (left < 0) { /* island */
- area = Vect_get_isle_area(Map, abs(left));
- G_debug(4, " left island -> area = %d", area);
- if (area > 0)
- dig_list_add(List, area);
- }
-
- if (right > 0) {
- dig_list_add(List, right);
- }
- else if (right < 0) { /* island */
- area = Vect_get_isle_area(Map, abs(right));
- G_debug(4, " right island -> area = %d", area);
- if (area > 0)
- dig_list_add(List, area);
- }
- }
-
- /* But the Polygon may be completely inside the area (only one), in that case
- * we find the area by one polygon point and add it to the list */
- area = Vect_find_area(Map, Polygon->x[0], Polygon->y[0]);
- if (area > 0)
- dig_list_add(List, area);
-
- G_debug(3, " %d areas selected by polygon", List->n_values);
-
- return List->n_values;
+ return (list->n_values);
}
Modified: grass/trunk/lib/vector/Vlib/sindex.c
===================================================================
--- grass/trunk/lib/vector/Vlib/sindex.c 2009-07-13 12:26:17 UTC (rev 38384)
+++ grass/trunk/lib/vector/Vlib/sindex.c 2009-07-13 12:27:16 UTC (rev 38385)
@@ -1,7 +1,7 @@
/*!
\file Vlib/sindex.c
- \brief Vector library - spatial index
+ \brief Vector library - select vector features
Higher level functions for reading/writing/manipulating vectors.
@@ -11,263 +11,349 @@
(>=v2). Read the file COPYING that comes with GRASS for details.
\author Radim Blazek
- \author Martin Landa <landa.martin gmail.com> (storing sidx to file)
*/
#include <grass/config.h>
#include <stdlib.h>
-#include <string.h>
-#include <grass/glocale.h>
#include <grass/gis.h>
#include <grass/vector.h>
-#include <grass/glocale.h>
/*!
- \brief Initialize spatial index structure
+ \brief Select lines by box.
- \param si pointer to spatial index structure
+ Select lines whose boxes overlap specified box!!! It means that
+ selected line may or may not overlap the box.
- \return void
+ \param Map vector map
+ \param Box bounding box
+ \param type line type
+ \param[out] list output list, must be initialized
+
+ \return number of lines
*/
-void Vect_spatial_index_init(SPATIAL_INDEX *si)
+int
+Vect_select_lines_by_box(struct Map_info *Map, const BOUND_BOX * Box,
+ int type, struct ilist *list)
{
- G_debug(1, "Vect_spatial_index_init()");
+ int i, line, nlines;
+ struct Plus_head *plus;
+ P_LINE *Line;
+ static struct ilist *LocList = NULL;
- si->root = RTreeNewIndex();
-}
+ G_debug(3, "Vect_select_lines_by_box()");
+ G_debug(3, " Box(N,S,E,W,T,B): %e, %e, %e, %e, %e, %e", Box->N, Box->S,
+ Box->E, Box->W, Box->T, Box->B);
+ plus = &(Map->plus);
-/*!
- \brief Destroy existing spatial index
+ if (!(plus->Spidx_built)) {
+ G_debug(3, "Building spatial index.");
+ Vect_build_sidx_from_topo(Map);
+ }
- Vect_spatial_index_init() must be call before new use.
+ list->n_values = 0;
+ if (!LocList)
+ LocList = Vect_new_list();
- \param si pointer to spatial index structure
+ nlines = dig_select_lines(plus, Box, LocList);
+ G_debug(3, " %d lines selected (all types)", nlines);
- \return void
- */
-void Vect_spatial_index_destroy(SPATIAL_INDEX *si)
-{
- G_debug(1, "Vect_spatial_index_destroy()");
+ /* Remove lines of not requested types */
+ for (i = 0; i < nlines; i++) {
+ line = LocList->value[i];
+ if (plus->Line[line] == NULL)
+ continue; /* Should not happen */
+ Line = plus->Line[line];
+ if (!(Line->type & type))
+ continue;
+ dig_list_add(list, line);
+ }
- RTreeDestroyNode(si->root);
+ G_debug(3, " %d lines of requested type", list->n_values);
+
+ return list->n_values;
}
/*!
- \brief Add a new item to spatial index structure
+ \brief Select areas by box.
- \param[in,out] si pointer to spatial index structure
- \param id item identifier
- \param box pointer to item bounding box
+ Select areas whose boxes overlap specified box!!!
+ It means that selected area may or may not overlap the box.
- \return void
- */
-void Vect_spatial_index_add_item(SPATIAL_INDEX *si, int id, const BOUND_BOX *box)
-{
- struct Rect rect;
+ \param Map vector map
+ \param Box bounding box
+ \param[out] output list, must be initialized
- G_debug(3, "Vect_spatial_index_add_item(): id = %d", id);
-
- rect.boundary[0] = box->W;
- rect.boundary[1] = box->S;
- rect.boundary[2] = box->B;
- rect.boundary[3] = box->E;
- rect.boundary[4] = box->N;
- rect.boundary[5] = box->T;
- RTreeInsertRect(&rect, id, &(si->root), 0);
-}
-
-/*!
- \brief Delete item from spatial index structure
-
- \param[in,out] si pointer to spatial index structure
- \param id item identifier
-
- \return void
+ \return number of areas
*/
-void Vect_spatial_index_del_item(SPATIAL_INDEX * si, int id)
+int
+Vect_select_areas_by_box(struct Map_info *Map, const BOUND_BOX * Box,
+ struct ilist *list)
{
- int ret;
- struct Rect rect;
+ int i;
+ const char *dstr;
+ int debug_level;
- G_debug(3, "Vect_spatial_index_del_item(): id = %d", id);
+ dstr = G__getenv("DEBUG");
- /* TODO */
- G_fatal_error("Vect_spatial_index_del_item() %s", _("not implemented"));
+ if (dstr != NULL)
+ debug_level = atoi(dstr);
+ else
+ debug_level = 0;
- /* Bounding box of item would be needed, which is not stored in si. */
+ G_debug(3, "Vect_select_areas_by_box()");
+ G_debug(3, "Box(N,S,E,W,T,B): %e, %e, %e, %e, %e, %e", Box->N, Box->S,
+ Box->E, Box->W, Box->T, Box->B);
- /*
- rect.boundary[0] = ; rect.boundary[1] = ; rect.boundary[2] = ;
- rect.boundary[3] = ; rect.boundary[4] = ; rect.boundary[5] = ;
- */
+ if (!(Map->plus.Spidx_built)) {
+ G_debug(3, "Building spatial index.");
+ Vect_build_sidx_from_topo(Map);
+ }
- ret = RTreeDeleteRect(&rect, id, &(si->root));
+ dig_select_areas(&(Map->plus), Box, list);
+ G_debug(3, " %d areas selected", list->n_values);
+ /* avoid loop when not debugging */
+ if (debug_level > 2) {
+ for (i = 0; i < list->n_values; i++) {
+ G_debug(3, " area = %d pointer to area structure = %lx",
+ list->value[i],
+ (unsigned long)Map->plus.Area[list->value[i]]);
+ }
+ }
- if (ret)
- G_fatal_error(_("Unable to delete item %d from spatial index"), id);
+ return list->n_values;
}
+
/*!
- \brief Create spatial index if necessary.
+ \brief Select isles by box.
- To be used in modules.
- Map must be opened on level 2.
+ Select isles whose boxes overlap specified box!!!
+ It means that selected isle may or may not overlap the box.
- \param[in,out] Map pointer to vector map
-
- \return 0 OK
- \return 1 error
+ \param Map vector map
+ \param Box bounding box
+ \param[out] list output list, must be initialized
+
+ \return number of isles
*/
-int Vect_build_spatial_index(struct Map_info *Map)
+int
+Vect_select_isles_by_box(struct Map_info *Map, const BOUND_BOX * Box,
+ struct ilist *list)
{
- if (Map->level < 2) {
- G_fatal_error(_("Unable to build spatial index from topology, "
- "vector map is not opened at topology level 2"));
- }
+ G_debug(3, "Vect_select_isles_by_box()");
+ G_debug(3, "Box(N,S,E,W,T,B): %e, %e, %e, %e, %e, %e", Box->N, Box->S,
+ Box->E, Box->W, Box->T, Box->B);
+
if (!(Map->plus.Spidx_built)) {
- return (Vect_build_sidx_from_topo(Map));
+ G_debug(3, "Building spatial index.");
+ Vect_build_sidx_from_topo(Map);
}
- return 0;
+
+ dig_select_isles(&(Map->plus), Box, list);
+ G_debug(3, " %d isles selected", list->n_values);
+
+ return list->n_values;
}
/*!
- \brief Create spatial index from topology if necessary
+ \brief Select nodes by box.
- \param Map pointer to vector map
+ \param Map vector map
+ \param Box bounding box
+ \param[out] list output list, must be initialized
- \return 0 OK
- \return 1 error
+ \return number of nodes
*/
-int Vect_build_sidx_from_topo(struct Map_info *Map)
+int
+Vect_select_nodes_by_box(struct Map_info *Map, const BOUND_BOX * Box,
+ struct ilist *list)
{
- int i, total, done;
struct Plus_head *plus;
- BOUND_BOX box;
- P_LINE *Line;
- P_NODE *Node;
- P_AREA *Area;
- P_ISLE *Isle;
- G_debug(3, "Vect_build_sidx_from_topo()");
+ G_debug(3, "Vect_select_nodes_by_box()");
+ G_debug(3, "Box(N,S,E,W,T,B): %e, %e, %e, %e, %e, %e", Box->N, Box->S,
+ Box->E, Box->W, Box->T, Box->B);
plus = &(Map->plus);
- dig_spidx_init(plus);
+ if (!(plus->Spidx_built)) {
+ G_debug(3, "Building spatial index.");
+ Vect_build_sidx_from_topo(Map);
+ }
- total = plus->n_nodes + plus->n_lines + plus->n_areas + plus->n_isles;
+ list->n_values = 0;
- /* Nodes */
- for (i = 1; i <= plus->n_nodes; i++) {
- G_percent(i, total, 3);
+ dig_select_nodes(plus, Box, list);
+ G_debug(3, " %d nodes selected", list->n_values);
- Node = plus->Node[i];
- if (!Node)
- G_fatal_error(_("BUG (Vect_build_sidx_from_topo): node does not exist"));
+ return list->n_values;
+}
- dig_spidx_add_node(plus, i, Node->x, Node->y, Node->z);
- }
+/*!
+ \brief Select lines by Polygon with optional isles.
- /* Lines */
- done = plus->n_nodes;
- for (i = 1; i <= plus->n_lines; i++) {
- G_percent(done + i, total, 3);
+ Polygons should be closed, i.e. first and last points must be identical.
- Line = plus->Line[i];
- if (!Line)
- G_fatal_error(_("BUG (Vect_build_sidx_from_topo): line does not exist"));
+ \param Map vector map
+ \param Polygon outer ring
+ \param nisles number of islands or 0
+ \param Isles array of islands or NULL
+ \param type line type
+ \param[out] list output list, must be initialised
- box.N = Line->N;
- box.S = Line->S;
- box.E = Line->E;
- box.W = Line->W;
- box.T = Line->T;
- box.B = Line->B;
+ \return number of lines
+ */
+int
+Vect_select_lines_by_polygon(struct Map_info *Map, struct line_pnts *Polygon,
+ int nisles, struct line_pnts **Isles, int type,
+ struct ilist *List)
+{
+ int i;
+ BOUND_BOX box;
+ static struct line_pnts *LPoints = NULL;
+ static struct ilist *LocList = NULL;
- dig_spidx_add_line(plus, i, &box);
- }
+ /* TODO: this function was not tested with isles */
+ G_debug(3, "Vect_select_lines_by_polygon() nisles = %d", nisles);
- /* Areas */
- done += plus->n_lines;
- for (i = 1; i <= plus->n_areas; i++) {
- G_percent(done + i, total, 3);
+ List->n_values = 0;
+ if (!LPoints)
+ LPoints = Vect_new_line_struct();
+ if (!LocList)
+ LocList = Vect_new_list();
- Area = plus->Area[i];
- if (!Area)
- G_fatal_error(_("BUG (Vect_build_sidx_from_topo): area does not exist"));
+ /* Select first all lines by box */
+ dig_line_box(Polygon, &box);
+ Vect_select_lines_by_box(Map, &box, type, LocList);
+ G_debug(3, " %d lines selected by box", LocList->n_values);
- box.N = Area->N;
- box.S = Area->S;
- box.E = Area->E;
- box.W = Area->W;
- box.T = Area->T;
- box.B = Area->B;
+ /* Check all lines if intersect the polygon */
+ for (i = 0; i < LocList->n_values; i++) {
+ int j, line, intersect = 0;
- dig_spidx_add_area(plus, i, &box);
- }
+ line = LocList->value[i];
+ /* Read line points */
+ Vect_read_line(Map, LPoints, NULL, line);
- /* Isles */
- done += plus->n_areas;
- for (i = 1; i <= plus->n_isles; i++) {
- G_percent(done + i, total, 3);
+ /* Check if any of line vertices is within polygon */
+ for (j = 0; j < LPoints->n_points; j++) {
+ if (Vect_point_in_poly(LPoints->x[j], LPoints->y[j], Polygon) >= 1) { /* inside polygon */
+ int k, inisle = 0;
- Isle = plus->Isle[i];
- if (!Isle)
- G_fatal_error(_("BUG (Vect_build_sidx_from_topo): isle does not exist"));
+ for (k = 0; k < nisles; k++) {
+ if (Vect_point_in_poly(LPoints->x[j], LPoints->y[j], Isles[k]) >= 1) { /* in isle */
+ inisle = 1;
+ break;
+ }
+ }
- box.N = Isle->N;
- box.S = Isle->S;
- box.E = Isle->E;
- box.W = Isle->W;
- box.T = Isle->T;
- box.B = Isle->B;
+ if (!inisle) { /* inside polygon, outside isles -> select */
+ intersect = 1;
+ break;
+ }
+ }
+ }
+ if (intersect) {
+ dig_list_add(List, line);
+ continue;
+ }
- dig_spidx_add_isle(plus, i, &box);
+ /* Check intersections of the line with area/isles boundary */
+ /* Outer boundary */
+ if (Vect_line_check_intersection(LPoints, Polygon, 0)) {
+ dig_list_add(List, line);
+ continue;
+ }
+
+ /* Islands */
+ for (j = 0; j < nisles; j++) {
+ if (Vect_line_check_intersection(LPoints, Isles[j], 0)) {
+ intersect = 1;
+ break;
+ }
+ }
+ if (intersect) {
+ dig_list_add(List, line);
+ }
}
- Map->plus.Spidx_built = 1;
+ G_debug(4, " %d lines selected by polygon", List->n_values);
- G_debug(3, "Spatial index was built");
-
- return 0;
+ return List->n_values;
}
-/************************* SELECT BY BOX *********************************/
-/* This function is called by RTreeSearch() to add selected item to the list */
-static int _add_item(int id, struct ilist *list)
-{
- dig_list_add(list, id);
- return 1;
-}
/*!
- \brief Select items by bounding box to list
+ \brief Select areas by Polygon with optional isles.
- \param si pointer to spatial index structure
- \param box bounding box
- \param[out] list pointer to list where selected items are stored
+ Polygons should be closed, i.e. first and last points must be identical.
- \return number of selected items
+ Warning : values in list may be duplicate!
+
+ \param Map vector map
+ \param Polygon outer ring
+ \param nisles number of islands or 0
+ \param Isles array of islands or NULL
+ \param[out] list output list, must be initialised
+
+ \return number of areas
*/
int
-Vect_spatial_index_select(const SPATIAL_INDEX * si, const BOUND_BOX * box,
- struct ilist *list)
+Vect_select_areas_by_polygon(struct Map_info *Map, struct line_pnts *Polygon,
+ int nisles, struct line_pnts **Isles,
+ struct ilist *List)
{
- struct Rect rect;
+ int i, area;
+ static struct ilist *BoundList = NULL;
- G_debug(3, "Vect_spatial_index_select()");
+ /* TODO: this function was not tested with isles */
+ G_debug(3, "Vect_select_areas_by_polygon() nisles = %d", nisles);
- Vect_reset_list(list);
+ List->n_values = 0;
+ if (!BoundList)
+ BoundList = Vect_new_list();
- rect.boundary[0] = box->W;
- rect.boundary[1] = box->S;
- rect.boundary[2] = box->B;
- rect.boundary[3] = box->E;
- rect.boundary[4] = box->N;
- rect.boundary[5] = box->T;
- RTreeSearch(si->root, &rect, (void *)_add_item, list);
+ /* Select boundaries by polygon */
+ Vect_select_lines_by_polygon(Map, Polygon, nisles, Isles, GV_BOUNDARY,
+ BoundList);
- G_debug(3, " %d items selected", list->n_values);
+ /* Add areas on left/right side of selected boundaries */
+ for (i = 0; i < BoundList->n_values; i++) {
+ int line, left, right;
- return (list->n_values);
+ line = BoundList->value[i];
+
+ Vect_get_line_areas(Map, line, &left, &right);
+ G_debug(4, "boundary = %d left = %d right = %d", line, left, right);
+
+ if (left > 0) {
+ dig_list_add(List, left);
+ }
+ else if (left < 0) { /* island */
+ area = Vect_get_isle_area(Map, abs(left));
+ G_debug(4, " left island -> area = %d", area);
+ if (area > 0)
+ dig_list_add(List, area);
+ }
+
+ if (right > 0) {
+ dig_list_add(List, right);
+ }
+ else if (right < 0) { /* island */
+ area = Vect_get_isle_area(Map, abs(right));
+ G_debug(4, " right island -> area = %d", area);
+ if (area > 0)
+ dig_list_add(List, area);
+ }
+ }
+
+ /* But the Polygon may be completely inside the area (only one), in that case
+ * we find the area by one polygon point and add it to the list */
+ area = Vect_find_area(Map, Polygon->x[0], Polygon->y[0]);
+ if (area > 0)
+ dig_list_add(List, area);
+
+ G_debug(3, " %d areas selected by polygon", List->n_values);
+
+ return List->n_values;
}
-
More information about the grass-commit
mailing list