I'm in a puzzle of figuring out how to create correct coordinate transformations using GDAL C/C++ API.
I have referenced code from tutorial – https://gdal.org/tutorials/osr_api_tut.html#coordinate-transformation
Input coordinate references would be based on Sentinel-2 dataset over Estonia. Hence source SRS – EPSG:32635 (WGS 84 / UTM 35N). Destination SRS – EPSG:3301 (Estonian Coordinate System of 1997).
When running the tutorial example as (args – {source EPSG} {dest EPSG} {lon} {lat}):
gdal-srs 32635 3301 499980 6500040
Result would be
(499980.000000,6500040.000000) -> (6503984.251806,674129.197242)
I tried the same with the following python tutorial example – https://pcjericks.github.io/py-gdalogr-cookbook/projection.html#reproject-a-geometry Results are the same. Also tried specifying options with SetAreaOfInterest()
– still no change.
To make sure the installation is not messed up (GDAL 3.0.4 official apt package for Ubuntu 20.04), the following gdal_retile command produces correct coordinates for output tiles (Same tool is used for reference too)
gdal_retile.py -ps 512 512 -overlap 32 -of Gtiff -s_srs 'EPSG:3301' -targetDir tiled_multi T35VNE_20211102T093049_B8A.jp2
I have debugged gdal_retile.py source code but I have not understood what it does differently. In practice this transformation between these EPSGs should not actually change the input coordinates, but the functionality should expect all the projections to be supported that GDAL/PROJ supports.
Best Answer
The gdal_retile python script https://gdal.org/programs/gdal_retile.html does not warp. It creates a new tiled file structure in the same coordinate system than the original. The
-s_srs
option is used for describing the coordinate system of the source data if it is otherwise unknown, or to override the internal coordinate system information of the source data if user knows that it is wrong in the metadata.What you have effectively done is to tell gdal_retile that source image is in EPSG:3301, that is wrong because you seem to know that it is actually in EPSG:32635. Gdal_retile does not warp but the result files get labeled to be in EPSG:3301.
If you want to warp and retile do it in two steps:
gdalwarp -of VRT -s_srs epsg:32635 -t_srs epsg:3301 T35VNE_20211102T093049_B8A.jp2 T35VNE_20211102T093049_B8A.vrt
gdal_retile.py -ps 512 512 -overlap 32 -of Gtiff -targetDir tiled_multi T35VNE_20211102T093049_B8A.vrt
EPSG:32635 and EPSG:3301 are not at all the same and your results with C++ are correct. Notice that in the result from gdaltransform X and Y coordinates are flipped because the Estonian EPSG:3301 is using coordinate axis order Northing-Easting https://epsg.org/crs_3301/Estonian-Coordinate-System-of-1997.html.
Projinfo shows details about the steps that are needed for that conversion. You can see for example that EPSG:326535 is a UTM projection but the Estonian one is Lambert Conic Conformal.