Exploring the estimation of areas on the Earth’s surface.

projection

python

pydeck

Author

Tyler Erickson

Published

August 17, 2024

In recent years, there has been an explosion of high-resolution raster datasets that characterize the Earth’s surface. An example of this is the Global map of Forest Cover for year 2020 (GFC2020).

The GFC2020 dataset is described as serving to “define the extent of tree cover for year 2020 at 10 m spatial resolution”^{1}, and it shows variations in forested land in great detail. However, using the dataset to accurately answer a seemingly simple question can be quite challenging. For example:

What is the total area of forest in Iceland?

To determine the total area, can’t we just simply count up the forest pixels and multiply by 10m squared?

\[A_{forest} = n_{forest pixels} \cdot (10m)^2 \]

Short answer:

Nope.

Longer answer: (click to expand)

There are at least two problems with that approach.

The GFC2020 dataset uses a coordinate reference system (CRS) of EPSG:4326. This is a fairly common CRS used by global datasets. EPSG:4326 has axes of geodetic latitude and geodetic longitude measured in units of degrees, which is an angle measurement. Typically, the area of something on the Earth’s surface is not measured in unit of degrees squared, but rather a unit of length squared^{2}.

With a little sleuthing^{3} it can be determined that the unit of measure is an angle in degrees:

AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],

and that the pixel size is much smaller than a degree:

The surface area of a sphere (or ellipsoid) bulges out, so it is greater than a planar area between the same boundaries. The difference becomes more signficant as the curvature increases. For more details, read up on spherical excess.

What about using a library to reproject the raster image to an equal-area projection? That may be a viable solution, but the data will be altered by resampling (usually nearest neighbor resampling for categorical data), which changes the resulting area calculation^{4}. And it may be challenging from a computation and storage standpoint, depending on the size of dataset. And you will need to install a projection library, which a complication that I would like to avoid.

So how can you go about accurately determining the areas of regions on the Earth’s surface? That’s what this post will cover, at least for a certain class of rectangular areas. We’ll start with the relevant mathematical equations and show how they can be implemented with standard Python libraries.

Note: Alternative options

There are many software libraries and tools that can help with calculating areas on the Earth. Many are based on the open source PROJ software, which provides generic coordinate transformation functionality via both command line applications and an API to power higher-level tools such as GDAL, pyproj, Rasterio, GeoPandas, PostGIS, QGIS, and many, many more… use them if you can! These tools are far more flexible (and rigorously tested) than the approach I describe here. However, that flexibility results in a complicated codebase which few people understand (and can contribute to). Even if open source, they are black box systems for most people. This post attempts to transparently explain the “how”, but for a very narrow slice of functionality provided by PROJ-based tools.

Since calculating for a sphere is more straightforward than calculations for an ellipsoid, we solve for the radius of a spherical model that has the same surface area as the ellipsoid model. This radius is called the authalic radius (\(R_2\))

For area segments between two bands of constant latitudes (\(\phi_1\) and \(\phi_2\)), the Spherical segment formula can be used. The curved surface area of the spherical zone between two latitudes is given by: \[
A_{\text{zone}} = 2 \pi R_2 h
\]

where \(h\) is the height of the segment: \[
h = R_2 (\sin \phi_2 - \sin \phi_1)
\tag{1}\]

Constraining the area to be a patch (\(A_{\text{patch}}\)) between longitudes (\(\lambda_1 <= \lambda <= \lambda_2\)) gives:

print(f'{(patch_area_single_rectangle(R=R_auth, lon1=-180, lon2=180, lat1=-90, lat2=90)*1e-6):.3e} km^2 is the area of the globe')print(f'{(patch_area_single_rectangle(R=R_auth, lon1=0, lon2=180, lat1=-90, lat2=90)*1e-6):.3e} km^2 is the area of a (eastern) hemisphere')print(f'{(patch_area_single_rectangle(R=R_auth, lon1=0, lon2=10, lat1=-90, lat2=90)*1e-6):.3e} km^2 is the area of a 10 degree longitude band')print(f'{(patch_area_single_rectangle(R=R_auth, lon1=0, lon2=10, lat1=0, lat2=10)*1e-6):.3e} km^2 is the area of a 10 x 10 degree patch')

5.101e+08 km^2 is the area of the globe
2.550e+08 km^2 is the area of a (eastern) hemisphere
1.417e+07 km^2 is the area of a 10 degree longitude band
1.230e+06 km^2 is the area of a 10 x 10 degree patch

This works for single patches, but we can adapt the function to efficiently estimate a (potentially irregular) grid of patches by accepting lists of latitude edges and longitudes edges:

The globe is interactive. Roll your cursor over the patches to see the calculation results. Click and drag to rotate the globe.

Wrapping up

In this post I’ve presented an approach for calculating patch areas on the Earth’s surface. In a future post I’ll expand on this to calculate the areas of every pixel in a projected image.

Global map of Forest Cover for year 2020, produced by the European Union Joint Research Centre. Click here to interactively explore dataset on the live website.An ellipsoid / spheroid. Source: Wikimedia CommonsAn ellipsoid / spheroid. Source: Wikimedia CommonsAn ellipsoid / spheroid. Source: Wikimedia Commons

Footnotes

I really dislike when the term “resolution” is used instead of pixel size, but I’ll leave that discussion for another time…↩︎

Strangely, neither the website nor the publication describing the dataset list the actual pixel size. To determine it, I ended up using the command line program gdalinfo to inspect a sample tile obtained via the download page.↩︎

The sensitivity of accuracy to reprojection might be worth exploring in a future post.↩︎

a.k.a. the lat/lon projection or unprojected projection↩︎

the de facto standard for web mapping applications↩︎