[gdal-dev] Location change on gdalwarp reprojection
Evert Etienne (SITEMARK)
evert.etienne at sitemark.com
Thu Nov 12 00:02:41 PST 2020
Hi
Thanks for your response Jukka!
Your first explanation about the bounds make sense, but shouldn’t they still be in the same area (overlapping)?
When specifying the source epsg as EPSG:27573, there is indeed no changes anymore as visible in the following logs:
(Am I correct for thinking this is actually just equivalent to overwriting the CRS included in the tif?)
```
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]
gdalwarp -s_srs EPSG:27573 -t_srs EPSG:27573 orig.tif from_27573_into_27573.tif
from_27573_into_27573.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]
```
So you’re suggesting the included CRS in orig.tif does not match EPSG:27573 perfectly then? Would this be a mismatch in the WKT string then?
I am still a bit lost how this would cause that in QGIS the orig.tif and repro.tif show up in a totally different location where orig.tif is where I would expect it to be.
And as well why rasterio reads the CRS as being EPSG:27573 (see logs), but that might be a separate issue?
Many thanks
On 11 Nov 2020, at 22:38, jratike80 <jukka.rahkonen at maanmittauslaitos.fi<mailto:jukka.rahkonen at maanmittauslaitos.fi>> wrote:
Hi,
The bounds do not mean exactly what you think. Re-projecting a rectangular
image from EPSG:27573 into EPSG:4326 rotates the image somewhat
counter-clockwise. This image is from another context
https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.usgs.gov%2Fmedia%2Fimages%2Flandsat-1-8-landsatlook-image-examples&data=04%7C01%7Cevert.etienne%40sitemark.com%7C3300bad68a154491a2d308d8868a300f%7Cfc89adff07ac47008853b7b7e906068e%7C0%7C0%7C637407275358259724%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=IIj2gHAu18%2Fd%2F%2BJyxHQxutxFg%2F0cqzNwSgIhNkEyNF8%3D&reserved=0 but
you can get the idea. The new bounds report the min-max coordinates of the
envelope that contains that rotated image. The corners of the original image
are not at the corners of the warped image. That explains why WGS 84 bounds
in EPSG:3857 are different.
The difference when you warp to EPSG:27573 may mean that the coordinate
system of the original image is not EPSG:27573, at least not for the Proj
library. What happens with command:
gdalwarp -s_srs EPSG:27573 -t_srs EPSG:27573 orig.tif
from_27573_into_27573.tif
-Jukka Rahkonen-
Evert Etienne (SITEMARK) wrote
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@
> 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 .osgeo
https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Flists.osgeo.org%2Fmailman%2Flistinfo%2Fgdal-dev&data=04%7C01%7Cevert.etienne%40sitemark.com%7C3300bad68a154491a2d308d8868a300f%7Cfc89adff07ac47008853b7b7e906068e%7C0%7C0%7C637407275358259724%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=WcvmUQ4g9%2FEBwzFxF4d4vWX3OlQwQmVvyNUAVCvrj4A%3D&reserved=0
_______________________________________________
gdal-dev mailing list
gdal-dev at .osgeo
https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Flists.osgeo.org%2Fmailman%2Flistinfo%2Fgdal-dev&data=04%7C01%7Cevert.etienne%40sitemark.com%7C3300bad68a154491a2d308d8868a300f%7Cfc89adff07ac47008853b7b7e906068e%7C0%7C0%7C637407275358259724%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=WcvmUQ4g9%2FEBwzFxF4d4vWX3OlQwQmVvyNUAVCvrj4A%3D&reserved=0
--
Sent from: https://eur03.safelinks.protection.outlook.com/?url=http%3A%2F%2Fosgeo-org.1560.x6.nabble.com%2FGDAL-Dev-f3742093.html&data=04%7C01%7Cevert.etienne%40sitemark.com%7C3300bad68a154491a2d308d8868a300f%7Cfc89adff07ac47008853b7b7e906068e%7C0%7C0%7C637407275358259724%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=3FP%2Fo4GIWEzFS9CUz0cjEVbT02XPgKCoTnI%2F6EC%2Bhzw%3D&reserved=0
_______________________________________________
gdal-dev mailing list
gdal-dev at lists.osgeo.org<mailto: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%7C3300bad68a154491a2d308d8868a300f%7Cfc89adff07ac47008853b7b7e906068e%7C0%7C0%7C637407275358259724%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=WcvmUQ4g9%2FEBwzFxF4d4vWX3OlQwQmVvyNUAVCvrj4A%3D&reserved=0
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20201112/173e148d/attachment-0001.html>
More information about the gdal-dev
mailing list