Vlib new spatial index

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

lib/vector/Vlib/build.c

@@ -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 @@ int Vect_build_partial(struct Map_info *Map, int build)
     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 @@ int Vect_build_partial(struct Map_info *Map, int build)
 
 	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 @@ int Vect_build_partial(struct Map_info *Map, int build)
 
 	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 @@ int Vect_topo_dump(const struct Map_info *Map, FILE * out)
 	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,43 +383,179 @@ int Vect_topo_dump(const struct Map_info *Map, FILE * out)
 }
 
 /*!
-   \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) {
+
+	/*  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;
+	}
 
-    if (0 > dig_write_spidx(&fp, plus)) {
-	G_warning(_("Error writing out spatial index file"));
-	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 @@ int Vect_save_spatial_index(struct Map_info *Map)
    \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);

lib/vector/Vlib/build_nat.c

@@ -73,16 +73,17 @@ int Vect_build_line_area(struct Map_info *Map, int iline, int side)
 	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 @@ int Vect_isle_find_area(struct Map_info *Map, int isle)
 		    /* 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 @@ int Vect_isle_find_area(struct Map_info *Map, int isle)
 
 		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 @@ int Vect_attach_centroids(struct Map_info *Map, const BOUND_BOX * box)
 	 * 3) it's faster
 	 */
 	if (Line->left > 0)
-	    continue; 
+	    continue;
 
 	orig_area = Line->left;
 
@@ -539,7 +540,7 @@ int Vect_build_nat(struct Map_info *Map, int build)
 
 	    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 @@ int Vect_build_nat(struct Map_info *Map, int build)
 		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);

lib/vector/Vlib/close.c

@@ -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 @@ int Vect_close(struct Map_info *Map)
 
 	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 @@ int Vect_close(struct Map_info *Map)
 	    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));

lib/vector/Vlib/intersect.c

@@ -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>
-                 \  /    \  /  \  /
-                  \/      \/    \/
-                  /\             \
-                 /  \             \
-           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>
+   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 @@ Vect_line_intersection(struct line_pnts *APoints,
     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 @@ Vect_line_intersection(struct line_pnts *APoints,
      *  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 @@ Vect_line_intersection(struct line_pnts *APoints,
 	    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 @@ Vect_line_intersection(struct line_pnts *APoints,
 	    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 @@ static int find_cross(int id, int *arg)
 
     if (!IPnts)
 	IPnts = Vect_new_line_struct();
-    
+
     switch (ret) {
     case 0:
     case 5:
@@ -1132,7 +1134,7 @@ Vect_line_check_intersection(struct line_pnts *APoints,
     int i;
     double dist, rethresh;
     struct Rect rect;
-    struct Node *RTree;
+    struct RTree *MyRTree;
 
     rethresh = 0.000001;	/* TODO */
     APnts = APoints;
@@ -1212,7 +1214,7 @@ Vect_line_check_intersection(struct line_pnts *APoints,
      *  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 @@ Vect_line_check_intersection(struct line_pnts *APoints,
 	    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 @@ Vect_line_check_intersection(struct line_pnts *APoints,
 	    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 @@ Vect_line_check_intersection(struct line_pnts *APoints,
     }
 
     /* Free RTree */
-    RTreeDestroyNode(RTree);
+    RTreeFreeIndex(MyRTree);
 
     return cross_found;
 }

lib/vector/Vlib/open.c

@@ -212,7 +212,8 @@ Vect__open_old(struct Map_info *Map, const char *name, const char *mapset,
 		_("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 @@ Vect__open_old(struct Map_info *Map, const char *name, const char *mapset,
 	/* 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 @@ Vect__open_old(struct Map_info *Map, const char *name, const char *mapset,
 	    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 @@ Vect__open_old(struct Map_info *Map, const char *name, const char *mapset,
 	    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 @@ Vect_open_update(struct Map_info *Map, const char *name, const char *mapset)
 	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 @@ Vect_open_update_head(struct Map_info *Map, const char *name,
 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 @@ int Vect_open_new(struct Map_info *Map, const char *name, int with_z)
 
     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 @@ int Vect_coor_info(const struct Map_info *Map, struct Coor_info *Info)
 	    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 @@ int Vect_coor_info(const struct Map_info *Map, struct Coor_info *Info)
 	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 @@ int Vect_open_topo(struct Map_info *Map, int head_only)
 	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 @@ int Vect_open_topo(struct Map_info *Map, int head_only)
  * \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);
 
-    /* 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 (stat(file_path, &info) != 0)	/* does not exist */
+	    return 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;
-       }
-     */
+	Map->plus.spidx_fp.file =
+	    G_fopen_old(buf, GV_SIDX_ELEMENT, Map->mapset);
 
-    /* load file to the memory */
-    /* dig_file_load ( &fp); */
+	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;
+	}
 
-    /* load topo to memory */
-    dig_spidx_init(Plus);
-    dig_read_spidx(&fp, Plus);
+	/* get coor info */
+	Vect_coor_info(Map, &CInfo);
 
-    fclose(fp.file);
-    /* dig_file_free ( &fp); */
+	/* 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;
+	}
+    }
+
+    if (mode) {
+	/* open new spatial index */
+	Map->plus.Spidx_new = 1;
+
+	dig_spidx_init(Plus);
+
+	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;
+	    }
+	}
+    }
+
+    Map->plus.Spidx_built = 1;
 
     return 0;
 }

lib/vector/Vlib/select.c

@@ -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(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;
+    G_debug(1, "Vect_spatial_index_init()");
+
+    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;
-
-    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;
-}
+    G_debug(1, "Vect_spatial_index_destroy()");
 
+    RTreeFreeIndex(si->si_tree);
+}
 
 /*!
-   \brief Select isles by box.
-
-   Select isles whose boxes overlap specified box!!!
-   It means that selected isle may or may not overlap the box.
+   \brief Add a new item to 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
+   \param box pointer to item bounding box
 
-   \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);
+    if (ret)
+	G_fatal_error(_("Unable to delete item %d from spatial index"), id);
+}
 
-    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 lines 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
 
-   \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
+   \return number of selected items
  */
 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)
+Vect_spatial_index_select(const SPATIAL_INDEX * si, const BOUND_BOX * box,
+			  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;
-}
+    struct Rect rect;
 
+    G_debug(3, "Vect_spatial_index_select()");
 
-/*!
-   \brief Select areas by Polygon with optional isles. 
+    Vect_reset_list(list);
 
-   Polygons should be closed, i.e. first and last points must be identical.
+    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);
 
-   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
+    G_debug(3, "  %d items selected", list->n_values);
 
-   \return number of areas
- */
-int
-Vect_select_areas_by_polygon(struct Map_info *Map, struct line_pnts *Polygon,
-			     int nisles, struct line_pnts **Isles,
-			     struct ilist *List)
-{
-    int i, area;
-    static struct ilist *BoundList = NULL;
-
-    /* TODO: this function was not tested with isles */
-    G_debug(3, "Vect_select_areas_by_polygon() nisles = %d", nisles);
-
-    List->n_values = 0;
-    if (!BoundList)
-	BoundList = Vect_new_list();
-
-    /* Select boundaries by polygon */
-    Vect_select_lines_by_polygon(Map, Polygon, nisles, Isles, GV_BOUNDARY,
-				 BoundList);
-
-    /* Add areas on left/right side of selected boundaries */
-    for (i = 0; i < BoundList->n_values; i++) {
-	int line, left, right;
-
-	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);
 }

lib/vector/Vlib/sindex.c

@@ -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
- */
-void Vect_spatial_index_init(SPATIAL_INDEX *si)
-{
-    G_debug(1, "Vect_spatial_index_init()");
-
-    si->root = RTreeNewIndex();
-}
-
-/*!
-   \brief Destroy existing spatial index
-
-   Vect_spatial_index_init() must be call before new use.
-
-   \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 void
+   \return number of lines
  */
-void Vect_spatial_index_destroy(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_destroy()");
-
-    RTreeDestroyNode(si->root);
-}
+    int i, line, nlines;
+    struct Plus_head *plus;
+    P_LINE *Line;
+    static struct ilist *LocList = NULL;
 
-/*!
-   \brief Add a new item to spatial index structure
+    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);
 
-   \param[in,out] si pointer to spatial index structure
-   \param id item identifier
-   \param box pointer to item bounding box
+    if (!(plus->Spidx_built)) {
+	G_debug(3, "Building spatial index.");
+	Vect_build_sidx_from_topo(Map);
+    }
 
-   \return void
- */
-void Vect_spatial_index_add_item(SPATIAL_INDEX *si, int id, const BOUND_BOX *box)
-{
-    struct Rect rect;
+    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, "Vect_spatial_index_add_item(): id = %d", id);
+    G_debug(3, "  %d lines of requested type", 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->root), 0);
+    return list->n_values;
 }
 
 /*!
-   \brief Delete item from spatial index structure
+   \brief Select areas by box.
 
-   \param[in,out] si pointer to spatial index structure
-   \param id item identifier
+   Select areas whose boxes overlap specified box!!!
+   It means that selected area may or may not overlap the box.
 
-   \return void
+   \param Map vector map
+   \param Box bounding box
+   \param[out] output list, must be initialized
+
+   \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);
-
-    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);
+    if (!(plus->Spidx_built)) {
+	G_debug(3, "Building spatial index.");
+	Vect_build_sidx_from_topo(Map);
     }
 
-    /* Lines */
-    done = plus->n_nodes;
-    for (i = 1; i <= plus->n_lines; i++) {
-	G_percent(done + i, total, 3);
+    list->n_values = 0;
 
-	Line = plus->Line[i];
-	if (!Line)
-	    G_fatal_error(_("BUG (Vect_build_sidx_from_topo): line does not exist"));
+    dig_select_nodes(plus, Box, list);
+    G_debug(3, "  %d nodes selected", list->n_values);
 
-	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);
-    }
+    return list->n_values;
+}
 
