（二）GDAL开发&与esri的wkt

• 问题来源

最近在使用gdal来处理一些遥感数据，其中我的数据是一个波段的tif高程影像数据。然后在程序中需要判断当然影像数据的投影是否符合要求。我们都知道，如何使用某一个坐标系需要一个约定，而这个公认的，为大家所知wkt算是一个之一。关于wkt，其为空间参考文字描述。我们最为知晓的应该就是编号id为4326和3957这个了。前者为WGS84的空间参考，而3857为一种墨卡托投影参考。而我这里使用的ArcGIS Desktop将我们CGCS2000的椭球坐标定义为WGS84椭球坐标（在精度要求不是非常严格的情况下，可视这两种椭球坐标相当）。

• 过程处理

于是在ArcGIS Desktop中将上面的数据按照如下图的定义。注意其中的WKT为3857。

接下来就是使用GDAL将投影信息给读取出来了。代码实现非常简单。

const char* projection_wkt = dataset->GetProjectionRef();
if(projection_wkt == NULL){
       
        return false;
}

获取到的projextion_wkt的字符信息为：

PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],
PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0],
EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m [email protected] +wktext  +no_defs"]]

从这里可以看得出，ArcGIS Desktop声明了WGS_1984_Web_Mercator_Auxiliary_Sphere椭球，地理坐标椭球为GCS_WGS_1984，椭球信息SPHEROID长半轴为：6378137.0,短半轴为：298.257223563,起始子午线为格林尼治经线。单位为米。

好了我们现在使用GDAL的OGRSpatialReference类来创建一个参考信息。

	char    *EPSGID_WKT = NULL;
	OGRSpatialReference   myespEPSGIDosr;
	myespEPSGIDosr.importFromEPSG((int)3857);

	myespEPSGIDosr.exportToWkt(&EPSGID_WKT);

来看一下EPSGID_WKT带给我们什么样的信息。

PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m [email protected] +wktext  +no_defs"],AUTHORITY["EPSG","3857"]]

显然，同样这里的WKT是3857，但输出的字符信息包含却是完全不同。主要区别是椭球的标记一个是GCS_WGS_1984，另外一个是WGS_1984。

即便GDAL提供从esri投影转换相关函数。比如说morphToESRI、和morphFromESRI但是，最后的wkt描述字符串还是不完全一致。我这里为了判定两个坐标参考是否相等，只能通过其他逻辑来判定了。

---

                             更多内容，请微信扫下面二维码关注公众号，或者加入arcpy开发qq学习群：487352121