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.roua...@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. > >
_______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev