GDAL坐标转换

1、引言

最近研究了一下GIS、测绘学的坐标转换的问题，感受大部分资料专业性太强，上来就是一通专业性论述；但感受对于相关从业者来讲，其实没必要了解那么多背景知识的；就经过GDAL这个工具，来简单总结下坐标转换相关的问题。
GDAL坐标转换其实也是调用proj4来实现，可是proj4有个特别麻烦的地方，就是坐标系描述的部分特别繁复，须要对专业知识有必定的了解。使用GDAL则相对简单不少。

2、地理坐标系

地理坐标系就是常说的经纬度坐标系，好比用GPS直接获取的坐标就是在地理坐标系下获取的。一个真实坐标不管怎么变换，必定会有地理坐标系做为基准，也必定有能够转换出来的经纬度坐标。以国内的状况来讲，经常使用到的地理坐标系有WGS84，beijing54，xian80，CGCS2000这四种。
GDAL能够像proj4那样自定义坐标系，也能够仅经过字符串定义一些经常使用的坐标系，但本文认为最方便的仍是经过EPSG数据库定义的代码来定义一个地理坐标系统；毕竟不少时候须要兼容的地理坐标系不少，所有一个个自定义坐标系太麻烦。

OGRSpatialReference spatialReference;
//spatialReference.importFromEPSG(4326);//WGS84
//spatialReference.importFromEPSG(4214);//BeiJing54
spatialReference.importFromEPSG(4610);//XIAN80
//spatialReference.importFromEPSG(4490);//CGCS2000

char *pszWKT = nullptr;
spatialReference.exportToPrettyWkt(&pszWKT);
cout << pszWKT << endl;
CPLFree(pszWKT);
pszWKT = nullptr;

如上演示了GDAL定义地理坐标系并输出坐标系的信息。四种不一样的坐标系只须要改变相应的代码编号，用起来十分方便。最终打印输出了的xian80坐标系信息：

在这里必定要注意，使用这种方式定义地理坐标系必定要经过配置GDAL_DATA路径，不然控制台会输出错误信息：

CPLSetConfigOption("GDAL_DATA", "gdaldata");

“gdaldata”表示一个路径（这里用的是相对路径，固然也能够设置成绝对路径），是GDAL编译完成后会生成的一个目录，里面记录了各类坐标系的参数文件。例如进入这个文件夹，用Excel打开pcs.csv这个文件，就能够看到各类坐标系的参数了。

除了这种方法，也能够在环境变量中设置GDAL_DATA变量来实现。

3、投影平面坐标系

经纬度坐标是曲面上的坐标，曲面上的坐标投影到平面，不一样的投影方式就会产生不一样的平面坐标；即便是同一种投影方式，不一样的参数获得的平面坐标也会不一样。也就是说，地理坐标系下的经纬度坐标与投影坐标系下的平面坐标，是一对多的关系而不是一对一的关系。以国内的状况来讲，常常用到的投影有横轴墨卡托投影，高斯-克吕格投影和UTM投影。
在Global Mapper中三种投影设置方式以下：



能够看出，高斯-克吕格投影和UTM投影其实都是横轴墨卡托投影，二者都是经过设置带号（Zone）来实现设置横轴墨卡托投影的具体参数（Parameters）。在GDAL里面，高斯-克吕格投影就是经过设置横轴墨卡托投影来实现的。以下演示了一个xian80坐标系，3度带带号38的横轴墨卡托投影。

OGRSpatialReference spatialReference;
spatialReference.importFromEPSG(4610);          //XIAN80
spatialReference.SetTM(0, 114, 1.0, 38500000, 0);

设置横轴墨卡托投影的函数SetTM()有六个参数：

参数	设置
dfCenterLat	通常为0
dfCenterLong	中央经线，决定了是哪一投影带
dfScale	通常为1.0，UTM设置为0.9996
dfFalseEasting	通常为500000，高斯-克吕格前面加上带号
dfFalseNorthing	通常为0

用以前方法获得坐标系信息并输出，信息以下：

4、坐标转换

定义好坐标系以后，就能够进行坐标转换了。以下为xian80地理坐标系下某点（113.6，38.8）用高斯-克吕格投影到带平面坐标系：

OGRSpatialReference* pLonLat = spatialReference.CloneGeogCS();
OGRCoordinateTransformation* LonLat2XY = OGRCreateCoordinateTransformation(pLonLat, &spatialReference);
if (!LonLat2XY)
{       
    return 1;
}
    
double x = 113.6;
double y = 38.8;
printf("经纬度坐标：%.9lf\t%.9lf\n", x, y);
if (!LonLat2XY->Transform(1, &x, &y))
{
    return 1;
}
printf("平面坐标：%.9lf\t%.9lf\n", x, y);

OGRCoordinateTransformation::DestroyCT(LonLat2XY);
LonLat2XY = nullptr;

这里经过复制以前定义的高斯-克吕格投影平面坐标系获得相同的地理坐标系（固然也能够自定义新坐标系），而后使用OGRCoordinateTransformation::Transform()方法来进行坐标转换。最后的输出结果为：

经过Global Mapper的坐标转换工具来验证结果是否正确：


比较能够发现二者转换的结果基本一致。除此以外，将平面坐标逆投影到地理坐标也是能够的，只须要在OGRCreateCoordinateTransformation（）的时候颠倒下顺序便可。

5、其余

1.GDAL默认不编译proj.4，使用的时候须要添加proj.4的支持。
2.同一地理坐标系的投影转换是严密的，但不一样地理坐标系之间须要先转换到地心立体坐标系，而后经过七参数转换。
3.能够根据坐标值选择正确的分带，使用这个分带的上下几个分带进行投影问题也不是很大。可是分带差距太大可能有问题，由于投影公式只能在必定范围内有效；即不一样的软件对的时候结果比较一致，错的时候结果可能不一样。
以上三个问题有时间再作进一步的研究和总结。

6、参考文献

1.GDAL源码剖析（十一）之OGR投影说明
2.墨卡托投影、高斯-克吕格投影、UTM投影及我国分带方法
3.GDAL库学习笔记(五):坐标系之间的转化
4.GIS坐标转换库Proj.4的使用
5.GDAL影像投影转换