Wrong projection from Geographic to Projected CS in gdal 3.x

Question

I have the following trivial code, which transforms coordinates from Geographic CS (EPSG 4230, ED50) to Projected CS (EPSG 23021 ED50 / UTM zone 31N)

#include <ogr_spatialref.h>
#include <ogr_core.h>
#include <ogr_srs_api.h>

int main(int argc, char *argv[])
{
  double coord[3] = {11.0, 60.0, 0.0};
  int res;
  
  OGRCoordinateTransformationH transform = NULL;
  OGRSpatialReference oSRS, iSRS;
  oSRS.importFromEPSG(23031);
  iSRS.SetWellKnownGeogCS( "EPSG:4230" );
  
  transform = OCTNewCoordinateTransformation(OGRSpatialReference::ToHandle(&iSRS), OGRSpatialReference::ToHandle(&oSRS));
  OCTTransform(transform, 1, coord+0, coord+1, coord+2);

  printf("%lf %lf %lf\n", coord[0], coord[1]);
}

Linking this against gday version below 3.x prints the correct result:

945515.849694 6678568.489194

However, if I use any gdal 3.x, I get the following values, which are wrong:

7952467.873008 2186283.396945

I can't see what I'm doing wrong. Is this a bug, or is there a mistake in my code?

Answer 1

Read the "CRS and axis order" section of the OSR API tutorial.

Basically, until GDAL 3.0 the OGRSpatialReference class ignored the axis order (lon, lat v.s lat, lon) specified by the CRS and as of 3.0 GDAL honours the specified axis order.

To go back to the pre-3.0 behaviour, just set OGRSpatialReference::SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER)

I'm not sure of the syntax for C/C++ so below is a python example:

from osgeo import osr

coord = (11.0, 60.0, 0.0)
oSRS = osr.SpatialReference()
iSRS = osr.SpatialReference()
iSRS.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
iSRS.ImportFromEPSG(4230)
oSRS.ImportFromEPSG(23031)

transform = osr.CoordinateTransformation(iSRS, oSRS)
p = transform.TransformPoint(coord)
print(p)

Output:

(945515.8496942197, 6678568.489194222, 0.0)