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?
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
area(f) = pi * b^2 * (log(zp/zm) / (2*e) + sin(f) / (zp*zm))
b are the lengths of the major and minor axes of the generating ellipse, respectively,
e = sqrt(1 - (b/a)^2)
is its eccentricity, and
zm = 1 - e*sin(f); zp = 1 + e*sin(f)
(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
a = 6 378 137 meters, b = 6 356 752.3142 meters,
e = 0.08181919084296
(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 then
zp as 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:
area(f) = 2 * pi * b^2 * z * (1 + (4/3)y + (6/5)y^2 + (8/7)y^3 + ...)
z = sin(f), y = (e*z)^2.
(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.
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.
The answer to radouxju's question depends on the shape of the pixel when projected onto the ellipsoid. If the coordinate system of the raster is longitude and latitude, then the pixel is a rhumb line rectangle and whuber's answer can be used, or, more generally, you can use the formula for a polygon whose edges are rhumb lines. If the coordinate system is a large scale conformal projection (UTM, state plane, etc.), it would be more accurate to approximate the edges by geodesics and use the formula for a geodesic polygon. Geodesic polygons are probably best for general use, since, unlike rhumb line polygons, they are "well behaved" close to the poles.
Implementations of the formulas for geodesic and rhumb line polygons are provided by my library GeographicLib. The geodesic area is available in several languages; the rhumb line area is C++ only. There's a online version (geodesic + rhumb line) available here. The accuracy of these calculations is typically better than 0.1 square meter.
You'll have to judge on credible/official... The geodesic formulas are derived in The area under the geodesic (Danielsen, 1989, subscription required), and Algorithms for geodesics (Karney, 2013, open access). The rhumb line formulas are given here.
I ran across this question when trying to determine a formula for the area of a WGS84 pixel. While @whuber's answer does contain this information, it was still some work to get a formula for the area of a square degree pixel at a given latitude. I've included a Python function I wrote below that abstracts this into a single call. While it does not directly answer the poster's question about the area of an ENTIRE raster (though one could sum the areas of all the pixels), I think it is still useful information for someone who may be searching for a similar calculation.
def area_of_pixel(pixel_size, center_lat): """Calculate m^2 area of a wgs84 square pixel. Adapted from: https://gis.stackexchange.com/a/127327/2397 Parameters: pixel_size (float): length of side of pixel in degrees. center_lat (float): latitude of the center of the pixel. Note this value +/- half the `pixel-size` must not exceed 90/-90 degrees latitude or an invalid area will be calculated. Returns: Area of square pixel of side length `pixel_size` centered at `center_lat` in m^2. """ a = 6378137 # meters b = 6356752.3142 # meters e = math.sqrt(1 - (b/a)**2) area_list =  for f in [center_lat+pixel_size/2, center_lat-pixel_size/2]: zm = 1 - e*math.sin(math.radians(f)) zp = 1 + e*math.sin(math.radians(f)) area_list.append( math.pi * b**2 * ( math.log(zp/zm) / (2*e) + math.sin(math.radians(f)) / (zp*zm))) return pixel_size / 360. * (area_list - area_list)