[gdal-dev] gdal2tiles and EPSG:27700

[Javier Jimenez Shaw]

Hi Even,

Using "gdalwarp -ct" was a good solution. It worked fine. Thanks!
However, when I tested it with another GeoTIFF, this case is from Japan,
gdalwarp is producing something strange.

The image is in EPSG:2455 "JGD2000 / Japan Plane Rectangular CS XIII" ,
that in EPSG is defined as northing-easting.
It is located here:
$ gdalinfo input.tif | grep Center
Center      (  -88109.342, -129144.325) (143d10'20.17"E, 42d49'56.66"N)

I generate the pipeline this way
$ projinfo -s EPSG:2455 -t EPSG:3857 -o PROJ --single-line -q
+proj=pipeline +step +proj=axisswap +order=2,1 +step +inv +proj=tmerc
+lat_0=44 +lon_0=144.25 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +step
+proj=webmerc +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84

But gdalwarp is producing strange outputs with that pipeline.
Letting gdalwarp read the CRS, it works fine:
$ gdalwarp -t_srs EPSG:3857 input.tif to3857.tif
$ gdalinfo to3857.tif | grep Center
Center      (15937864.003, 5286496.999) (143d10'20.17"E, 42d49'56.66"N)

Setting both CRSs works fine as well:
$ gdalwarp -s_srs EPSG:2455 -t_srs EPSG:3857 input.tif from2455to3857.tif
$ gdalinfo from2455to3857.tif | grep Center
Center      (15937864.003, 5286496.999) (143d10'20.17"E, 42d49'56.66"N)

However using "-ct", swaps the output:
$ gdalwarp -ct '+proj=pipeline +step +proj=axisswap +order=2,1 +step +inv
+proj=tmerc +lat_0=44 +lon_0=144.25 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80
+step +proj=webmerc +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84' input.tif
pipeline.tif
$ gdal_edit.py -a_srs EPSG:3857 pipeline.tif
$ gdalinfo pipeline.tif | grep Center
Center      ( 5286496.998,15937864.005) ( 47d29'21.88"E, 80d36'13.79"N)

Am I doing something wrong, or is there any problem in "gdalwarp"?

Thanks.
.___ ._ ..._ .. . ._.  .___ .. __ . _. . __..  ... .... ._ .__
Entre dos pensamientos racionales
hay infinitos pensamientos irracionales.



On Wed, 3 Nov 2021 at 15:03, Even Rouault <even.rouault at spatialys.com>
wrote:

>
> 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.
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20211115/b702220b/attachment.html>