search

OSGeo / gdal

GDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.
https://gdal.org
Other
4.8k stars 2.51k forks source link

gdalwarp - warped image with unexpected small shift #10510

Closed jacmendt closed 1 month ago

jacmendt commented 1 month ago

What is the bug?

We have a couple of historic images with the coordinate system EPSG:4314. We publish the images via WMS as well as TMS. For the TMS we use gdal2tiles. Recently we did a reprocessing of some of the TMS and observe now shifts for historic maps of poland, which were not there before. On the image below you can see the shifts on the edges of the map, were it does not fit the roads.

pic_1

To track down our issues, we start to first warp the image via gdalwarp (GDAL-Version: 3.4.1 | Proj-Bin: 8.2.1-1 | Ubuntu: 22.04.4) from EPSG:4314 to EPSG:3857. After that we displayed the results in QGIS (in EPSG:3857). We could again see the shift.

gdalwarp -t_srs EPSG:3857 historic_map_4314.tif historic_map_3857_standard.tif

pic_2

Because the warping was working in the past, I tried it with an older gdal version.

docker run --rm -v $(pwd):/data osgeo/gdal:alpine-normal-v2.4.1 gdalwarp -t_srs EPSG:3857 /data/historic_map_4314.tif /data/historic_map_3857_docker.tif

pic_3

I also detected that somehow, when I load in the historic_map_4314.tif directly into QGIS, QGIS does some other kind of transformation and the image again is displayed properly.

pic_4 pic_5

I tried to use the QGIS transformation directly with gdalwarp, but this leads me to the following error:

gdalwarp -ct '+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=de_adv_BETA2007.tif +step +proj=webmerc +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84' historic_map_4314.tif historic_map_3857_ct.tif
ERROR 1: Too many points (9801 out of 9801) failed to transform, unable to compute output bounds.

I'm not an expert in coordinate systems, but according to [https://epsg.io/4314] there are different transformations. By selection from this page one of the more specific transformation, I can also produce a proper result:

gdalwarp -s_srs '+proj=longlat +ellps=bessel +towgs84=612.4,77,440.2,-0.054,0.057,-2.797,2.55 +no_defs +type=crs' -t_srs EPSG:3857 historic_map_4314.tif historic_map_3857_specific.tif

pic_6

I'm not sure if this a bug or just something with the change of definition of the coordinate system. We already have found some ways to work around it, so do not hesitate to close the issue, if you think this is not a GDAL issue. If it is a regression in the GDAL library we would also be willing to help fund the fix. Of course, we are also happy to receive feedback in order to better address the topic.

Steps to reproduce the issue

With gdalversion 3.4.1 and Ubuntu 22.04 use the following command:

gdalwarp -t_srs EPSG:3857 historic_map_4314.tif historic_map_3857_standard.tif

Test-Data: historic_map_4314.zip

Versions and provenance

GDAL: 3.4.1 Ubuntu: 22.04

Additional context

No response

rouault commented 1 month ago

I tried to use the QGIS transformation directly with gdalwarp, but this leads me to the following error:

gdalwarp -ct '+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=de_adv_BETA2007.tif +step +proj=webmerc +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84' historic_map_4314.tif historic_map_3857_ct.tif

yes this is expected as the BETA2007 grid only goes up to longitude 15.75, and your raster is around 18.5 degree. More generally the DHDN datum of EPSG:4314 is documented to be only valid for former West Germany, up to longitude 13.84. So PROJ auto guessing methods will not work very weel for such data that is outside the normal area of use. Your solution of manually specifying the transformation that works for you is the best one I can recommend.