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.
_______________________________________________
gdal-dev mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
