Raster math with latitude-specific correction factor for each pixel

Question

I'm working on a problem that requires creating a new raster from a large (~100 GB) global basemap. A different correction factor must be applied for every degree latitude of the basemap, making an otherwise simple problem less straightforward. This has proven quite clunky and inefficient to do in ArcGIS Pro via Python scripts.

I have tried these two methods in ArcGIS Pro:

1. Converting the raster to points (one point per pixel) -> calculating the geometry attributes (lat) of each point -> applying the correction factor for each latitude accordingly and calculating the new value -> converting the points back to a raster. This method takes a very long time to calculate the latitude of each point, and the resulting point files are very large.

2. Splitting the raster into 0.5 degree latitude thick strips and applying the correction factor to each strip using the raster calculator / raster math functionality, then mosaicking the strips back together. This method seems to introduce slight numerical errors, as well as taking a long time.

Is there a faster way to do this processing using Python?

Answer 1

You could try xarray. It's a very useful, broad set of tools. Something along these lines should work:

import xarray as xr

raster = xr.open_rasterio(<your dataset>)
groupby = raster.groupby_bins(....)
calculated = groupby.apply(your_function)

Uou could also try iterating over a dimension