Raster analysis workflow in numpy

I like to do most raster analysis in Python for its flexible visualisation options and the ease of implementing more advanced algorithms than those supported by gdal_calc.

Usually this means working with multiple rasters covering the same extent: this is certainly true for gpxz where I’m merging multiple DEM rasters. But even when working with single rasters it can be helpful to load satellite imagery over the same extent.

My workflow looks like this:

  • Pick a reference CRS, resolution, and bounds (usually these will be taken from the highest-resolution raster, or they’ll be the desired parameters of the output).
  • Use gdalwarp to project, sample, and crop all rasters to the reference parameters.
  • Load each raster band into a numpy array using rasterio.
  • Now you have an array for each raster that are all the same shape, making it easy to do vectorised operations on multiple bands, or to apply the same function to all bands.
  • Arrays can be simply plotted with matplotlib.pyplot.imshow without worrying (too much) about projections.