Land surface area including elevation

23.08.2020 | Niel Wagensommer

The question

On this episode of the "A Problem Squared" podcast, a viewer asked Bec Hill and Matt Parker to investigate how the estimated surface area of a country changes when taking into account topography. This turned out to be a much harder problem than expected, as he explains in this video:

If you want to figure out the actual area when taking into account topography, you run into the fractal problem. As Matt explains, the higher your resolution gets, the more distorting features you get in your data. First mountains, then gullies and at some point ridges between stones. So the results I figured out are all relative to a resolution.

The data

If you are looking for any kind of geographical data, NASA has you covered. Unfortunately, those datasets can get quite large. The ALOS_PALSAR_RTC_HIGH_RES with 1 arcsecond resolution and worldwide coverage is a nice 119.1TB. For our purposes, we'll stick with the GEBCO dataset which you can get as a nice simply 21600x10800 PNG file. The coordinates for each pixel represent latitude and longitude. The brightness between 0-255 represents heights from 0m to 6400m. This gives us a resolution at the equator of 40075km/21600px = 1.85km/px horizontally and 6400m/255 = 25m vertically. At 60° latitude, resolution in longitude gets even better at 40075km*cos(60°)/21600px = 0.93km/px. Resolution in latitude stays constant at any position.

The process

To find out how much topography distorts surface area, I estimated the local distortion around each pixel. First the code constructs a region of 5x5 samples around the current pixel. These samples are then transformed into 3D positions based on latitude, longitude and elevation assuming a round Earth - once at sea level and once including topography from the dataset. The surface area of this patch is approximated by turning it into a mesh of 16 quadrilaterals split into a total of 32 triangles and taking the sum of their surface area. I did this to avoid calculating the actual surface ares which leads down a rabbit hole of mathematical questions, and I don't need to be dealing with Chebyshev polynomial expansion on my day off. The result is good enough. Once you have the results for the flat and distorted patch, finding out the rate of distortion is as easy as dividing one by the other.

The result

I ran the whole process at 100%, 20% and 10% resolution versions of the base image, roughly 2km, 10km and 20km. All the code can be found on GitHub. Here are the links to the results at full resolution: 10% 20% 100%

Some observations: