[gdal-dev] gdal2tiles and EPSG:27700

[Even Rouault]

Le 03/11/2021 à 14:49, Javier Jimenez Shaw a écrit :
> Hi
>
> I have a GeoTIFF in EPSG:27700, "OSGB36 / British National Grid"
> For consistency reasons, I want to run gdal2tiles with an equivalent 
> transformation than the one I will use later to transform points using 
> a PROJ pipeline. For consistency reasons I will use that pipeline 
> along the time to transform the points, so an update in EPSG/PROJ does 
> not change the results (I generate the gdal2tiles at project creation 
> time, but run the transformation on the marker points along the time).
>
> Getting the PROJ pipeline is easy with something like
> projinfo -s EPSG:27700 -t EPSG:3857 --spatial-test intersects -o PROJ 
> --single-line -q
> (the option --spatial-test intersects is very important. Otherwise you 
> get a ballpark. I am not using grids. I can use the area of use of the 
> GeoTIFF as well) It produces
> +proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 
> +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +step +proj=push 
> +v_3 +step +proj=cart +ellps=airy +step +proj=helmert +x=446.448 
> +y=-125.157 +z=542.06 +rx=0.15 +ry=0.247 +rz=0.842 +s=-20.489 
> +convention=position_vector +step +inv +proj=cart +ellps=WGS84 +step 
> +proj=pop +v_3 +step +proj=webmerc +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 
> +ellps=WGS84
>
> However this pipeline does not work in gdal2tiles with the "-s" 
> option, complaining with these errors, and placing the data in the 
> wrong place (not very far from Null Island)
> ERROR 1: PROJ: proj_crs_get_coordinate_system: Object is not a SingleCRS
> ERROR 1: PROJ: proj_as_wkt: PROJBasedOperation can only be exported to 
> WKT2
> ERROR 1: PROJ: proj_crs_get_coordinate_system: Object is not a SingleCRS
> ERROR 1: PROJ: proj_as_wkt: PROJBasedOperation can only be exported to 
> WKT2
This is expected. A coordinate operation is not a source CRS.
>
> Running projinfo on EPSG:27700 gives me the ballpark one
> projinfo EPSG:27700 -o PROJ --single-line
> PROJ.4 string:
> +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 
> +y_0=-100000 +ellps=airy +units=m +no_defs +type=crs
>
> The difference using this pipeline (without Helmert parameters) is 
> about 100 m near London. Perfectly noticeable ;)
>
> (I think it does not include the Helmert transformation, because the 
> area of use of that transformation is a bit smaller than the 
> geographic crs, EPSG:4277)
>
> Of course I could create the command manually, but I would like to 
> have a generic solution for multiple projects.
>
> What can I do?

Enhance gdal2tiles to be able to specify coordinate transformations: 
https://github.com/OSGeo/gdal/issues/3998

You could also possibly use gdalwarp -ct to directly generate the input 
dataset in EPSG:3857 with the pipeline of your choice

Or in that instance just use -s "+proj=tmerc +lat_0=49 +lon_0=-2 
+k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs 
+towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489"

Even

-- 
http://www.spatialys.com
My software is free, but my time generally not.