GDAL地理坐标投影坐标转换类OGRCoordinateTransformation的成员函数Transform的参数顺序
OGRCoordinateTransformation类的坐标转换方法Transform
OGRCoordinateTransformation::Transform的函数原型为：
int Transform(size_t nCount, double *x, double *y, double *z = nullptr,int *pabSuccess = nullptr);
可以一次性转换nCount个坐标。X、Y、Z坐标分别存储在x,y,z三个数组中。转换完成后，x,y,z数组中的值都被改变，所以如果要保留之前的值，需要用户自行注意。
请看下面两段示例代码：
(656090.245245055, 2900041.074837993, 0.0)是WGS 84椭球UTM zone49投影带（中央经线111°E）下的坐标，对应的WGS 84椭球上的经纬度坐标为(112.562380,26.211450)。
1. 将投影坐标转为经纬度坐标：
int main()
{OGRSpatialReference oUTM, * poLongLat;OGRCoordinateTransformation* poTransform;oUTM.SetProjCS("UTM 49 / WGS84");oUTM.SetWellKnownGeogCS("WGS84");oUTM.SetUTM(49);poLongLat = oUTM.CloneGeogCS();// WGS84poTransform = OGRCreateCoordinateTransformation(&oUTM, poLongLat);int nPoints = 1;double x = 656090.245245055, y = 2900041.074837993, z=0.0;printf("Input UTM (x,y):(%f,%f)\n", x, y);if (!poTransform->Transform(nPoints, &x, &y, &z))printf("Transformation failed.\n");else{printf("Output (lon,lat):(%f,%f)\n", x, y);}return 0;
}2.将经纬度坐标转为UTM投影坐标
OGRSpatialReference oUTM, * poLongLat;OGRCoordinateTransformation* poTransform;oUTM.SetProjCS("UTM 49 / WGS84");oUTM.SetWellKnownGeogCS("WGS84");oUTM.SetUTM(49);poLongLat = oUTM.CloneGeogCS();// WGS84poTransform = OGRCreateCoordinateTransformation(poLongLat, &oUTM);// 从经纬度转到UTM投影坐标系int nPoints = 1;double x = 656090.245245055, y = 2900041.074837993, z=0.0;double lon = 112.562380, lat = 26.211450;if (!poTransform->Transform(nPoints, &lon, &lat, &z)){printf("Transformation failed.\n");}else{printf("Output UTM (x,y):(%f,%f)\n", lon, lat);}
问题解析
这两段看起来毫无问题，然而，第一段代码的输出是：
Input UTM (x,y):(656090.245245,2900041.074838)
Output (lon,lat):(26.211450,112.562380)
第二段代码的输出是：
ERROR 1: PROJ: utm: Invalid latitude
Transformation failed.
其实在GDAL 3.0以前的版本中，上面两段代码的输出都是正常的；然而从GDAL 3.0开始，输出就会产生上述异常。
为什么会这样呢？这与OGRSpatialReference的变动有关。从GDAL 3.0开始，对于最常用的空间参考系统WGS 84（“EPSG:4326”或“WGS84”），OGRSpatialReference类的对象默认接受经纬度的顺序发生了改变。这导致调用Transform函数输入/输出经纬度坐标时，参数x、y不再按以前的约定（x(East,经度)-y(North,纬度)）对应经纬度，而是变为x(North,纬度)-y(East,经度)；
实际应用时，需要特别注意。所以，在第一段代码中，我们认为变量x,y被转换后应该分别是经度和纬度，所以输出的lon,lat对应反了；而第二段代码中，我们调用Transform的传参顺序是 &lon, &lat，GDAL认为我们传入的经度值竟然是112.52380，自然会报错。
要解决这个问题，需要在坐标转换前加一句，将OGRSpatialReference类的对象的经纬度的顺序设置为x(East,经度)-y(North,纬度)：
poLongLat->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
或者老老实实的按照新的顺序修改程序：
第一段代码：
上略if (!poTransform->Transform(nPoints, &x, &y, &z))printf("Transformation failed.\n");else{printf("Output (lon,lat):(%f,%f)\n", y, x);}
第二段代码：
上略if (!poTransform->Transform(nPoints, &lat, &lon, &z)){printf("Transformation failed.\n");}else{printf("Output UTM (x,y):(%f,%f)\n", lat, lon);}
SetAxisMappingStrategy
参考网页
