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
G_begin_cell_area_calculations(void)46 int G_begin_cell_area_calculations(void)
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:
54 G_get_ellipsoid_parameters(&a, &e2);
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;
69 factor = G_database_units_to_meters_factor();
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 */
G_area_of_cell_at_row(int row)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 */
G_begin_polygon_area_calculations(void)120 int G_begin_polygon_area_calculations(void)
121 {
122 double a, e2;
123 double factor;
124
125 if ((st->projection = G_projection()) == PROJECTION_LL) {
126 G_get_ellipsoid_parameters(&a, &e2);
127 G_begin_ellipsoid_polygon_area(a, e2);
128 return 2;
129 }
130 factor = G_database_units_to_meters_factor();
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 */
G_area_of_polygon(const double * x,const double * y,int n)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 }
170