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
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Reply via email to