GRASS GIS 7 Programmer's Manual  7.7.svn(2018)-r73574
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 {
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(const struct Map_info *Map,
69  int isle, 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 &&
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 {
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(const 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 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;
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(const 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 /*!
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 {
239  const struct Plus_head *Plus;
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)
248  G_fatal_error(_("Attempt to read topo for dead area (%d)"), area);
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
262  \return 0 area not found
263  */
264 int Vect_get_isle_area(const struct Map_info *Map, int isle)
265 {
266  const struct Plus_head *Plus;
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)
275  G_fatal_error(_("Attempt to read topo for dead isle (%d)"), isle);
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 {
293  const struct Plus_head *Plus;
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;
339  const struct Plus_head *Plus;
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 {
378  const struct Plus_head *Plus;
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) {
432  Vect_read_line(Map, NULL, Cats, centroid);
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
Lines cache for reading feature.
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.
struct Plus_head plus
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