I am trying to estimate how much area of emerged landmasses is covered by each of the 12 great groups of soils defined in the USDA Soil Taxonomy with rasterio/rioxarray. My results differ quite a bit from previous attempts and I would like to know if anyone can double check my approach. I am suspicious that my approach to introduce minimum distortion in area estimates might be flawed.

The dataset's original projection is EPSG:4326 with a resolution of ~250 m. I first reproject it to an equal-earth projection (i.e., EPSG:8857) using rio.reproject for a given resolution res using the nearest resampling method, then I count unique values with np.unique, I extract the number of pixels for the class of interest and I multiply that by res**2.

# Load data
soilsrc = rio.open_rasterio('data.tif')
# Set resolution 
res = 1000
# Set classes of interest
andIdx = [50, 58, 59, 61,63,64,74,75,76,77,80]
# Reprojection
soilsrcP = soilsrc.rio.reproject(CRS.from_epsg(8857), resolution=res, resampling=Resampling.nearest)
# Get unique values
unique, count = np.unique(soilsrcP, return_counts = True)
# convert to a DataFrame
cts = pd.DataFrame()
cts['values'] = unique
cts['counts'] = count
# Get the area of the class of interest
andosols = cts[cts['values'].isin(andIdx)]['counts']*res**2

Does that make sense conceptually? Besides the resampling bit (which I do for now so I can test the code on my laptop), am I introducing any unnecessary distortion? Can anyone suggest a better way to preserve the area, preferably using Python? I don't see a polygonisation as an option as I think this global dataset is too large for that.

1 Answer 1


Conceptually at least, I think what you have here is ok. You are basically finding the area of a given pixel and multiplying it by the number of pixels in that class. Now, this can get a little more complicated for global models as pixel area changes with latitude for geographic projections and reprojecting to a projected system won't be perfect, but an equal area is a good choice.

You can always double check by manually transforming the pixel area in degrees to m and summing those (using code from rsgislib):

ellipse = [6378137.0, 6356752.314245]

radlat = numpy.deg2rad(latitude)

Rsq = (ellipse[0] * numpy.cos(radlat)) ** 2 + (ellipse[1] * numpy.sin(radlat)) ** 2
Mlat = (ellipse[0] * ellipse[1]) ** 2 / (Rsq**1.5)
Nlon = ellipse[0] ** 2 / numpy.sqrt(Rsq)
x_size = numpy.pi / 180 * numpy.cos(radlat) * Nlon * lon_size
y_size = numpy.pi / 180 * Mlat * lat_size

x_size, y_size

However, you will first need an array of the image pixel latitude which you can work out using the geotransform and image size

  • Thank loads for double checking my sanity. I am still puzzled by the lack of sensitivity/error analysis in the majority of the literature using global data that I am familiar with. Do you know of any approach/textbooks that treat on assessing the uncertainties associated with these reprojections? Surely this must have been assessed somewhere...
    – e5k
    Nov 28, 2022 at 9:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.