GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
sindex.c
Go to the documentation of this file.
1 /*!
2  \file lib/vector/Vlib/sindex.c
3 
4  \brief Vector library - select vector features
5 
6  Higher level functions for reading/writing/manipulating vectors.
7 
8  (C) 2001-2011 by the GRASS Development Team
9 
10  This program is free software under the GNU General Public License
11  (>=v2). Read the file COPYING that comes with GRASS for details.
12 
13  \author Radim Blazek, Markus Metz
14  */
15 
16 #include <stdlib.h>
17 #include <grass/vector.h>
18 
19 
20 /*!
21  \brief Select lines with bounding boxes by box.
22 
23  Select lines whose boxes overlap specified box!!! It means that
24  selected line may or may not overlap the box.
25 
26  \param Map vector map
27  \param Box bounding box
28  \param type line type
29  \param[out] list output list, must be initialized
30 
31  \return number of lines
32  */
33 int
34 Vect_select_lines_by_box(struct Map_info *Map, const struct bound_box *Box,
35  int type, struct boxlist *list)
36 {
37  int i, line, nlines, ntypes, mtype;
38  struct Plus_head *plus;
39  struct P_line *Line;
40  static struct boxlist *LocList = NULL;
41 
42  G_debug(3, "Vect_select_lines_by_box()");
43  G_debug(3, " Box(N,S,E,W,T,B): %e, %e, %e, %e, %e, %e", Box->N, Box->S,
44  Box->E, Box->W, Box->T, Box->B);
45  plus = &(Map->plus);
46 
47  Vect_reset_boxlist(list);
48 
49  ntypes = mtype = 0;
50  /* count the number of different primitives in Map */
51  if (plus->n_plines != 0) {
52  ntypes++;
53  mtype |= GV_POINT;
54  }
55  if (plus->n_llines != 0) {
56  ntypes++;
57  mtype |= GV_LINE;
58  }
59  if (plus->n_blines != 0) {
60  ntypes++;
61  mtype |= GV_BOUNDARY;
62  }
63  if (plus->n_clines != 0) {
64  ntypes++;
65  mtype |= GV_CENTROID;
66  }
67  if (plus->n_flines != 0) {
68  ntypes++;
69  mtype |= GV_FACE;
70  }
71  if (plus->n_klines != 0) {
72  ntypes++;
73  mtype |= GV_KERNEL;
74  }
75 
76  if (ntypes == 1) {
77  /* there is only one type in Map */
78  if (mtype & type)
79  return dig_select_lines(plus, Box, list);
80  return 0;
81  }
82 
83  if (ntypes == 0)
84  /* empty vector */
85  return 0;
86 
87  if (!LocList) {
88  LocList = (struct boxlist *)G_malloc(sizeof(struct boxlist));
89  dig_init_boxlist(LocList, 1);
90  }
91 
92  nlines = dig_select_lines(plus, Box, LocList);
93  G_debug(3, " %d lines selected (all types)", nlines);
94 
95  /* Remove lines of not requested types */
96  for (i = 0; i < nlines; i++) {
97  line = LocList->id[i];
98  if (plus->Line[line] == NULL)
99  continue; /* Should not happen */
100  Line = plus->Line[line];
101  if (!(Line->type & type))
102  continue;
103  dig_boxlist_add(list, line, &LocList->box[i]);
104  }
105 
106  G_debug(3, " %d lines of requested type", list->n_values);
107 
108  return list->n_values;
109 }
110 
111 
112 /*!
113  \brief Select areas with bounding boxes by box.
114 
115  Select areas whose boxes overlap specified box!!!
116  It means that selected area may or may not overlap the box.
117 
118  \param Map vector map
119  \param Box bounding box
120  \param[out] output list, must be initialized
121 
122  \return number of areas
123  */
124 int
125 Vect_select_areas_by_box(struct Map_info *Map, const struct bound_box * Box,
126  struct boxlist *list)
127 {
128  int i;
129  static int debug_level = -1;
130 
131  if (debug_level == -1) {
132  const char *dstr = G_getenv_nofatal("DEBUG");
133 
134  if (dstr != NULL)
135  debug_level = atoi(dstr);
136  else
137  debug_level = 0;
138  }
139 
140  G_debug(3, "Vect_select_areas_by_box()");
141  G_debug(3, "Box(N,S,E,W,T,B): %e, %e, %e, %e, %e, %e", Box->N, Box->S,
142  Box->E, Box->W, Box->T, Box->B);
143 
144  dig_select_areas(&(Map->plus), Box, list);
145  G_debug(3, " %d areas selected", list->n_values);
146  /* avoid loop when not debugging */
147  if (debug_level > 2) {
148  for (i = 0; i < list->n_values; i++) {
149  G_debug(3, " area = %d pointer to area structure = %p",
150  list->id[i], (void *)Map->plus.Area[list->id[i]]);
151  }
152  }
153 
154  return list->n_values;
155 }
156 
157 
158 /*!
159  \brief Select isles with bounding boxes by box.
160 
161  Select isles whose boxes overlap specified box!!!
162  It means that selected isle may or may not overlap the box.
163 
164  \param Map vector map
165  \param Box bounding box
166  \param[out] list output list, must be initialized
167 
168  \return number of isles
169  */
170 int
171 Vect_select_isles_by_box(struct Map_info *Map, const struct bound_box * Box,
172  struct boxlist *list)
173 {
174  G_debug(3, "Vect_select_isles_by_box()");
175  G_debug(3, "Box(N,S,E,W,T,B): %e, %e, %e, %e, %e, %e", Box->N, Box->S,
176  Box->E, Box->W, Box->T, Box->B);
177 
178  dig_select_isles(&(Map->plus), Box, list);
179  G_debug(3, " %d isles selected", list->n_values);
180 
181  return list->n_values;
182 }
183 
184 /*!
185  \brief Select nodes by box.
186 
187  \param Map vector map
188  \param Box bounding box
189  \param[out] list output list, must be initialized
190 
191  \return number of nodes
192  */
193 int
194 Vect_select_nodes_by_box(struct Map_info *Map, const struct bound_box * Box,
195  struct ilist *list)
196 {
197  struct Plus_head *plus;
198 
199  G_debug(3, "Vect_select_nodes_by_box()");
200  G_debug(3, "Box(N,S,E,W,T,B): %e, %e, %e, %e, %e, %e", Box->N, Box->S,
201  Box->E, Box->W, Box->T, Box->B);
202 
203  plus = &(Map->plus);
204 
205  Vect_reset_list(list);
206 
207  dig_select_nodes(plus, Box, list);
208  G_debug(3, " %d nodes selected", list->n_values);
209 
210  return list->n_values;
211 }
212 
213 /*!
214  \brief Select lines by Polygon with optional isles.
215 
216  Polygons should be closed, i.e. first and last points must be identical.
217 
218  \param Map vector map
219  \param Polygon outer ring
220  \param nisles number of islands or 0
221  \param Isles array of islands or NULL
222  \param type line type
223  \param[out] list output list, must be initialised
224 
225  \return number of lines
226  */
227 int
228 Vect_select_lines_by_polygon(struct Map_info *Map, struct line_pnts *Polygon,
229  int nisles, struct line_pnts **Isles, int type,
230  struct ilist *List)
231 {
232  int i;
233  struct bound_box box;
234  static struct line_pnts *LPoints = NULL;
235  static struct boxlist *LocList = NULL;
236 
237  /* TODO: this function was not tested with isles */
238  G_debug(3, "Vect_select_lines_by_polygon() nisles = %d", nisles);
239 
240  Vect_reset_list(List);
241  if (!LPoints)
242  LPoints = Vect_new_line_struct();
243  if (!LocList) {
244  LocList = Vect_new_boxlist(0);
245  }
246 
247  /* Select first all lines by box */
248  dig_line_box(Polygon, &box);
249  box.T = PORT_DOUBLE_MAX;
250  box.B = -PORT_DOUBLE_MAX;
251  Vect_select_lines_by_box(Map, &box, type, LocList);
252  G_debug(3, " %d lines selected by box", LocList->n_values);
253 
254  /* Check all lines if intersect the polygon */
255  for (i = 0; i < LocList->n_values; i++) {
256  int j, line, intersect = 0;
257 
258  line = LocList->id[i];
259  /* Read line points */
260  Vect_read_line(Map, LPoints, NULL, line);
261 
262  /* Check if any of line vertices is within polygon */
263  for (j = 0; j < LPoints->n_points; j++) {
264  if (Vect_point_in_poly(LPoints->x[j], LPoints->y[j], Polygon) >= 1) { /* inside polygon */
265  int k, inisle = 0;
266 
267  for (k = 0; k < nisles; k++) {
268  if (Vect_point_in_poly(LPoints->x[j], LPoints->y[j], Isles[k]) >= 1) { /* in isle */
269  inisle = 1;
270  break;
271  }
272  }
273 
274  if (!inisle) { /* inside polygon, outside isles -> select */
275  intersect = 1;
276  break;
277  }
278  }
279  }
280  if (intersect) {
281  Vect_list_append(List, line);
282  continue;
283  }
284 
285  /* Check intersections of the line with area/isles boundary */
286  /* Outer boundary */
287  if (Vect_line_check_intersection(LPoints, Polygon, 0)) {
288  Vect_list_append(List, line);
289  continue;
290  }
291 
292  /* Islands */
293  for (j = 0; j < nisles; j++) {
294  if (Vect_line_check_intersection(LPoints, Isles[j], 0)) {
295  intersect = 1;
296  break;
297  }
298  }
299  if (intersect) {
300  Vect_list_append(List, line);
301  }
302  }
303 
304  G_debug(4, " %d lines selected by polygon", List->n_values);
305 
306  return List->n_values;
307 }
308 
309 
310 /*!
311  \brief Select areas by Polygon with optional isles.
312 
313  Polygons should be closed, i.e. first and last points must be identical.
314 
315  \param Map vector map
316  \param Polygon outer ring
317  \param nisles number of islands or 0
318  \param Isles array of islands or NULL
319  \param[out] list output list, must be initialised
320 
321  \return number of areas
322  */
323 int
324 Vect_select_areas_by_polygon(struct Map_info *Map, struct line_pnts *Polygon,
325  int nisles, struct line_pnts **Isles,
326  struct ilist *List)
327 {
328  int i, area;
329  static struct ilist *BoundList = NULL;
330 
331  /* TODO: this function was not tested with isles */
332  G_debug(3, "Vect_select_areas_by_polygon() nisles = %d", nisles);
333 
334  Vect_reset_list(List);
335  if (!BoundList)
336  BoundList = Vect_new_list();
337 
338  /* Select boundaries by polygon */
339  Vect_select_lines_by_polygon(Map, Polygon, nisles, Isles, GV_BOUNDARY,
340  BoundList);
341 
342  /* Add areas on left/right side of selected boundaries */
343  for (i = 0; i < BoundList->n_values; i++) {
344  int line, left, right;
345 
346  line = BoundList->value[i];
347 
348  Vect_get_line_areas(Map, line, &left, &right);
349  G_debug(4, "boundary = %d left = %d right = %d", line, left, right);
350 
351  if (left > 0) {
352  Vect_list_append(List, left);
353  }
354  else if (left < 0) { /* island */
355  area = Vect_get_isle_area(Map, abs(left));
356  G_debug(4, " left island -> area = %d", area);
357  if (area > 0)
358  Vect_list_append(List, area);
359  }
360 
361  if (right > 0) {
362  Vect_list_append(List, right);
363  }
364  else if (right < 0) { /* island */
365  area = Vect_get_isle_area(Map, abs(right));
366  G_debug(4, " right island -> area = %d", area);
367  if (area > 0)
368  Vect_list_append(List, area);
369  }
370  }
371 
372  /* But the Polygon may be completely inside the area (only one), in that case
373  * we find the area by one polygon point and add it to the list */
374  area = Vect_find_area(Map, Polygon->x[0], Polygon->y[0]);
375  if (area > 0)
376  Vect_list_append(List, area);
377 
378  G_debug(3, " %d areas selected by polygon", List->n_values);
379 
380  return List->n_values;
381 }
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 lines by Polygon with optional isles.
Definition: sindex.c:228
#define G_malloc(n)
Definition: defs/gis.h:112
int Vect_select_areas_by_box(struct Map_info *Map, const struct bound_box *Box, struct boxlist *list)
Select areas with bounding boxes by box.
Definition: sindex.c:125
int Vect_select_lines_by_box(struct Map_info *Map, const struct bound_box *Box, int type, struct boxlist *list)
Select lines with bounding boxes by box.
Definition: sindex.c:34
int Vect_point_in_poly(double, double, const struct line_pnts *)
Determines if a point (X,Y) is inside a polygon.
Definition: Vlib/poly.c:826
Bounding box.
Definition: dig_structs.h:65
int Vect_select_areas_by_polygon(struct Map_info *Map, struct line_pnts *Polygon, int nisles, struct line_pnts **Isles, struct ilist *List)
Select areas by Polygon with optional isles.
Definition: sindex.c:324
int * id
Array of ids.
Definition: dig_structs.h:1755
plus_t n_klines
Current number of kernels.
Definition: dig_structs.h:922
double W
West.
Definition: dig_structs.h:82
int Vect_reset_list(struct ilist *)
Reset ilist structure.
Vector geometry.
Definition: dig_structs.h:1574
struct P_line ** Line
Array of vector geometries.
Definition: dig_structs.h:887
int Vect_select_nodes_by_box(struct Map_info *Map, const struct bound_box *Box, struct ilist *list)
Select nodes by box.
Definition: sindex.c:194
float Box[8][3]
Vertices for box.
Definition: gsd_objs.c:1421
plus_t n_plines
Current number of points.
Definition: dig_structs.h:902
int dig_select_isles(struct Plus_head *, const struct bound_box *, struct boxlist *)
Select isles with boxes by box.
Definition: spindex.c:956
struct P_area ** Area
Array of areas.
Definition: dig_structs.h:891
int dig_line_box(const struct line_pnts *, struct bound_box *)
int n_points
Number of points.
Definition: dig_structs.h:1692
#define GV_CENTROID
Definition: dig_defines.h:185
int n_values
Number of values in the list.
Definition: gis.h:709
double E
East.
Definition: dig_structs.h:78
#define NULL
Definition: ccmath.h:32
int n_values
Number of items in the list.
Definition: dig_structs.h:1767
#define GV_POINT
Feature types used in memory on run time (may change)
Definition: dig_defines.h:182
struct bound_box * box
Array of bounding boxes.
Definition: dig_structs.h:1759
int dig_boxlist_add(struct boxlist *, int, const struct bound_box *)
#define GV_LINE
Definition: dig_defines.h:183
double N
North.
Definition: dig_structs.h:70
double * x
Array of X coordinates.
Definition: dig_structs.h:1680
char type
Line type.
Definition: dig_structs.h:1586
Feature geometry info - coordinates.
Definition: dig_structs.h:1675
#define GV_FACE
Definition: dig_defines.h:186
Basic topology-related info.
Definition: dig_structs.h:784
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
int dig_select_nodes(struct Plus_head *, const struct bound_box *, struct ilist *)
Select nodes by bbox.
Definition: spindex.c:663
#define GV_BOUNDARY
Definition: dig_defines.h:184
double B
Bottom.
Definition: dig_structs.h:90
struct Plus_head plus
Plus info (topology, version, ...)
Definition: dig_structs.h:1286
int Vect_select_isles_by_box(struct Map_info *Map, const struct bound_box *Box, struct boxlist *list)
Select isles with bounding boxes by box.
Definition: sindex.c:171
double T
Top.
Definition: dig_structs.h:86
int dig_init_boxlist(struct boxlist *, int)
plus_t n_flines
Current number of faces.
Definition: dig_structs.h:918
int Vect_get_isle_area(const struct Map_info *, int)
Returns area id for isle.
int Vect_reset_boxlist(struct boxlist *)
Reset boxlist structure.
int Vect_find_area(struct Map_info *, double, double)
Find the nearest area.
Vector map info.
Definition: dig_structs.h:1259
#define PORT_DOUBLE_MAX
Limits for portable types.
Definition: dig_defines.h:66
double * y
Array of Y coordinates.
Definition: dig_structs.h:1684
plus_t n_llines
Current number of lines.
Definition: dig_structs.h:906
struct ilist * Vect_new_list(void)
Creates and initializes a struct ilist.
int dig_select_areas(struct Plus_head *, const struct bound_box *, struct boxlist *)
Select areas with boxes by box.
Definition: spindex.c:862
Definition: manage.h:4
int Vect_list_append(struct ilist *, int)
Append new item to the end of list if not yet present.
List of bounding boxes with id.
Definition: dig_structs.h:1750
double S
South.
Definition: dig_structs.h:74
List of integers.
Definition: gis.h:700
plus_t n_clines
Current number of centroids.
Definition: dig_structs.h:914
const char * G_getenv_nofatal(const char *)
Get environment variable.
Definition: env.c:398
int * value
Array of values.
Definition: gis.h:705
struct boxlist * Vect_new_boxlist(int)
Creates and initializes a struct boxlist.
plus_t n_blines
Current number of boundaries.
Definition: dig_structs.h:910
int Vect_read_line(const struct Map_info *, struct line_pnts *, struct line_cats *, int)
Read vector feature (topological level required)
int G_debug(int, const char *,...) __attribute__((format(printf
int Vect_get_line_areas(const struct Map_info *, int, int *, int *)
Get area id on the left and right side of the boundary.
Definition: level_two.c:350
int Vect_line_check_intersection(struct line_pnts *, struct line_pnts *, int)
Check if 2 lines intersect.
#define GV_KERNEL
Definition: dig_defines.h:187
int dig_select_lines(struct Plus_head *, const struct bound_box *, struct boxlist *)
Select lines with boxes by box.
Definition: spindex.c:750