<div dir="ltr"><div>Hi Even,</div><div><br></div><div>Using "gdalwarp -ct" was a good solution. It worked fine. Thanks!<br></div><div>However, when I tested it with another GeoTIFF, this case is from Japan, gdalwarp is producing something strange.</div><div><br></div><div>The image is in EPSG:2455 "JGD2000 / Japan Plane Rectangular CS XIII" , that in EPSG is defined as northing-easting.</div><div>It is located here:<br><span style="font-family:monospace">$ gdalinfo input.tif | grep Center<br>Center ( -88109.342, -129144.325) (143d10'20.17"E, 42d49'56.66"N)</span></div><div><br></div><div>I generate the pipeline this way<br></div><div><span style="font-family:monospace">$ projinfo -s EPSG:2455 -t EPSG:3857 -o PROJ --single-line -q<br>+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</span></div><div><br></div><div>But gdalwarp is producing strange outputs with that pipeline.</div><div>Letting gdalwarp read the CRS, it works fine:</div><div><span style="font-family:monospace">$ gdalwarp -t_srs EPSG:3857 input.tif to3857.tif</span></div><div><span style="font-family:monospace">$ gdalinfo to3857.tif | grep Center<br>Center (15937864.003, 5286496.999) (143d10'20.17"E, 42d49'56.66"N)</span></div><div><br></div><div>Setting both CRSs works fine as well:<br><span style="font-family:monospace">$ gdalwarp -s_srs EPSG:2455 -t_srs EPSG:3857 input.tif from2455to3857.tif<br>$ gdalinfo from2455to3857.tif | grep Center<br>Center (15937864.003, 5286496.999) (143d10'20.17"E, 42d49'56.66"N)</span></div><div><br></div><div>However using "-ct", swaps the output:</div><div><span style="font-family:monospace">$ 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<br>$ gdal_edit.py -a_srs EPSG:3857 pipeline.tif<br>$ gdalinfo pipeline.tif | grep Center<br>Center ( 5286496.998,15937864.005) ( 47d29'21.88"E, 80d36'13.79"N)</span></div><div><br></div><div>Am I doing something wrong, or is there any problem in "gdalwarp"?</div><div><br></div><div>Thanks.<br></div><div><div><div dir="ltr" class="gmail_signature" data-smartmail="gmail_signature">.___ ._ ..._ .. . ._. .___ .. __ . _. . __.. ... .... ._ .__<br>Entre dos pensamientos racionales <br>hay infinitos pensamientos irracionales.<br><br></div></div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, 3 Nov 2021 at 15:03, Even Rouault <<a href="mailto:even.rouault@spatialys.com">even.rouault@spatialys.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><br>
Le 03/11/2021 à 14:49, Javier Jimenez Shaw a écrit :<br>
> Hi<br>
><br>
> I have a GeoTIFF in EPSG:27700, "OSGB36 / British National Grid"<br>
> For consistency reasons, I want to run gdal2tiles with an equivalent <br>
> transformation than the one I will use later to transform points using <br>
> a PROJ pipeline. For consistency reasons I will use that pipeline <br>
> along the time to transform the points, so an update in EPSG/PROJ does <br>
> not change the results (I generate the gdal2tiles at project creation <br>
> time, but run the transformation on the marker points along the time).<br>
><br>
> Getting the PROJ pipeline is easy with something like<br>
> projinfo -s EPSG:27700 -t EPSG:3857 --spatial-test intersects -o PROJ <br>
> --single-line -q<br>
> (the option --spatial-test intersects is very important. Otherwise you <br>
> get a ballpark. I am not using grids. I can use the area of use of the <br>
> GeoTIFF as well) It produces<br>
> +proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 <br>
> +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +step +proj=push <br>
> +v_3 +step +proj=cart +ellps=airy +step +proj=helmert +x=446.448 <br>
> +y=-125.157 +z=542.06 +rx=0.15 +ry=0.247 +rz=0.842 +s=-20.489 <br>
> +convention=position_vector +step +inv +proj=cart +ellps=WGS84 +step <br>
> +proj=pop +v_3 +step +proj=webmerc +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 <br>
> +ellps=WGS84<br>
><br>
> However this pipeline does not work in gdal2tiles with the "-s" <br>
> option, complaining with these errors, and placing the data in the <br>
> wrong place (not very far from Null Island)<br>
> ERROR 1: PROJ: proj_crs_get_coordinate_system: Object is not a SingleCRS<br>
> ERROR 1: PROJ: proj_as_wkt: PROJBasedOperation can only be exported to <br>
> WKT2<br>
> ERROR 1: PROJ: proj_crs_get_coordinate_system: Object is not a SingleCRS<br>
> ERROR 1: PROJ: proj_as_wkt: PROJBasedOperation can only be exported to <br>
> WKT2<br>
This is expected. A coordinate operation is not a source CRS.<br>
><br>
> Running projinfo on EPSG:27700 gives me the ballpark one<br>
> projinfo EPSG:27700 -o PROJ --single-line<br>
> PROJ.4 string:<br>
> +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 <br>
> +y_0=-100000 +ellps=airy +units=m +no_defs +type=crs<br>
><br>
> The difference using this pipeline (without Helmert parameters) is <br>
> about 100 m near London. Perfectly noticeable ;)<br>
><br>
> (I think it does not include the Helmert transformation, because the <br>
> area of use of that transformation is a bit smaller than the <br>
> geographic crs, EPSG:4277)<br>
><br>
> Of course I could create the command manually, but I would like to <br>
> have a generic solution for multiple projects.<br>
><br>
> What can I do?<br>
<br>
Enhance gdal2tiles to be able to specify coordinate transformations: <br>
<a href="https://github.com/OSGeo/gdal/issues/3998" rel="noreferrer" target="_blank">https://github.com/OSGeo/gdal/issues/3998</a><br>
<br>
You could also possibly use gdalwarp -ct to directly generate the input <br>
dataset in EPSG:3857 with the pipeline of your choice<br>
<br>
Or in that instance just use -s "+proj=tmerc +lat_0=49 +lon_0=-2 <br>
+k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs <br>
+towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489"<br>
<br>
Even<br>
<br>
-- <br>
<a href="http://www.spatialys.com" rel="noreferrer" target="_blank">http://www.spatialys.com</a><br>
My software is free, but my time generally not.<br>
<br>
</blockquote></div>