GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-d2f347782a
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 {
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 {
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 {
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)
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;
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)
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;
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)
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
207  */
208 int Vect_get_area_num_isles(struct Map_info *Map, int area)
209 {
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)
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 {
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)
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
262  */
263 int Vect_get_isle_area(struct Map_info *Map, int isle)
264 {
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)
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 {
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;
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 {
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) {
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
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
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