[gdal-dev] Location change on gdalwarp reprojection
jratike80
jukka.rahkonen at maanmittauslaitos.fi
Wed Nov 11 13:38:46 PST 2020
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://www.usgs.gov/media/images/landsat-1-8-landsatlook-image-examples 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%7C4ba1c63db9e74e36d0e008d88678923c%7Cfc89adff07ac47008853b7b7e906068e%7C0%7C0%7C637407199695587428%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=5vBljKWRNYXuXdTYgISUX8G5imyZ40Fbwd5Ai%2FA635E%3D&reserved=0
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at .osgeo
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
--
Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html
More information about the gdal-dev
mailing list