Hi all, Many thanks for your help Jukka, much appreciated!
What you are saying is indeed very probably the cause. When using gdal_edit, the tif indeed remains on the location I expect it to be, after reprojecting and tiling. So this might actually also show bugs/weird behaviour in both QGIS and rasterio? - QGIS displaying the original file on the correct location -> it probably looks at the EPSG and not the full (kinda-off) WKT string then, correct? Should this be reported? - Rasterio reporting the EPSG as well wrongly, not sure if this is a confidence issue and if this should be reported as well? The issue is solved for me with this workaround, but if anyone on this list believes I should raise the above points somewhere else, let me know and I can take this up. Greetings On 12 Nov 2020, at 09:39, Rahkonen Jukka (MML) <[email protected]<mailto:[email protected]>> wrote: 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) <[email protected]<mailto:[email protected]>> Lähetetty: torstai 12. marraskuuta 2020 10.00 Vastaanottaja: Rahkonen Jukka (MML) <[email protected]<mailto:[email protected]>> 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 <[email protected]<mailto:[email protected]>> 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<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%7Cfeaf0ee3588b4e52d2b408d886e68280%7Cfc89adff07ac47008853b7b7e906068e%7C0%7C0%7C637407671880779985%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=wLYFUkTxcVuEAQ39eI0pG4ZDMvv8%2FCFilZOrrPBGymY%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 [email protected]<mailto:[email protected]> 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<https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Flists.osgeo.org%2Fmailman%2Flistinfo%2Fgdal-dev&data=04%7C01%7Cevert.etienne%40sitemark.com%7Cfeaf0ee3588b4e52d2b408d886e68280%7Cfc89adff07ac47008853b7b7e906068e%7C0%7C0%7C637407671880789981%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=GN7KZj2sVcynjI0PY9hLgKo3gY60ivKbne4HAs1kTtY%3D&reserved=0> _______________________________________________ gdal-dev mailing list [email protected]<mailto:[email protected]> 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<https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Flists.osgeo.org%2Fmailman%2Flistinfo%2Fgdal-dev&data=04%7C01%7Cevert.etienne%40sitemark.com%7Cfeaf0ee3588b4e52d2b408d886e68280%7Cfc89adff07ac47008853b7b7e906068e%7C0%7C0%7C637407671880789981%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=GN7KZj2sVcynjI0PY9hLgKo3gY60ivKbne4HAs1kTtY%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<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%7Cfeaf0ee3588b4e52d2b408d886e68280%7Cfc89adff07ac47008853b7b7e906068e%7C0%7C0%7C637407671880799978%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=J1bgclpj07UjPCNxgrKT5YFXE1rEyr7LUaeHBlNMx30%3D&reserved=0> _______________________________________________ gdal-dev mailing list [email protected]<mailto:[email protected]> 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<https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Flists.osgeo.org%2Fmailman%2Flistinfo%2Fgdal-dev&data=04%7C01%7Cevert.etienne%40sitemark.com%7Cfeaf0ee3588b4e52d2b408d886e68280%7Cfc89adff07ac47008853b7b7e906068e%7C0%7C0%7C637407671880799978%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=yVbmyGYQELLIeKqcepWFSvuSoSlrC9KRTW8%2BDIaN0t4%3D&reserved=0>
_______________________________________________ gdal-dev mailing list [email protected] https://lists.osgeo.org/mailman/listinfo/gdal-dev
