Hi Sean,

Ok thanks for the clarifications. I think I found GDAL's logic around here: https://github.com/OSGeo/gdal/blob/d01976cd4fc11144662f44ec4e658ce1871f1713/gdal/apps/gdalwarp_lib.cpp#L3105

It would indeed be useful to have a function that returns the best transform and shape given any bounds, resolution and shape. That's what I tried to do here in a Python implementation: https://github.com/GlacioHack/GeoUtils/blob/12bd607fe8024720ff264545abf0cc5ba0156ea6/geoutils/georaster.py#L879 (until L 932), if that can be of any help. In the meantime, I think I'll just keep my current implementation, which seems to provide similar results to what gdalwarp would do (I'll check the code in detail later).

Amaury