GDAL coordinate transformation using C/C++ API (OGRCoordinateTransformation)

Question

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.

Answer 1

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.

-s_srs <srs_def>

Source spatial reference to use. The coordinate systems that can be passed are anything supported by the OGRSpatialReference.SetFromUserInput() call, which includes EPSG, PCS, and GCSes (i.e. EPSG:4296), PROJ.4 declarations (as above), or the name of a .prj file containing well known text. If no srs_def is given, the srs_def of the source tiles is used (if there is any). The srs_def will be propagated to created tiles (if possible) and to the optional shape file(s)

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.

gdaltransform -s_srs epsg:32635 -t_srs epsg:3301
Enter X Y [Z [T]] values separated by space, and press Return.
499980 6500040
674129.197241796 6503984.25180551 0

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.

projinfo -s epsg:32635 -t epsg:3301
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Inverse of UTM zone 35N + Inverse of EST97 to WGS 84 (1) + Estonian National Grid, 1 m, Estonia - onshore and offshore.

PROJ string:
+proj=pipeline
  +step +inv +proj=utm +zone=35 +ellps=WGS84
  +step +proj=lcc +lat_0=57.5175539305556 +lon_0=24 +lat_1=59.3333333333333
        +lat_2=58 +x_0=500000 +y_0=6375000 +ellps=GRS80
  +step +proj=axisswap +order=2,1