Search code examples
c++gdalcoordinate-systemsmingw-w64

Reprojecting coordinates using GDAL library

I'm trying to reproject a lat/lon coordinate in the WGS84 system to an UTM one in the SIRGAS 2000 coord system.

Using GDAL, I started with changing a known utm coordinate to the lat/lon counterpart in the same coordinate system (29 N), just to check that I wrote the right code (I'm omitting error-checking here):

OGRSpatialReference monUtm;
monUtm.SetWellKnownGeogCS("WGS84");
monUtm.SetUTM(29, true);

OGRSpatialReference monGeo;
monGeo.SetWellKnownGeogCS("WGS84");

OGRCoordinateTransformation* coordTrans = OGRCreateCoordinateTransformation(&monUtm, &monGeo);

double x = 621921.3413490148;
double y = 4794536.070196861;

int reprojected = coordTrans->Transform(1, &x, &y);
// If OK, print the coords.

delete coordTrans;
coordTrans = OGRCreateCoordinateTransformation(&monGeo, &monUtm);
reprojected = coordTrans->Transform(1, &x, &y);

// If OK, Print the coords.
delete coordTrans;

The coordinates 621921.3413490148, 4794536.070196861 correspond to the Moncelos region in northern Galicia. The forth-and-back transformation seems to work right: the lat/lon coordiantes are the correct ones and when projecting back to UTM, I get the same as the originals:

UTM:          621921.34135 , 4794536.0702
Lat/lon:      43.293779579 , -7.4970160261
Back to UTM:  621921.34135 , 4794536.0702

Now, reprojecting from WGS84 lat/long to SIRGAS 2000 UTM:

// Rodovia dos Tamoios, Brazil:
//  - UTM          -> 23 S
//  - WGS 84       -> EPSG:32723
//  - SIRGAS 2000  -> EPSG:31983

OGRSpatialReference wgs;
wgs.SetWellKnownGeogCS("WGS84");

OGRSpatialReference sirgas;
sirgas.importFromEPSG(31983);

coordTrans = OGRCreateCoordinateTransformation(&wgs, &sirgas);

double x = -23.57014667;
double y = -45.49159617;

reprojected = coordTrans->Transform(1, &x, &y);
// If OK, print results
delete coordTrans;

coordTrans = OGRCreateCoordinateTransformation(&sirgas, &wgs);
reprojected = coordTrans->Transform(1, &x, &y);
// If OK, print results.

this doesn't give the same results:

WGS84 Lat/Lon input:      -23.57014667  ,  -45.49159617
SIRGAS 2000 UTM output:   2173024.0216  ,  4734004.2131
Back to WGS84 Lat/Lon:    -23.570633824 ,  -45.491627598

As you can see, the original WGS84 lat/lon and the back-to_WGS84 lat/lon coords aren't exactly the same, unlike the first test case. Also, the UTM x-coord has 7 digits (I thought it was limited to 6 (?) ).

In Google Maps, we can see that there's a difference of 27 meters between the two points (original point is represented by a circle. My "back-reprojected" point is represented by a dagger).

Finally, the question: am I doing the reprojection right? If so, why is there a 27 meter difference between reprojections in the second test case?

Solution

  • The problem is that you need to swap your axis order to use Cartesian X/Y space or Lon/Lat, and not "Lat/Lon" order.

    Setting this should work.

    double x = -45.49159617;  // Lon
double y = -23.57014667;  // Lat

    The difference you saw from your round-trip conversion was from projecting outside the bounds for the UTM zone due to a swapped axis order.