[GIS] How to calculate the bounding box in projected coordinates

coordinate systemgdal

I'm using GDAL in Python to reproject rasters. I'd like to be able to take an existing raster, and project it into a new raster to my projection and output pixel size of choice.

I'm worried about how to calculate the bounding box in projected coordinates that I will in turn base the geotransform off of. gdal has a function called gdal.AutoCreateWarpedVRT that apparently creates a memory raster that I can then directly read the geotransform from. But I sometimes have VERY large rasters that don't even fit in main memory, I'd hate to create a memory based raster just to get the geotransform.

Is it sufficient to transform the vertices of the original bounding box into the projected coordinate system and then somehow make that box a square? (I'm sure the box will be distorted.)

Best Answer

A safe way to get the exact new extent is to

  1. create a polygon from the raster extent with gdaltindex
  2. densify the polygon geometry
  3. reproject the polygon to the target CRS with ogr2ogr
  4. calculate the extent in the new CRS
Related Question