[gdal-dev] Location change on gdalwarp reprojection

Thu Nov 12 00:39:42 PST 2020

Hi,

Yes, envelopes should grow every time when the raster is rotated, but they should overlap. If they move like they do in your case, something seems to go wrong. I can only say that obviously Proj is seeing some difference between the projection that is stored into your GeoTIFF and what it thinks that EPSG:27573 means. By comparing the WKT text from your mail with the printout of â€œprojinfo epsg:27573â€ (with Proj version 6.3.2) I can see these differences:

Image
--------
CONVERSION["unnamed",
PARAMETER["Latitude of natural origin",54.4444444444445,

Proj
-----
CONVERSION["Lambert zone III",
PARAMETER["Latitude of natural origin",49,

Experts on the French coordinate systems and Proj may continue from this, I fear I canâ€™t help more than this. And you are right that using â€œ-s_srsâ€ in gdalwarp overrides the projection that is stored into GeoTIFF tags. If the result is good you can use that as a workaround. It is also possible to update the geotiff tags with â€œgdal_edit -a_srs epsg:27573â€ but before that you should be sure that it is the right thing to do.

-Jukka Rahkonen-

LÃ¤hettÃ¤jÃ¤: Evert Etienne (SITEMARK) <evert.etienne at sitemark.com>
LÃ¤hetetty: torstai 12. marraskuuta 2020 10.00
Vastaanottaja: Rahkonen Jukka (MML) <jukka.rahkonen at maanmittauslaitos.fi>
Aihe: Re: [gdal-dev] Location change on gdalwarp reprojection

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<mailto: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<mailto: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/43b6ef96/attachment-0001.html>