-    /* Isles */
-    done += plus->n_areas;
-    for (i = 1; i <= plus->n_isles; i++) {
-	G_percent(done + i, total, 3);
+/*!
+   \brief Select lines by Polygon with optional isles. 
 
-	Isle = plus->Isle[i];
-	if (!Isle)
-	    G_fatal_error(_("BUG (Vect_build_sidx_from_topo): isle does not exist"));
+   Polygons should be closed, i.e. first and last points must be identical.
 
-	box.N = Isle->N;
-	box.S = Isle->S;
-	box.E = Isle->E;
-	box.W = Isle->W;
-	box.T = Isle->T;
-	box.B = Isle->B;
+   \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
 
-	dig_spidx_add_isle(plus, i, &box);
+   \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;
+
+    /* 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);
+	}
     }
 
-    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. 
+
+   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!
 
-   \return number of selected items
+   \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;
-
-    G_debug(3, "Vect_spatial_index_select()");
-
-    Vect_reset_list(list);
+    int i, area;
+    static struct ilist *BoundList = NULL;
+
+    /* TODO: this function was not tested with isles */
+    G_debug(3, "Vect_select_areas_by_polygon() nisles = %d", nisles);
+
+    List->n_values = 0;
+    if (!BoundList)
+	BoundList = Vect_new_list();
+
+    /* Select boundaries by polygon */
+    Vect_select_lines_by_polygon(Map, Polygon, nisles, Isles, GV_BOUNDARY,
+				 BoundList);
+
+    /* Add areas on left/right side of selected boundaries */
+    for (i = 0; i < BoundList->n_values; i++) {
+	int line, left, right;
+
+	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);
+	}
+    }
 
-    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);
+    /* 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 items selected", list->n_values);
+    G_debug(3, "  %d areas selected by polygon", List->n_values);
 
-    return (list->n_values);
+    return List->n_values;
 }
-