[gdal-dev] Location change on gdalwarp reprojection

Wed Nov 11 11:32:37 PST 2020

Hey all,

I have some behaviour that I canâ€™t wrap my head around. When I reproject a tif using gdalwarp, the location (as visible in QGIS or after using gdal2tiles on a map) changes. It is visible when checking the bounds using QGIS.
It is also noticeable when checking the bounds with rasterio:

```
orig.tif
EPSG:27573
bounds [828662.0710931281, 3203193.9700861075, 829595.8710931281, 3204344.6700861077]
EPSG:27573 bounds [828662.0710931281, 3203193.9700861075, 829595.8710931281, 3204344.6700861077]
EPSG:3857 bounds [578066.9358574335, 5479808.142856716, 579420.0775091586, 5481457.907968456]
EPSG:4326 bounds [5.192863637248721, 44.09275933293886, 5.205019115521942, 44.103402376252255] <----- 1
gdalwarp -t_srs EPSG:3857 orig.tif repro.tif
repro.tif
EPSG:3857
bounds [608015.5294828439, 6272484.913981743, 609507.0864301398, 6274297.75045133]
EPSG:27573 bounds [829484.0038619096, 3748101.0103735547, 830512.7970793804, 3749330.248761157]
EPSG:3857 bounds [608015.5294828439, 6272484.913981743, 609507.0864301398, 6274297.75045133]
EPSG:4326 bounds [5.46189643116462, 48.98599229672266, 5.475295315193526, 48.99667806803407]  <----- 2
gdalwarp -t_srs EPSG:27573 repro.tif back.tif
back.tif
EPSG:27573
bounds [829484.0038619096, 3748101.034601565, 830512.7912119445, 3749330.2487611575]
EPSG:27573 bounds [829484.0038619096, 3748101.034601565, 830512.7912119445, 3749330.2487611575]
EPSG:3857 bounds [607947.003095218, 6272428.113526381, 609575.9130150724, 6274354.58877996]
EPSG:4326 bounds [5.461280848150926, 48.98565744919683, 5.475913594925499, 48.99701306471833]  <----- 3
```

The accompanying python code is three times like the following
```
input_file = folder + 'orig.tif'
img = rasterio.open(input_file)
print('orig.tif')
print(img.crs)
print('bounds', list(warp.transform_bounds(img.crs, img.crs, *img.bounds)))
print('EPSG:27573 bounds',list(warp.transform_bounds(img.crs, 'EPSG:27573', *img.bounds)))
print('EPSG:3857 bounds',list(warp.transform_bounds(img.crs, 'EPSG:3857', *img.bounds)))
print('EPSG:4326 bounds', list(warp.transform_bounds(img.crs, 'EPSG:4326', *img.bounds)))
```

As the arrows show, the bounds change between 1 and 2 (this is the unexpected behaviour for me). Yet they donâ€™t change when projecting back.

I am unsure if this is a bug in GDAL, something weird with this specific EPSG or the source tif. Any further steps for investigation or ideas would be very welcome

Many thanks

Evert