0

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.

4
  • 1
    What means "Results are the same"?
    – user30184
    May 2, 2022 at 13:33
  • @user30184 It means that the results of the C/C++ GDAL and Python GDAL examples are the same.
    – arapEST
    May 3, 2022 at 5:59
  • Can you show a reference to confirm that coordinates should be the same in EPSG:32635 and EPSG:3301? Or is that what you mean?
    – user30184
    May 3, 2022 at 6:17
  • @user30184 Thanks for discussion. Reference that I use is the output of gdal tools. For example the given gdal_retile.py -ps 512 512 -overlap 32 -of Gtiff -s_srs 'EPSG:3301' -targetDir tiled_multi T35VNE_20211102T093049_B8A.jp2 is converting the input dataset's EPSG:32635 coordinates to EPSG:3301 and they are not changed, due to the fact that input projection is in a UTM zone that is similar to EPSG:3301's (covering Estonia). Therefore I think gdal tools' output should be the template here.
    – arapEST
    May 3, 2022 at 13:16

1 Answer 1

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
1
  • Very much appreciated, how can I buy you a beer/coffee? In the midst of all the details I stupidly missed the definition of -s_srs. This is such a relief. Thanks for extra explanation and gdaltransform example.
    – arapEST
    May 3, 2022 at 17:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.