Loading...
Searching...
No Matches
area.c
Go to the documentation of this file.
/*!
* \file lib/gis/area.c
*
* \brief GIS Library - Area calculation functions.
*
* (C) 2001-2009 by the GRASS Development Team
*
* This program is free software under the GNU General Public License
* (>=v2). Read the file COPYING that comes with GRASS for details.
*
* \author Original author CERL
*/
#include <grass/gis.h>
static struct state {
struct Cell_head window;
double square_meters;
int projection;
double units_to_meters_squared;
/* these next are for lat-long only */
int next_row;
double north_value;
double north;
double (*darea0)(double);
} state;
static struct state *st = &state;
/*!
* \brief Begin cell area calculations.
*
* This routine must be called once before any call to
* G_area_of_cell_at_row(). It perform all initializations needed to
* do area calculations for grid cells, based on the current window
* "projection" field. It can be used in either planimetric
* projections or the latitude-longitude projection.
*
* \return 0 if the projection is not measurable (ie. imagery or xy)
* \return 1 if the projection is planimetric (ie. UTM or SP)
* \return 2 if the projection is non-planimetric (ie. latitude-longitude)
*/
int G_begin_cell_area_calculations(void)
{
double a, e2;
double factor;
G_get_set_window(&st->window);
switch (st->projection = st->window.proj) {
case PROJECTION_LL:
G_get_ellipsoid_parameters(&a, &e2);
if (e2) {
G_begin_zone_area_on_ellipsoid(a, e2, st->window.ew_res / 360.0);
st->darea0 = G_darea0_on_ellipsoid;
}
else {
G_begin_zone_area_on_sphere(a, st->window.ew_res / 360.0);
st->darea0 = G_darea0_on_sphere;
}
st->next_row = 0;
st->north = st->window.north;
st->north_value = st->darea0(st->north);
return 2;
default:
st->square_meters = st->window.ns_res * st->window.ew_res;
factor = G_database_units_to_meters_factor();
if (factor > 0.0)
st->square_meters *= (factor * factor);
return (factor > 0.0);
}
}
/*!
* \brief Cell area in specified row.
*
* This routine returns the area in square meters of a cell in the
* specified <i>row</i>. This value is constant for planimetric grids
* and varies with the row if the projection is latitude-longitude.
*
* \param row row number
*
* \return cell area
*/
double G_area_of_cell_at_row(int row)
{
double south_value;
double cell_area;
if (st->projection != PROJECTION_LL)
return st->square_meters;
if (row != st->next_row) {
st->north = st->window.north - row * st->window.ns_res;
st->north_value = st->darea0(st->north);
}
st->north -= st->window.ns_res;
south_value = st->darea0(st->north);
cell_area = st->north_value - south_value;
st->next_row = row + 1;
st->north_value = south_value;
return cell_area;
}
/*!
* \brief Begin polygon area calculations.
*
* This initializes the polygon area calculation routines. It is used
* both for planimetric and latitude-longitude projections.
*
* \return 0 if the projection is not measurable (ie. imagery or xy)
* \return 1 if the projection is planimetric (ie. UTM or SP)
* \return 2 if the projection is non-planimetric (ie. latitude-longitude)
*/
int G_begin_polygon_area_calculations(void)
{
double a, e2;
double factor;
if ((st->projection = G_projection()) == PROJECTION_LL) {
G_get_ellipsoid_parameters(&a, &e2);
G_begin_ellipsoid_polygon_area(a, e2);
return 2;
}
factor = G_database_units_to_meters_factor();
if (factor > 0.0) {
st->units_to_meters_squared = factor * factor;
return 1;
}
st->units_to_meters_squared = 1.0;
return 0;
}
/*!
* \brief Area in square meters of polygon.
*
* Returns the area in square meters of the polygon described by the
* <i>n</i> pairs of <i>x,y</i> coordinate vertices. It is used both for
* planimetric and latitude-longitude projections.
*
* You should call G_begin_polygon_area_calculations() function before
* calling this function.
*
* <b>Note:</b> If the database is planimetric with the non-meter grid,
* this routine performs the required unit conversion to produce square
* meters.
*
* \param x array of x coordinates
* \param y array of y coordinates
* \param n number of x,y coordinate pairs
*
* \return area in square meters of the polygon
*/
double G_area_of_polygon(const double *x, const double *y, int n)
{
double area;
if (st->projection == PROJECTION_LL)
area = G_ellipsoid_polygon_area(x, y, n);
else
area =
G_planimetric_polygon_area(x, y, n) * st->units_to_meters_squared;
return area;
}
G_area_of_cell_at_row
double G_area_of_cell_at_row(int row)
Cell area in specified row.
Definition area.c:86
G_begin_cell_area_calculations
int G_begin_cell_area_calculations(void)
Begin cell area calculations.
Definition area.c:45
G_area_of_polygon
double G_area_of_polygon(const double *x, const double *y, int n)
Area in square meters of polygon.
Definition area.c:158
G_begin_polygon_area_calculations
int G_begin_polygon_area_calculations(void)
Begin polygon area calculations.
Definition area.c:119
G_darea0_on_ellipsoid
double G_darea0_on_ellipsoid(double lat)
Calculate integral for area between two latitudes.
Definition area_ellipse.c:63
G_begin_zone_area_on_ellipsoid
void G_begin_zone_area_on_ellipsoid(double a, double e2, double s)
Begin area calculations for ellipsoid.
Definition area_ellipse.c:47
G_begin_ellipsoid_polygon_area
void G_begin_ellipsoid_polygon_area(double a, double e2)
Begin area calculations.
Definition area_poly1.c:65
G_ellipsoid_polygon_area
double G_ellipsoid_polygon_area(const double *lon, const double *lat, int n)
Area of lat-long polygon.
Definition area_poly1.c:133
G_planimetric_polygon_area
double G_planimetric_polygon_area(const double *x, const double *y, int n)
Calculates planimetric polygon area.
Definition area_poly2.c:25
G_begin_zone_area_on_sphere
void G_begin_zone_area_on_sphere(double r, double s)
Initialize calculations for sphere.
Definition area_sphere.c:35
G_darea0_on_sphere
double G_darea0_on_sphere(double lat)
Calculates integral for area between two latitudes.
Definition area_sphere.c:47
G_get_ellipsoid_parameters
int G_get_ellipsoid_parameters(double *a, double *e2)
get ellipsoid parameters
Definition get_ellipse.c:66
G_get_set_window
void G_get_set_window(struct Cell_head *window)
Get the current working window (region).
Definition gis/set_window.c:35
G_projection
int G_projection(void)
Query cartographic projection.
Definition proj1.c:32
G_database_units_to_meters_factor
double G_database_units_to_meters_factor(void)
Conversion to meters.
Definition proj3.c:146
x
#define x