How are calculate_default_transform and aligned_target different?


yon.davies@...
 

Greetings everyone,

I am trying to do reprojection in chunks and have been using rasterio.warp.calculate_default_transform for getting the projected extent and transform of the chunk. However, I noticed rasterio.warp.aligned_target also returns a transform and bounds, but the transform is a bit different from the one returned by calculate_default_transform. The documentation is a little vague – calculate_default_transform: Output dimensions and transform for a reprojection; aligned_target: Aligns target to specified resolution.

Can anyone illuminate this for me? And do you have any insight for which might be better for reprojecting a huge raster in chunks that will eventually be merged later? Thank you.


Sean Gillies
 


Hi Yon,

The intent is that rasterio.warp.reproject, given the values output from calculate_default_transform, will produce a result close to that of GDAL's gdalwarp utility *without* its "-tap" option. And that if you wrap the output of calculate_default_transform with aligned_target, adding a resolution, you will get a result close to that of gdalwarp *with* its "-tap" option.

You can certainly do your own warping in chunks, using Python concurrency features, but you may also rely on GDAL's warper for this. rasterio.warp.reproject will take rasterio datasets as inputs and output and can turn over all the details to GDAL. You can set the number of threads and a limit for memory utilization and GDAL's warper will handle the chunking.

Yours,

On Thu, Jun 10, 2021 at 3:12 PM <yon.davies@...> wrote:
Greetings everyone,

I am trying to do reprojection in chunks and have been using rasterio.warp.calculate_default_transform for getting the projected extent and transform of the chunk. However, I noticed rasterio.warp.aligned_target also returns a transform and bounds, but the transform is a bit different from the one returned by calculate_default_transform. The documentation is a little vague – calculate_default_transform: Output dimensions and transform for a reprojection; aligned_target: Aligns target to specified resolution.

Can anyone illuminate this for me? And do you have any insight for which might be better for reprojecting a huge raster in chunks that will eventually be merged later? Thank you.


--
Sean Gillies


amaury.dehecq@...
 

Hi all,

on the same topic, is there an easy way to get the destination parameters (transform, width, height) exactly as GDAL would do, considering the user might provide any of, or a combination of, the following: dst_bounds, dst_width/height, dst_res.

I understand that if dst_width/height or resolution are provided, the output transformation and width/height can easily be calculated using calculate_default_transform (and aligned_target if the resolution is to be strictly enforced). But what if dst_bounds are also set?

Let's consider the following cases:

- dst_bounds is set alone -> GDAL output strictly enforces the output bounds but keeps the output resolution as close as possible to the input resolution

- dst_bounds and res are set -> GDAL output enforces the output resolution and use the bounds that are as close as possible to the requested bounds

- dst_bounds and res are set+ tap option -> same as before but output bounds is a multiple of the resolution

- dst_bounds and width/height are set -> GDAL outputs shape and bounds are strictly enforced, the resolution is set accordingly.

Currently, I cannot find an easy way to do this with rasterio. At least a specific, non trivial order of rasterio functions have to be called, to get exactly this behavior. There could possibly be a specific rasterio function to do this.

Thanks in advance!

Amaury

 


Sean Gillies
 

Hi Amaury,

Rasterio's calculate_default_transform is only a thin wrapper around GDALSuggestedWarpOutput2: https://gdal.org/api/gdal_alg.html#_CPPv424GDALSuggestedWarpOutput212GDALDatasetH19GDALTransformerFuncPvPdPiPiPdi. The behavior you see in GDAL is implemented in the gdalwarp program and is only partially implemented in rasterio's rio-warp command. So, to get the same behavior, you would need to port some logic from C (the gdalwap_lib.cpp file, if I remember correctly) to Python.

We'd like to change rasterio's warp API to align better with GDAL https://github.com/mapbox/rasterio/issues/1990 but no work on this has started yet.

On Tue, Jul 6, 2021 at 4:20 AM <amaury.dehecq@...> wrote:

Hi all,

on the same topic, is there an easy way to get the destination parameters (transform, width, height) exactly as GDAL would do, considering the user might provide any of, or a combination of, the following: dst_bounds, dst_width/height, dst_res.

I understand that if dst_width/height or resolution are provided, the output transformation and width/height can easily be calculated using calculate_default_transform (and aligned_target if the resolution is to be strictly enforced). But what if dst_bounds are also set?

Let's consider the following cases:

- dst_bounds is set alone -> GDAL output strictly enforces the output bounds but keeps the output resolution as close as possible to the input resolution

- dst_bounds and res are set -> GDAL output enforces the output resolution and use the bounds that are as close as possible to the requested bounds

- dst_bounds and res are set+ tap option -> same as before but output bounds is a multiple of the resolution

- dst_bounds and width/height are set -> GDAL outputs shape and bounds are strictly enforced, the resolution is set accordingly.

Currently, I cannot find an easy way to do this with rasterio. At least a specific, non trivial order of rasterio functions have to be called, to get exactly this behavior. There could possibly be a specific rasterio function to do this.

Thanks in advance!

Amaury

 

_,_._,_



--
Sean Gillies


amaury.dehecq@...
 

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