GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-55b2a2bccc
vector/Vlib/area.c
Go to the documentation of this file.
1 /*!
2  \file lib/vector/Vlib/area.c
3 
4  \brief Vector library - area-related functions
5 
6  Higher level functions for reading/writing/manipulating vectors.
7 
8  (C) 2001-2009, 2011-2013 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 Original author CERL, probably Dave Gerdes or Mike Higgins.
14  \author Update to GRASS 5.7 Radim Blazek and David D. Gray.
15  */
16 
17 #include <stdlib.h>
18 #include <grass/vector.h>
19 #include <grass/glocale.h>
20 
21 #ifdef HAVE_POSTGRES
22 #include "pg_local_proto.h"
23 #endif
24 
25 #include "local_proto.h"
26 
27 /*!
28  \brief Returns polygon array of points (outer ring) of given area
29 
30  \param Map pointer to Map_info structure
31  \param area area id
32  \param[out] BPoints points array
33 
34  \return number of points
35  \return -1 on error
36  */
37 int Vect_get_area_points(struct Map_info *Map, int area,
38  struct line_pnts *BPoints)
39 {
40  const struct Plus_head *Plus;
41  struct P_area *Area;
42 
43  G_debug(3, "Vect_get_area_points(): area = %d", area);
44  Vect_reset_line(BPoints);
45 
46  Plus = &(Map->plus);
47  Area = Plus->Area[area];
48 
49  if (Area == NULL) { /* dead area */
50  G_warning(_("Attempt to read points of nonexistent area"));
51  return -1; /* error, because we should not read dead areas */
52  }
53 
54  G_debug(3, " n_lines = %d", Area->n_lines);
55  return Vect__get_area_points(Map, Area->lines, Area->n_lines, BPoints);
56 }
57 
58 /*!
59  \brief Returns polygon array of points for given isle
60 
61  \param Map pointer to Map_info structure
62  \param isle island id
63  \param[out] BPoints points array
64 
65  \return number of points
66  \return -1 on error
67  */
68 int Vect_get_isle_points(struct Map_info *Map, int isle,
69  struct line_pnts *BPoints)
70 {
71  const struct Plus_head *Plus;
72  struct P_isle *Isle;
73 
74  G_debug(3, "Vect_get_isle_points(): isle = %d", isle);
75  Vect_reset_line(BPoints);
76 
77  Plus = &(Map->plus);
78  Isle = Plus->Isle[isle];
79 
80  if (Isle == NULL) { /* dead isle */
81  G_warning(_("Attempt to read points of nonexistent isle"));
82  return -1; /* error, because we should not read dead isles */
83  }
84 
85  G_debug(3, " n_lines = %d", Isle->n_lines);
86 
87  if (Map->format == GV_FORMAT_POSTGIS && Map->fInfo.pg.toposchema_name &&
88  Map->fInfo.pg.cache.ctype != CACHE_MAP) {
89 #ifdef HAVE_POSTGRES
90  /* PostGIS Topology */
91  return Vect__get_area_points_pg(Map, Isle->lines, Isle->n_lines,
92  BPoints);
93 #else
94  G_fatal_error(_("GRASS is not compiled with PostgreSQL support"));
95 #endif
96  }
97  /* native format */
98  return Vect__get_area_points_nat(Map, Isle->lines, Isle->n_lines, BPoints);
99 }
100 
101 /*!
102  \brief Returns centroid id for given area
103 
104  \param Map pointer to Map_info structure
105  \param area area id
106 
107  \return centroid id of area
108  \return 0 if no centroid found
109  */
110 int Vect_get_area_centroid(struct Map_info *Map, int area)
111 {
112  const struct Plus_head *Plus;
113  struct P_area *Area;
114 
115  G_debug(3, "Vect_get_area_centroid(): area = %d", area);
116 
117  Plus = &(Map->plus);
118  Area = Plus->Area[area];
119 
120  if (Area == NULL)
121  G_fatal_error(_("Attempt to read topo for dead area (%d)"), area);
122 
123  return (Area->centroid);
124 }
125 
126 /*!
127  \brief Creates list of boundaries for given area
128 
129  Note that ids in <b>List</b> can be negative. The sign indicates in
130  which direction the boundary should be read (negative for
131  backward).
132 
133  \param Map pointer to Map_info structure
134  \param area area id
135  \param[out] List pointer to list of boundaries
136 
137  \return number of boundaries
138  */
139 int Vect_get_area_boundaries(struct Map_info *Map, int area, struct ilist *List)
140 {
141  int i, line;
142  const struct Plus_head *Plus;
143  struct P_area *Area;
144 
145  G_debug(3, "Vect_get_area_boundaries(): area = %d", area);
146 
147  Vect_reset_list(List);
148 
149  Plus = &(Map->plus);
150  Area = Plus->Area[area];
151 
152  if (Area == NULL)
153  G_fatal_error(_("Attempt to read topo for dead area (%d)"), area);
154 
155  for (i = 0; i < Area->n_lines; i++) {
156  line = Area->lines[i];
157  Vect_list_append(List, line);
158  }
159 
160  return (List->n_values);
161 }
162 
163 /*!
164  \brief Creates list of boundaries for given isle
165 
166  Note that ids in <b>List</b> can be negative. The sign indicates in
167  which direction the boundary should be read (negative for forward).
168 
169  \param Map pointer to Map_info structure
170  \param isle island number
171  \param[out] List pointer to list where boundaries are stored
172 
173  \return number of boundaries
174  */
175 int Vect_get_isle_boundaries(struct Map_info *Map, int isle, struct ilist *List)
176 {
177  int i, line;
178  const struct Plus_head *Plus;
179  struct P_isle *Isle;
180 
181  G_debug(3, "Vect_get_isle_boundaries(): isle = %d", isle);
182 
183  Vect_reset_list(List);
184 
185  Plus = &(Map->plus);
186  Isle = Plus->Isle[isle];
187 
188  if (Isle == NULL)
189  G_fatal_error(_("Attempt to read topo for dead isle (%d)"), isle);
190 
191  for (i = 0; i < Isle->n_lines; i++) {
192  line = Isle->lines[i];
193  Vect_list_append(List, line);
194  }
195 
196  return (List->n_values);
197 }
198 
199 /*!
200  \brief Returns number of isles for given area
201 
202  \param Map pointer to Map_info structure
203  \param area area id
204 
205  \return number of isles for area
206  \return 0 if area not found
207  */
208 int Vect_get_area_num_isles(struct Map_info *Map, int area)
209 {
210  const struct Plus_head *Plus;
211  struct P_area *Area;
212 
213  G_debug(3, "Vect_get_area_num_isles(): area = %d", area);
214 
215  Plus = &(Map->plus);
216  Area = Plus->Area[area];
217 
218  if (Area == NULL)
219  G_fatal_error(_("Attempt to read topo for dead area (%d)"), area);
220 
221  G_debug(3, " n_isles = %d", Area->n_isles);
222 
223  return (Area->n_isles);
224 }
225 
226 /*!
227  \brief Returns isle id for area
228 
229  \param Map pointer to Map_info structure
230  \param area area id
231  \param isle isle index (0 .. nisles - 1)
232 
233  \return isle id
234  \return 0 if no isle found
235  */
236 int Vect_get_area_isle(struct Map_info *Map, int area, int isle)
237 {
238  const struct Plus_head *Plus;
239  struct P_area *Area;
240 
241  G_debug(3, "Vect_get_area_isle(): area = %d isle = %d", area, isle);
242 
243  Plus = &(Map->plus);
244  Area = Plus->Area[area];
245 
246  if (Area == NULL)
247  G_fatal_error(_("Attempt to read topo for dead area (%d)"), area);
248 
249  G_debug(3, " -> isle = %d", Area->isles[isle]);
250 
251  return (Area->isles[isle]);
252 }
253 
254 /*!
255  \brief Returns area id for isle
256 
257  \param Map vector
258  \param isle isle number (0 .. nisles - 1)
259 
260  \return area id
261  \return 0 area not found
262  */
263 int Vect_get_isle_area(struct Map_info *Map, int isle)
264 {
265  const struct Plus_head *Plus;
266  struct P_isle *Isle;
267 
268  G_debug(3, "Vect_get_isle_area(): isle = %d", isle);
269 
270  Plus = &(Map->plus);
271  Isle = Plus->Isle[isle];
272 
273  if (Isle == NULL)
274  G_fatal_error(_("Attempt to read topo for dead isle (%d)"), isle);
275 
276  G_debug(3, " -> area = %d", Isle->area);
277 
278  return (Isle->area);
279 }
280 
281 /*!
282  \brief Returns perimeter of area with perimeter of isles
283 
284  \param Map pointer to Map_info structure
285  \param area area id
286 
287  \return perimeter of area with perimeters of isles in meters
288  */
289 
290 double Vect_get_area_perimeter(struct Map_info *Map, int area)
291 {
292  const struct Plus_head *Plus;
293  struct P_area *Area;
294  struct line_pnts *Points;
295  double d;
296  int i;
297 
298  G_debug(3, "Vect_get_area_perimeter(): area = %d", area);
299 
300  Points = Vect_new_line_struct();
301  Plus = &(Map->plus);
302  Area = Plus->Area[area];
303 
304  Vect_get_area_points(Map, area, Points);
305  Vect_line_prune(Points);
306  d = Vect_line_geodesic_length(Points);
307 
308  /* adding island perimeters */
309  for (i = 0; i < Area->n_isles; i++) {
310  Vect_get_isle_points(Map, Area->isles[i], Points);
311  Vect_line_prune(Points);
312  d += Vect_line_geodesic_length(Points);
313  }
314 
315  Vect_destroy_line_struct(Points);
316 
317  G_debug(3, " perimeter = %f", d);
318 
319  return (d);
320 }
321 
322 /*!
323  \brief Check if point is in area
324 
325  \param x,y point coordinates
326  \param Map pointer to Map_info structure
327  \param area area id
328  \param box area bounding box
329 
330  \return 0 if point is outside area
331  \return 1 if point is inside area
332  \return 2 if point is on the area's outer ring
333  */
334 int Vect_point_in_area(double x, double y, struct Map_info *Map, int area,
335  struct bound_box *box)
336 {
337  int i, isle;
338  const struct Plus_head *Plus;
339  struct P_area *Area;
340  struct bound_box ibox;
341  int poly;
342 
343  Plus = &(Map->plus);
344  Area = Plus->Area[area];
345  if (Area == NULL)
346  return 0;
347 
348  poly = Vect_point_in_area_outer_ring(x, y, Map, area, box);
349  if (poly == 0)
350  return 0;
351 
352  if (poly == 2) /* includes area boundary, OK? */
353  return 2;
354 
355  /* check if in islands */
356  for (i = 0; i < Area->n_isles; i++) {
357  isle = Area->isles[i];
358  Vect_get_isle_box(Map, isle, &ibox);
359  poly = Vect_point_in_island(x, y, Map, isle, &ibox);
360  if (poly >= 1)
361  return 0; /* excludes island boundary (poly == 2), OK? */
362  }
363 
364  return 1;
365 }
366 
367 /*!
368  \brief Returns area of area without areas of isles
369 
370  \param Map pointer to Map_info structure
371  \param area area id
372 
373  \return area of area without areas of isles
374  */
375 double Vect_get_area_area(struct Map_info *Map, int area)
376 {
377  const struct Plus_head *Plus;
378  struct P_area *Area;
379  struct line_pnts *Points;
380  double size;
381  int i;
382  static int first_time = 1;
383 
384  G_debug(3, "Vect_get_area_area(): area = %d", area);
385 
386  if (first_time == 1) {
388  first_time = 0;
389  }
390 
391  Points = Vect_new_line_struct();
392  Plus = &(Map->plus);
393  Area = Plus->Area[area];
394 
395  Vect_get_area_points(Map, area, Points);
396  Vect_line_prune(Points);
397  size = G_area_of_polygon(Points->x, Points->y, Points->n_points);
398 
399  /* subtracting island areas */
400  for (i = 0; i < Area->n_isles; i++) {
401  Vect_get_isle_points(Map, Area->isles[i], Points);
402  Vect_line_prune(Points);
403  size -= G_area_of_polygon(Points->x, Points->y, Points->n_points);
404  }
405 
406  Vect_destroy_line_struct(Points);
407 
408  G_debug(3, " area = %f", size);
409 
410  return (size);
411 }
412 
413 /*!
414  \brief Get area categories
415 
416  \param Map pointer to Map_info structure
417  \param area area id
418  \param[out] Cats list of categories
419 
420  \return 0 centroid found (but may be without categories)
421  \return 1 no centroid found
422  */
423 int Vect_get_area_cats(struct Map_info *Map, int area, struct line_cats *Cats)
424 {
425  int centroid;
426 
427  Vect_reset_cats(Cats);
428 
429  centroid = Vect_get_area_centroid(Map, area);
430  if (centroid > 0) {
431  Vect_read_line(Map, NULL, Cats, centroid);
432  }
433  else {
434  return 1; /* no centroid */
435  }
436 
437  return 0;
438 }
439 
440 /*!
441  \brief Find FIRST category of given field and area
442 
443  \param Map pointer to Map_info structure
444  \param area area id
445  \param field layer number
446 
447  \return first found category of given field
448  \return -1 no centroid or no category found
449  */
450 int Vect_get_area_cat(struct Map_info *Map, int area, int field)
451 {
452  int i;
453  static struct line_cats *Cats = NULL;
454 
455  if (!Cats)
456  Cats = Vect_new_cats_struct();
457  else
458  Vect_reset_cats(Cats);
459 
460  if (Vect_get_area_cats(Map, area, Cats) == 1 || Cats->n_cats == 0) {
461  return -1;
462  }
463 
464  for (i = 0; i < Cats->n_cats; i++) {
465  if (Cats->field[i] == field) {
466  return Cats->cat[i];
467  }
468  }
469 
470  return -1;
471 }
472 
473 /*!
474  \brief Get area boundary points (internal use only)
475 
476  For PostGIS Topology calls Vect__get_area_points_pg() otherwise
477  Vect__get_area_points_nat(),
478 
479  \param Map pointer to Map_info struct
480  \param lines array of boundary lines
481  \param n_lines number of lines in array
482  \param[out] APoints pointer to output line_pnts struct
483 
484  \return number of points
485  \return -1 on error
486  */
487 int Vect__get_area_points(struct Map_info *Map, const plus_t *lines,
488  int n_lines, struct line_pnts *BPoints)
489 {
490  if (Map->format == GV_FORMAT_POSTGIS && Map->fInfo.pg.toposchema_name &&
491  Map->fInfo.pg.cache.ctype != CACHE_MAP) {
492 #ifdef HAVE_POSTGRES
493  /* PostGIS Topology */
494  return Vect__get_area_points_pg(Map, lines, n_lines, BPoints);
495 #else
496  G_fatal_error(_("GRASS is not compiled with PostgreSQL support"));
497 #endif
498  }
499  /* native format */
500  return Vect__get_area_points_nat(Map, lines, n_lines, BPoints);
501 }
502 
503 /*!
504  \brief Get area boundary points (native format)
505 
506  Used by Vect_build_line_area() and Vect_get_area_points().
507 
508  \param Map pointer to Map_info struct
509  \param lines array of boundary lines
510  \param n_lines number of lines in array
511  \param[out] APoints pointer to output line_pnts struct
512 
513  \return number of points
514  \return -1 on error
515  */
516 int Vect__get_area_points_nat(struct Map_info *Map, const plus_t *lines,
517  int n_lines, struct line_pnts *BPoints)
518 {
519  int i, line, aline, dir;
520  static struct line_pnts *Points;
521 
522  if (!Points)
523  Points = Vect_new_line_struct();
524 
525  Vect_reset_line(BPoints);
526  for (i = 0; i < n_lines; i++) {
527  line = lines[i];
528  aline = abs(line);
529  G_debug(5, " append line(%d) = %d", i, line);
530 
531  if (0 > Vect_read_line(Map, Points, NULL, aline))
532  return -1;
533 
534  dir = line > 0 ? GV_FORWARD : GV_BACKWARD;
535  Vect_append_points(BPoints, Points, dir);
536  BPoints->n_points--; /* skip last point, avoids duplicates */
537  }
538  BPoints->n_points++; /* close polygon */
539 
540  return BPoints->n_points;
541 }
int Vect__get_area_points_pg(struct Map_info *Map, const plus_t *lines, int n_lines, struct line_pnts *APoints)
Get area boundary points (PostGIS Topology)
Definition: area_pg.c:37
#define NULL
Definition: ccmath.h:32
double G_area_of_polygon(const double *, const double *, int)
Area in square meters of polygon.
Definition: gis/area.c:159
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void G_warning(const char *,...) __attribute__((format(printf
int G_debug(int, const char *,...) __attribute__((format(printf
int G_begin_polygon_area_calculations(void)
Begin polygon area calculations.
Definition: gis/area.c:120
void Vect_destroy_line_struct(struct line_pnts *)
Frees all memory associated with a line_pnts structure, including the structure itself.
Definition: line.c:77
int Vect_reset_cats(struct line_cats *)
Reset category structure to make sure cats structure is clean to be re-used.
int Vect_point_in_island(double, double, struct Map_info *, int, struct bound_box *)
Determines if a point (X,Y) is inside an island.
Definition: Vlib/poly.c:923
double Vect_line_geodesic_length(const struct line_pnts *)
Calculate line length.
Definition: line.c:602
int Vect_list_append(struct ilist *, int)
Append new item to the end of list if not yet present.
int Vect_read_line(struct Map_info *, struct line_pnts *, struct line_cats *, int)
Read vector feature (topological level required)
int Vect_point_in_area_outer_ring(double, double, struct Map_info *, int, struct bound_box *)
Determines if a point (X,Y) is inside an area outer ring. Islands are not considered.
Definition: Vlib/poly.c:854
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
int Vect_get_isle_box(struct Map_info *, int, struct bound_box *)
Get bounding box of isle.
void Vect_reset_line(struct line_pnts *)
Reset line.
Definition: line.c:129
int Vect_line_prune(struct line_pnts *)
Remove duplicate points, i.e. zero length segments.
Definition: line.c:279
struct line_cats * Vect_new_cats_struct(void)
Creates and initializes line_cats structure.
int Vect_reset_list(struct ilist *)
Reset ilist structure.
int Vect_append_points(struct line_pnts *, const struct line_pnts *, int)
Appends points to the end of a line.
Definition: line.c:335
#define GV_FORMAT_POSTGIS
PostGIS format.
Definition: dig_defines.h:89
#define GV_FORWARD
Line direction indicator forward/backward.
Definition: dig_defines.h:179
#define GV_BACKWARD
Definition: dig_defines.h:180
int plus_t
plus_t size
Definition: dig_structs.h:41
#define _(str)
Definition: glocale.h:10
int ctype
Cache type.
Definition: dig_structs.h:496
char * toposchema_name
Topology schema name and id.
Definition: dig_structs.h:686
struct Format_info_cache cache
Lines cache for reading feature.
Definition: dig_structs.h:670
struct Format_info_pg pg
PostGIS info.
Definition: dig_structs.h:712
Vector map info.
Definition: dig_structs.h:1243
int format
Map format (native, ogr, postgis)
Definition: dig_structs.h:1255
struct Format_info fInfo
Format info for non-native formats.
Definition: dig_structs.h:1400
struct Plus_head plus
Plus info (topology, version, ...)
Definition: dig_structs.h:1270
Area (topology) info.
Definition: dig_structs.h:1583
plus_t n_isles
Number of islands inside.
Definition: dig_structs.h:1609
plus_t * isles
1st generation interior islands
Definition: dig_structs.h:1617
plus_t n_lines
Number of boundary lines.
Definition: dig_structs.h:1587
plus_t * lines
List of boundary lines.
Definition: dig_structs.h:1598
plus_t centroid
Number of first centroid within area.
Definition: dig_structs.h:1605
Isle (topology) info.
Definition: dig_structs.h:1623
plus_t * lines
List of boundary lines.
Definition: dig_structs.h:1638
plus_t n_lines
Number of boundary lines.
Definition: dig_structs.h:1627
plus_t area
Area it exists w/in, if any.
Definition: dig_structs.h:1645
Basic topology-related info.
Definition: dig_structs.h:769
struct P_area ** Area
Array of areas.
Definition: dig_structs.h:875
struct P_isle ** Isle
Array of isles.
Definition: dig_structs.h:879
Bounding box.
Definition: dig_structs.h:64
List of integers.
Definition: gis.h:708
int n_values
Number of values in the list.
Definition: gis.h:716
Feature category info.
Definition: dig_structs.h:1677
int * field
Array of layers (fields)
Definition: dig_structs.h:1681
int * cat
Array of categories.
Definition: dig_structs.h:1685
int n_cats
Number of categories attached to element.
Definition: dig_structs.h:1689
Feature geometry info - coordinates.
Definition: dig_structs.h:1651
double * y
Array of Y coordinates.
Definition: dig_structs.h:1659
double * x
Array of X coordinates.
Definition: dig_structs.h:1655
int n_points
Number of points.
Definition: dig_structs.h:1667
double Vect_get_area_perimeter(struct Map_info *Map, int area)
Returns perimeter of area with perimeter of isles.
int Vect_get_area_isle(struct Map_info *Map, int area, int isle)
Returns isle id for area.
int Vect_point_in_area(double x, double y, struct Map_info *Map, int area, struct bound_box *box)
Check if point is in area.
int Vect_get_area_points(struct Map_info *Map, int area, struct line_pnts *BPoints)
Returns polygon array of points (outer ring) of given area.
int Vect_get_area_cat(struct Map_info *Map, int area, int field)
Find FIRST category of given field and area.
int Vect_get_area_boundaries(struct Map_info *Map, int area, struct ilist *List)
Creates list of boundaries for given area.
int Vect_get_area_num_isles(struct Map_info *Map, int area)
Returns number of isles for given area.
int Vect__get_area_points(struct Map_info *Map, const plus_t *lines, int n_lines, struct line_pnts *BPoints)
Get area boundary points (internal use only)
int Vect_get_area_centroid(struct Map_info *Map, int area)
Returns centroid id for given area.
int Vect_get_isle_points(struct Map_info *Map, int isle, struct line_pnts *BPoints)
Returns polygon array of points for given isle.
int Vect_get_area_cats(struct Map_info *Map, int area, struct line_cats *Cats)
Get area categories.
int Vect_get_isle_boundaries(struct Map_info *Map, int isle, struct ilist *List)
Creates list of boundaries for given isle.
int Vect__get_area_points_nat(struct Map_info *Map, const plus_t *lines, int n_lines, struct line_pnts *BPoints)
Get area boundary points (native format)
double Vect_get_area_area(struct Map_info *Map, int area)
Returns area of area without areas of isles.
int Vect_get_isle_area(struct Map_info *Map, int isle)
Returns area id for isle.
#define x