GRASS GIS 7 Programmer's Manual  7.5.svn(2018)-r72645
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gis/area.c
Go to the documentation of this file.
1 /*!
2  * \file lib/gis/area.c
3  *
4  * \brief GIS Library - Area calculation functions.
5  *
6  * (C) 2001-2009 by the GRASS Development Team
7  *
8  * This program is free software under the GNU General Public License
9  * (>=v2). Read the file COPYING that comes with GRASS for details.
10  *
11  * \author Original author CERL
12  */
13 
14 #include <grass/gis.h>
15 
16 static struct state {
17  struct Cell_head window;
18  double square_meters;
19  int projection;
20 
21  double units_to_meters_squared;
22 
23  /* these next are for lat-long only */
24  int next_row;
25  double north_value;
26  double north;
27  double (*darea0) (double);
28 } state;
29 
30 static struct state *st = &state;
31 
32 /*!
33  * \brief Begin cell area calculations.
34  *
35  * This routine must be called once before any call to
36  * G_area_of_cell_at_row(). It perform all inititalizations needed to
37  * do area calculations for grid cells, based on the current window
38  * "projection" field. It can be used in either planimetric
39  * projections or the latitude-longitude projection.
40  *
41  * \return 0 if the projection is not measurable (ie. imagery or xy)
42  * \return 1 if the projection is planimetric (ie. UTM or SP)
43  * \return 2 if the projection is non-planimetric (ie. latitude-longitude)
44  */
45 
47 {
48  double a, e2;
49  double factor;
50 
51  G_get_set_window(&st->window);
52  switch (st->projection = st->window.proj) {
53  case PROJECTION_LL:
55  if (e2) {
56  G_begin_zone_area_on_ellipsoid(a, e2, st->window.ew_res / 360.0);
57  st->darea0 = G_darea0_on_ellipsoid;
58  }
59  else {
60  G_begin_zone_area_on_sphere(a, st->window.ew_res / 360.0);
61  st->darea0 = G_darea0_on_sphere;
62  }
63  st->next_row = 0;
64  st->north = st->window.north;
65  st->north_value = st->darea0(st->north);
66  return 2;
67  default:
68  st->square_meters = st->window.ns_res * st->window.ew_res;
70  if (factor > 0.0)
71  st->square_meters *= (factor * factor);
72  return (factor > 0.0);
73  }
74 }
75 
76 /*!
77  * \brief Cell area in specified row.
78  *
79  * This routine returns the area in square meters of a cell in the
80  * specified <i>row</i>. This value is constant for planimetric grids
81  * and varies with the row if the projection is latitude-longitude.
82  *
83  * \param row row number
84  *
85  * \return cell area
86  */
87 double G_area_of_cell_at_row(int row)
88 {
89  double south_value;
90  double cell_area;
91 
92  if (st->projection != PROJECTION_LL)
93  return st->square_meters;
94 
95  if (row != st->next_row) {
96  st->north = st->window.north - row * st->window.ns_res;
97  st->north_value = st->darea0(st->north);
98  }
99 
100  st->north -= st->window.ns_res;
101  south_value = st->darea0(st->north);
102  cell_area = st->north_value - south_value;
103 
104  st->next_row = row + 1;
105  st->north_value = south_value;
106 
107  return cell_area;
108 }
109 
110 /*!
111  * \brief Begin polygon area calculations.
112  *
113  * This initializes the polygon area calculation routines. It is used
114  * both for planimetric and latitude-longitude projections.
115  *
116  * \return 0 if the projection is not measurable (ie. imagery or xy)
117  * \return 1 if the projection is planimetric (ie. UTM or SP)
118  * \return 2 if the projection is non-planimetric (ie. latitude-longitude)
119  */
121 {
122  double a, e2;
123  double factor;
124 
125  if ((st->projection = G_projection()) == PROJECTION_LL) {
128  return 2;
129  }
131  if (factor > 0.0) {
132  st->units_to_meters_squared = factor * factor;
133  return 1;
134  }
135  st->units_to_meters_squared = 1.0;
136  return 0;
137 }
138 
139 /*!
140  * \brief Area in square meters of polygon.
141  *
142  * Returns the area in square meters of the polygon described by the
143  * <i>n</i> pairs of <i>x,y</i> coordinate vertices. It is used both for
144  * planimetric and latitude-longitude projections.
145  *
146  * You should call G_begin_polygon_area_calculations() function before
147  * calling this function.
148  *
149  * <b>Note:</b> If the database is planimetric with the non-meter grid,
150  * this routine performs the required unit conversion to produce square
151  * meters.
152  *
153  * \param x array of x coordinates
154  * \param y array of y coordinates
155  * \param n number of x,y coordinate pairs
156  *
157  * \return area in square meters of the polygon
158  */
159 double G_area_of_polygon(const double *x, const double *y, int n)
160 {
161  double area;
162 
163  if (st->projection == PROJECTION_LL)
164  area = G_ellipsoid_polygon_area(x, y, n);
165  else
166  area = G_planimetric_polygon_area(x, y, n) * st->units_to_meters_squared;
167 
168  return area;
169 }
double G_darea0_on_sphere(double lat)
Calculates integral for area between two latitudes.
Definition: area_sphere.c:48
void G_begin_zone_area_on_ellipsoid(double a, double e2, double s)
Begin area calculations for ellipsoid.
Definition: area_ellipse.c:47
double G_area_of_polygon(const double *x, const double *y, int n)
Area in square meters of polygon.
Definition: gis/area.c:159
2D/3D raster map header (used also for region)
Definition: gis.h:390
double G_area_of_cell_at_row(int row)
Cell area in specified row.
Definition: gis/area.c:87
int G_get_ellipsoid_parameters(double *a, double *e2)
get ellipsoid parameters
Definition: get_ellipse.c:67
double G_planimetric_polygon_area(const double *x, const double *y, int n)
Calculates planimetric polygon area.
Definition: area_poly2.c:25
#define x
double G_darea0_on_ellipsoid(double lat)
Calculate integral for area between two latitudes.
Definition: area_ellipse.c:63
struct state * st
Definition: parser.c:103
double north
Extent coordinates (north)
Definition: gis.h:436
int G_begin_polygon_area_calculations(void)
Begin polygon area calculations.
Definition: gis/area.c:120
int G_begin_cell_area_calculations(void)
Begin cell area calculations.
Definition: gis/area.c:46
void G_begin_zone_area_on_sphere(double r, double s)
Initialize calculations for sphere.
Definition: area_sphere.c:36
#define PROJECTION_LL
Projection code - Latitude-Longitude.
Definition: gis.h:99
double G_ellipsoid_polygon_area(const double *lon, const double *lat, int n)
Area of lat-long polygon.
Definition: area_poly1.c:129
void G_get_set_window(struct Cell_head *window)
Get the current working window (region)
double G_database_units_to_meters_factor(void)
Conversion to meters.
Definition: proj3.c:135
struct state state
Definition: parser.c:102
int G_projection(void)
Query cartographic projection.
Definition: proj1.c:32
void G_begin_ellipsoid_polygon_area(double a, double e2)
Begin area calculations.
Definition: area_poly1.c:61