GRASS GIS 7 Programmer's Manual  7.7.svn(2018)-r73574
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(const struct Map_info *Map,
38  int area, 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(const struct Map_info *Map,
69  int isle, 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 &&
88  Map->fInfo.pg.toposchema_name &&
89  Map->fInfo.pg.cache.ctype != CACHE_MAP) {
90 #ifdef HAVE_POSTGRES
91  /* PostGIS Topology */
92  return Vect__get_area_points_pg(Map, Isle->lines, Isle->n_lines, 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(const 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(const 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 structur
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(const 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(const 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 /*!
228  \brief Returns isle id for area
229
230  \param Map pointer to Map_info structure
231  \param area area id
232  \param isle isle index (0 .. nisles - 1)
233
234  \return isle id
235  \return 0 if no isle found
236  */
237 int Vect_get_area_isle(const struct Map_info *Map, int area, int isle)
238 {
240  struct P_area *Area;
241
242  G_debug(3, "Vect_get_area_isle(): area = %d isle = %d", area, isle);
243
244  Plus = &(Map->plus);
245  Area = Plus->Area[area];
246
247  if (Area == NULL)
249
250  G_debug(3, " -> isle = %d", Area->isles[isle]);
251
252  return (Area->isles[isle]);
253 }
254
255 /*!
256  \brief Returns area id for isle
257
258  \param Map vector
259  \param isle isle number (0 .. nisles - 1)
260
261  \return area id
263  */
264 int Vect_get_isle_area(const struct Map_info *Map, int isle)
265 {
267  struct P_isle *Isle;
268
269  G_debug(3, "Vect_get_isle_area(): isle = %d", isle);
270
271  Plus = &(Map->plus);
272  Isle = Plus->Isle[isle];
273
274  if (Isle == NULL)
276
277  G_debug(3, " -> area = %d", Isle->area);
278
279  return (Isle->area);
280 }
281
282 /*!
283  \brief Returns perimeter of area with perimeter of isles
284
285  \param Map pointer to Map_info structure
286  \param area area id
287
288  \return perimeter of area with perimeters of isles in meters
289  */
290
291 double Vect_get_area_perimeter(const struct Map_info *Map, int area)
292 {
294  struct P_area *Area;
295  struct line_pnts *Points;
296  double d;
297  int i;
298
299  G_debug(3, "Vect_get_area_perimeter(): area = %d", area);
300
301  Points = Vect_new_line_struct();
302  Plus = &(Map->plus);
303  Area = Plus->Area[area];
304
305  Vect_get_area_points(Map, area, Points);
306  Vect_line_prune(Points);
307  d = Vect_line_geodesic_length(Points);
308
309  /* adding island perimeters */
310  for (i = 0; i < Area->n_isles; i++) {
311  Vect_get_isle_points(Map, Area->isles[i], Points);
312  Vect_line_prune(Points);
313  d += Vect_line_geodesic_length(Points);
314  }
315
316  Vect_destroy_line_struct(Points);
317
318  G_debug(3, " perimeter = %f", d);
319
320  return (d);
321 }
322
323 /*!
324  \brief Check if point is in area
325
326  \param x,y point coordinates
327  \param Map pointer to Map_info structure
328  \param area area id
329  \param box area bounding box
330
331  \return 0 if point is outside area
332  \return 1 if point is inside area
333  \return 2 if point is on the area's outer ring
334  */
335 int Vect_point_in_area(double x, double y, const struct Map_info *Map,
336  int area, struct bound_box *box)
337 {
338  int i, isle;
340  struct P_area *Area;
341  struct bound_box ibox;
342  int poly;
343
344  Plus = &(Map->plus);
345  Area = Plus->Area[area];
346  if (Area == NULL)
347  return 0;
348
349  poly = Vect_point_in_area_outer_ring(x, y, Map, area, box);
350  if (poly == 0)
351  return 0;
352
353  if (poly == 2) /* includes area boundary, OK? */
354  return 2;
355
356  /* check if in islands */
357  for (i = 0; i < Area->n_isles; i++) {
358  isle = Area->isles[i];
359  Vect_get_isle_box(Map, isle, &ibox);
360  poly = Vect_point_in_island(x, y, Map, isle, &ibox);
361  if (poly >= 1)
362  return 0; /* excludes island boundary (poly == 2), OK? */
363  }
364
365  return 1;
366 }
367
368 /*!
369  \brief Returns area of area without areas of isles
370
371  \param Map pointer to Map_info structure
372  \param area area id
373
374  \return area of area without areas of isles
375  */
376 double Vect_get_area_area(const struct Map_info *Map, int area)
377 {
379  struct P_area *Area;
380  struct line_pnts *Points;
381  double size;
382  int i;
383  static int first_time = 1;
384
385  G_debug(3, "Vect_get_area_area(): area = %d", area);
386
387  if (first_time == 1) {
389  first_time = 0;
390  }
391
392  Points = Vect_new_line_struct();
393  Plus = &(Map->plus);
394  Area = Plus->Area[area];
395
396  Vect_get_area_points(Map, area, Points);
397  Vect_line_prune(Points);
398  size = G_area_of_polygon(Points->x, Points->y, Points->n_points);
399
400  /* subtracting island areas */
401  for (i = 0; i < Area->n_isles; i++) {
402  Vect_get_isle_points(Map, Area->isles[i], Points);
403  Vect_line_prune(Points);
404  size -= G_area_of_polygon(Points->x, Points->y, Points->n_points);
405  }
406
407  Vect_destroy_line_struct(Points);
408
409  G_debug(3, " area = %f", size);
410
411  return (size);
412 }
413
414 /*!
415  \brief Get area categories
416
417  \param Map pointer to Map_info structure
418  \param area area id
419  \param[out] Cats list of categories
420
421  \return 0 centroid found (but may be without categories)
422  \return 1 no centroid found
423  */
424 int Vect_get_area_cats(const struct Map_info *Map, int area, struct line_cats *Cats)
425 {
426  int centroid;
427
428  Vect_reset_cats(Cats);
429
430  centroid = Vect_get_area_centroid(Map, area);
431  if (centroid > 0) {
433  }
434  else {
435  return 1; /* no centroid */
436  }
437
438
439  return 0;
440 }
441
442 /*!
443  \brief Find FIRST category of given field and area
444
445  \param Map pointer to Map_info structure
446  \param area area id
447  \param field layer number
448
449  \return first found category of given field
450  \return -1 no centroid or no category found
451  */
452 int Vect_get_area_cat(const struct Map_info *Map, int area, int field)
453 {
454  int i;
455  static struct line_cats *Cats = NULL;
456
457  if (!Cats)
458  Cats = Vect_new_cats_struct();
459  else
460  Vect_reset_cats(Cats);
461
462  if (Vect_get_area_cats(Map, area, Cats) == 1 || Cats->n_cats == 0) {
463  return -1;
464  }
465
466  for (i = 0; i < Cats->n_cats; i++) {
467  if (Cats->field[i] == field) {
468  return Cats->cat[i];
469  }
470  }
471
472  return -1;
473 }
474
475 /*!
476  \brief Get area boundary points (internal use only)
477
478  For PostGIS Topology calls Vect__get_area_points_pg() otherwise
479  Vect__get_area_points_nat(),
480
481  \param Map pointer to Map_info struct
482  \param lines array of boundary lines
483  \param n_lines number of lines in array
484  \param[out] APoints pointer to output line_pnts struct
485
486  \return number of points
487  \return -1 on error
488 */
489 int Vect__get_area_points(const struct Map_info *Map, const plus_t *lines, int n_lines,
490  struct line_pnts *BPoints)
491 {
492  if (Map->format == GV_FORMAT_POSTGIS &&
493  Map->fInfo.pg.toposchema_name &&
494  Map->fInfo.pg.cache.ctype != CACHE_MAP) {
495 #ifdef HAVE_POSTGRES
496  /* PostGIS Topology */
497  return Vect__get_area_points_pg(Map, lines, n_lines, BPoints);
498 #else
499  G_fatal_error(_("GRASS is not compiled with PostgreSQL support"));
500 #endif
501  }
502  /* native format */
503  return Vect__get_area_points_nat(Map, lines, n_lines, BPoints);
504 }
505
506
507 /*!
508  \brief Get area boundary points (native format)
509
510  Used by Vect_build_line_area() and Vect_get_area_points().
511
512  \param Map pointer to Map_info struct
513  \param lines array of boundary lines
514  \param n_lines number of lines in array
515  \param[out] APoints pointer to output line_pnts struct
516
517  \return number of points
518  \return -1 on error
519 */
520 int Vect__get_area_points_nat(const struct Map_info *Map, const plus_t *lines, int n_lines,
521  struct line_pnts *BPoints)
522 {
523  int i, line, aline, dir;
524  static struct line_pnts *Points;
525
526  if (!Points)
527  Points = Vect_new_line_struct();
528
529  Vect_reset_line(BPoints);
530  for (i = 0; i < n_lines; i++) {
531  line = lines[i];
532  aline = abs(line);
533  G_debug(5, " append line(%d) = %d", i, line);
534
535  if (0 > Vect_read_line(Map, Points, NULL, aline))
536  return -1;
537
538  dir = line > 0 ? GV_FORWARD : GV_BACKWARD;
539  Vect_append_points(BPoints, Points, dir);
540  BPoints->n_points--; /* skip last point, avoids duplicates */
541  }
542  BPoints->n_points++; /* close polygon */
543
544  return BPoints->n_points;
545 }
int Vect_get_area_cat(const struct Map_info *Map, int area, int field)
Find FIRST category of given field and area.
int Vect__get_area_points_nat(const struct Map_info *Map, const plus_t *lines, int n_lines, struct line_pnts *BPoints)
Get area boundary points (native format)
char * toposchema_name
Topology schema name and id.
Definition: dig_structs.h:699
int Vect_get_isle_box(const struct Map_info *Map, int isle, struct bound_box *Box)
Get bounding box of isle.
Bounding box.
Definition: dig_structs.h:65
#define GV_FORWARD
Line direction indicator forward/backward.
Definition: dig_defines.h:178
int plus_t
plus_t size
Definition: dig_structs.h:41
int Vect_get_area_cats(const struct Map_info *Map, int area, struct line_cats *Cats)
Get area categories.
int Vect_point_in_area(double x, double y, const struct Map_info *Map, int area, struct bound_box *box)
Check if point is in area.
int Vect_get_area_boundaries(const struct Map_info *Map, int area, struct ilist *List)
Creates list of boundaries for given area.
double G_area_of_polygon(const double *x, const double *y, int n)
Area in square meters of polygon.
Definition: gis/area.c:159
#define GV_BACKWARD
Definition: dig_defines.h:179
Isle (topology) info.
Definition: dig_structs.h:1646
struct line_pnts * Vect_new_line_struct()
Creates and initializes a line_pnts structure.
Definition: line.c:45
struct P_area ** Area
Array of areas.
Definition: dig_structs.h:891
int Vect__get_area_points(const struct Map_info *Map, const plus_t *lines, int n_lines, struct line_pnts *BPoints)
Get area boundary points (internal use only)
int n_points
Number of points.
Definition: dig_structs.h:1692
int n_values
Number of values in the list.
Definition: gis.h:676
void Vect_destroy_line_struct(struct line_pnts *p)
Frees all memory associated with a line_pnts structure, including the structure itself.
Definition: line.c:77
struct Format_info fInfo
Format info for non-native formats.
Definition: dig_structs.h:1415
struct P_isle ** Isle
Array of isles.
Definition: dig_structs.h:895
plus_t * isles
1st generation interior islands
Definition: dig_structs.h:1640
int Vect_reset_cats(struct line_cats *Cats)
Reset category structure to make sure cats structure is clean to be re-used.
plus_t centroid
Number of first centroid within area.
Definition: dig_structs.h:1628
#define NULL
Definition: ccmath.h:32
plus_t n_lines
Number of boundary lines.
Definition: dig_structs.h:1651
#define x
int Vect_get_area_centroid(const struct Map_info *Map, int area)
Returns centroid id for given area.
int Vect_point_in_area_outer_ring(double X, double Y, const struct Map_info *Map, int area, struct bound_box *box)
Determines if a point (X,Y) is inside an area outer ring. Islands are not considered.
Definition: Vlib/poly.c:856
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:160
Feature category info.
Definition: dig_structs.h:1702
int Vect_append_points(struct line_pnts *Points, const struct line_pnts *APoints, int direction)
Appends points to the end of a line.
Definition: line.c:337
double * x
Array of X coordinates.
Definition: dig_structs.h:1680
Feature geometry info - coordinates.
Definition: dig_structs.h:1675
int Vect_point_in_island(double X, double Y, const struct Map_info *Map, int isle, struct bound_box *box)
Determines if a point (X,Y) is inside an island.
Definition: Vlib/poly.c:925
plus_t * lines
List of boundary lines.
Definition: dig_structs.h:1621
double Vect_get_area_area(const struct Map_info *Map, int area)
Returns area of area without areas of isles.
int ctype
Cache type.
Definition: dig_structs.h:507
Basic topology-related info.
Definition: dig_structs.h:784
int Vect_get_isle_boundaries(const struct Map_info *Map, int isle, struct ilist *List)
Creates list of boundaries for given isle.
struct Format_info_cache cache
Definition: dig_structs.h:683
int G_begin_polygon_area_calculations(void)
Begin polygon area calculations.
Definition: gis/area.c:120
struct Format_info_pg pg
PostGIS info.
Definition: dig_structs.h:726
int n_cats
Number of categories attached to element.
Definition: dig_structs.h:1715
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
int * cat
Array of categories.
Definition: dig_structs.h:1711
int Vect_get_isle_area(const struct Map_info *Map, int isle)
Returns area id for isle.
plus_t * lines
List of boundary lines.
Definition: dig_structs.h:1662
int Vect_reset_list(struct ilist *list)
Reset ilist structure.
Plus info (topology, version, ...)
Definition: dig_structs.h:1286
int Vect_list_append(struct ilist *list, int val)
Append new item to the end of list if not yet present.
plus_t n_isles
Number of islands inside.
Definition: dig_structs.h:1632
int Vect_get_area_num_isles(const struct Map_info *Map, int area)
Returns number of isles for given area.
Vector map info.
Definition: dig_structs.h:1259
#define GV_FORMAT_POSTGIS
PostGIS format.
Definition: dig_defines.h:89
double * y
Array of Y coordinates.
Definition: dig_structs.h:1684
Area (topology) info.
Definition: dig_structs.h:1605
struct line_cats * Vect_new_cats_struct()
Creates and initializes line_cats structure.
int Vect_line_prune(struct line_pnts *Points)
Remove duplicate points, i.e. zero length segments.
Definition: line.c:281
double Vect_get_area_perimeter(const struct Map_info *Map, int area)
Returns perimeter of area with perimeter of isles.
#define _(str)
Definition: glocale.h:13
List of integers.
Definition: gis.h:667
int Vect_read_line(const struct Map_info *Map, struct line_pnts *line_p, struct line_cats *line_c, int line)
Read vector feature (topological level required)
int Vect_get_area_points(const struct Map_info *Map, int area, struct line_pnts *BPoints)
Returns polygon array of points (outer ring) of given area.
int format
Map format (native, ogr, postgis)
Definition: dig_structs.h:1271
void Vect_reset_line(struct line_pnts *Points)
Reset line.
Definition: line.c:130
int * field
Array of layers (fields)
Definition: dig_structs.h:1707
double Vect_line_geodesic_length(const struct line_pnts *Points)
Calculate line length.
Definition: line.c:604
int Vect_get_isle_points(const struct Map_info *Map, int isle, struct line_pnts *BPoints)
Returns polygon array of points for given isle.
plus_t n_lines
Number of boundary lines.
Definition: dig_structs.h:1610
plus_t area
Area it exists w/in, if any.
Definition: dig_structs.h:1669
int Vect_get_area_isle(const struct Map_info *Map, int area, int isle)
Returns isle id for area.
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204