[gdal-dev] Location change on gdalwarp reprojection
Evert Etienne (SITEMARK)
evert.etienne at sitemark.com
Wed Nov 11 11:46:18 PST 2020
It does not happen when just warping, but it does occur when warping to the same EPSG as can be seen in the follow logs
```
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] <----- OK
gdalwarp orig.tif justwarp.tif
justwarp.tif
EPSG:27573
bounds [828662.0710931281, 3203193.9700861075, 829595.8710931282, 3204344.6700861077]
EPSG:27573 bounds [828662.0710931281, 3203193.9700861075, 829595.8710931282, 3204344.6700861077]
EPSG:3857 bounds [578066.9358574335, 5479808.142856716, 579420.0775091586, 5481457.907968456]
EPSG:4326 bounds [5.192863637248721, 44.09275933293886, 5.205019115521945, 44.103402376252255] <----- SAME
gdalwarp -t_srs EPSG:27573 orig.tif noreproj.tif
noreproj.tif
EPSG:27573
bounds [829527.7905651039, 3748136.6066527413, 830468.8074871361, 3749294.6346238866]
EPSG:27573 bounds [829527.7905651039, 3748136.6066527413, 830468.8074871361, 3749294.6346238866]
EPSG:3857 bounds [608015.3560148155, 6272484.6841915725, 609507.2374494626, 6274297.96416575]
EPSG:4326 bounds [5.461894872874811, 48.98599094207838, 5.475296671823179, 48.99667932763942] <----- CHANGED
```
GDAL info of orig.tif:
```
Driver: GTiff/GeoTIFF
Files: orig.tif
Size is 18676, 23014
Coordinate System is:
PROJCRS["NTF (Paris) / Lambert zone III",
BASEGEOGCRS["NTF (Paris)",
DATUM["Nouvelle Triangulation Francaise (Paris)",
ELLIPSOID["Clarke 1880 (IGN)",6378249.2,293.466021293627,
LENGTHUNIT["metre",1]]],
PRIMEM["Paris",2.5969213,
ANGLEUNIT["grad",0.0157079632679489]],
ID["EPSG",4807]],
CONVERSION["unnamed",
METHOD["Lambert Conic Conformal (1SP)",
ID["EPSG",9801]],
PARAMETER["Latitude of natural origin",54.4444444444445,
ANGLEUNIT["grad",0.0157079632679489],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",0,
ANGLEUNIT["grad",0.0157079632679489],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.999877499,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",600000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",3200000,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing",north,
ORDER[2],
LENGTHUNIT["metre",1]],
ID["EPSG",27573]]
Data axis to CRS axis mapping: 1,2
Origin = (828662.071093128062785,3204344.670086107682437)
Pixel Size = (0.050000000000002,-0.050000000000008)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 828662.071, 3204344.670) ( 3d 7'33.28"E, 48d59'48.23"N)
Lower Left ( 828662.071, 3203193.970) ( 3d 7'30.95"E, 48d59'11.01"N)
Upper Right ( 829595.871, 3204344.670) ( 3d 8'19.19"E, 48d59'46.98"N)
Lower Right ( 829595.871, 3203193.970) ( 3d 8'16.85"E, 48d59' 9.76"N)
Center ( 829128.971, 3203769.320) ( 3d 7'55.07"E, 48d59'28.99"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
Overviews: 9338x11507, 4669x5754, 2335x2877, 1168x1439, 584x720, 292x360, 146x180
Mask Flags: PER_DATASET ALPHA
Overviews of mask band: 9338x11507, 4669x5754, 2335x2877, 1168x1439, 584x720, 292x360, 146x180
Unit Type: metre
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
Overviews: 9338x11507, 4669x5754, 2335x2877, 1168x1439, 584x720, 292x360, 146x180
Mask Flags: PER_DATASET ALPHA
Overviews of mask band: 9338x11507, 4669x5754, 2335x2877, 1168x1439, 584x720, 292x360, 146x180
Unit Type: metre
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
Overviews: 9338x11507, 4669x5754, 2335x2877, 1168x1439, 584x720, 292x360, 146x180
Mask Flags: PER_DATASET ALPHA
Overviews of mask band: 9338x11507, 4669x5754, 2335x2877, 1168x1439, 584x720, 292x360, 146x180
Unit Type: metre
Band 4 Block=256x256 Type=Byte, ColorInterp=Alpha
Overviews: 9338x11507, 4669x5754, 2335x2877, 1168x1439, 584x720, 292x360, 146x180
Unit Type: metre
```
> On 11 Nov 2020, at 20:32, Evert Etienne (SITEMARK) <evert.etienne at sitemark.com> wrote:
>
>
> 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 at lists.osgeo.org
> https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Flists.osgeo.org%2Fmailman%2Flistinfo%2Fgdal-dev&data=04%7C01%7Cevert.etienne%40sitemark.com%7C4ba1c63db9e74e36d0e008d88678923c%7Cfc89adff07ac47008853b7b7e906068e%7C0%7C0%7C637407199695587428%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=5vBljKWRNYXuXdTYgISUX8G5imyZ40Fbwd5Ai%2FA635E%3D&reserved=0
More information about the gdal-dev
mailing list