In my daily work, I am constantly asked to calculate areas of global raster datasets in geographic projection at 30 arc second resolution. These datasets are normally the result of a Combine operation (a typical example is a vegetation classes combined with country layer). To do this, our unit created a raster dataset with the area of each pixel in geographic projection at 30 arc second. With this area grid, a zonalstat is performed to sum the areas for each class. Since I am not sure how this area grid was created, I always wondered whether this approach is more accurate of just reprojecting the raster in a equal-area projection (from simple tests the results of the two methods are similar). Does anybody experienced a similar situation?
[GIS] More accurate way to calculate area of rasters
areacoordinate systemraster
Best Answer
There is a relatively simple exact formula for the area of any spherical quadrangle bounded by parallels (lines of latitude) and meridians (lines of longitude). It can be derived straightforwardly using basic properties of the ellipse (of major axis a and minor axis b) that is rotated around its minor axis to produce the ellipsoid. (The derivation makes a nice integral Calculus exercise but I believe would be of little interest on this site.)
The formula is simplified by breaking the calculation down into basic steps.
First, the distance between the east and west boundaries--the meridians l0 and l1--is a fraction of a whole circle equal to q = (l1 - l0)/360 (when the meridians are measured in degrees) or 1 = (l1 - l0) / (2*pi) (when the meridians are measured in radians). Find the area of the entire slice located between parallels f0 and f1 and just multiply that by q.
Second, we will employ a formula for the area of a horizontal slice of the ellipsoid bounded by the Equator (at f0=0) and a parallel at latitude f (=f1). The area of the slice between any two latitudes f0 and f1 (lying on the same hemisphere) will be the difference between the larger and smaller area.
Finally, provided the model is truly an ellipsoid (and not a sphere), the area of such a slice between the Equator and the parallel at latitude f is given by
where
aandbare the lengths of the major and minor axes of the generating ellipse, respectively,is its eccentricity, and
(This is much simpler than computing with geodesics, which are only approximations to the parallels anyway. Please note the comment by @cffk concerning a way to compute
log(zp/zm)in a way that avoids loss of precision at low latitudes.)area(f)is the area of the opaque slice from the equator up to latitude f (around 30 degrees north in the illustration. X and Y are geocentric Cartesian coordinate axes shown for reference.For the WGS 84 ellipsoid, use the constant values
entailing
(For a spherical model with a = b, the formula becomes indefinite. You have to take a limit as e --> 0 from above, which then reduces to the standard formula
2 * pi * a^2 * sin(f).)According to these formulas, a 30' by 30' quadrangle based on the Equator has an area of 3077.2300079129 square kilometers, while a 30' by 30' quadrangle touching a pole (which is really just a triangle) has an area of only 13.6086152 square kilometers.
As a check, the formulas applied to all cells of a 720 by 360 grid covering the earth's surface give a total surface area of 4 * pi * (6371.0071809)^2 square kilometers, indicating that the earth's authalic radius should be 6371.0071809 kilometers. This differs from the Wikipedia value only in the last significant figure (about a tenth of a millimeter). (I think Wikipedia's calculations are off a little :-).
As additional checks, I used versions of these formulas to reproduce Appendices 4 and 5 in Lev M. Bugayevskiy & John P. Snyder, Map Projections: A Reference Manual (Taylor & Francis, 1995). Appendix 4 shows arc lengths of 30'-long sections of meridians and parallels, given to the nearest meter. A spot check of results showed perfect agreement. I then recreated the table with 0.0005' increments, rather than 0.5' increments, and numerically integrated the quadrangle areas as estimated with these arc lengths. The total area of the ellipsoid was accurately reproduced to better than eight significant figures. Appendix 5 shows the values of
area(f)for f = 0, 1/2, 1, ..., 90 degrees, multiplied by 1/(2*pi). These values are given to the nearest square kilometer. A visual check of values near 0, 45 and 90 degrees showed perfect agreement.This exact formula can be applied using raster algebra beginning with a grid giving the latitudes of the upper limits of each cell and another giving the latitudes of the lower limits. Each of these is, essentially, a y-coordinate grid. (In each case you might want to create
sin(f)and thenzmandzpas intermediate results.) Subtract the two results, take the absolute value of that, and multiply by the fraction q obtained in the first step (equal to 0.5/360 = 1/720 for a 30' cell width, for instance). This will be a grid whose values contain the exact areas of each cell (up to the grid's own numerical precision). Just make sure to express the latitudes in the form expected by the sine function: many raster calculators will give you coordinates in degrees but expect radians for their trig functions!For the record, here are the exact areas of 30' by 30' cells on the WGS 84 ellipsoid from the Equator up to a pole, in intervals of 30', to 11 figures (the same number used for the minor radius b):
The values are in square kilometers.
If you would like to approximate these areas or simply understand their behavior better, the formula reduces to a power series following this pattern:
where
(An equivalent formula appears in Bugayevskiy & Snyder, op. cit., equation (2.1).)
Since e^2 is so small (around 1/150 for all ellipsoidal models of the earth) and z lies between 0 and 1, y is small, too. Thus the terms y^2, y^3, ... rapidly get smaller, adding over two more decimal places precision with each term. If we were to ignore y altogether, the formula would be that of the area of a sphere of radius b. The remaining terms can be understood as correcting for the earth's equatorial bulge.
Edit
Some questions have been raised concerning how a geodesic-distance calculation of area compares with these exact formulas. The geodesic distance method approximates each quadrangle by the geodesics, rather than parallels, that connect its corners horizontally, and applies the Euclidean formula for a trapezoid. For small quadrangles, such as 30' quads, this is biased slightly low and has a relative accuracy between 6 and 10 parts per million. Here is a plot of the error for WGS 84 (or any reasonable earth ellipsoid, for that matter):
Thus, if (1) you have easy access to geodesic distance calculations and (2) can tolerate ppm-level error, you might consider using those geodesic calculations and multiplying their results by 1.00000791 to correct the bias. For two more decimal places of precision, subtract pi/2 * cos(2f) / 10^6 from the correction factor: the result will be accurate to within 0.04 ppm.