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.