1. 坐标转化
坐标转化总是一个令人头疼的问题。不知道你有没有写过坐标转换的程序代码,如果写过就会有体会。那一长串的公式和成群出现的大括号,真是太可怕了!只有写过这种东西,你才会对我们生活的这个家园肃然起敬。才会发现原来地球居然如此复杂而显得如此伟大,人类的数学和地理学家居然是如此地会折腾而显得如此伟大。敬畏会让你从此以后不敢乱丢垃圾和随地吐痰~~呵呵,话说大了。反正我是不敢碰那些符号和公式了。不过既然踏入GIS的门槛,就踏进了坐标的深渊,不得不面对它。算了,还是找一个人家已经做出来的轮子去拼自己的车吧,至于轮子有几根钢线,什么材质,怎么才能搞的比较圆,我就不费神去想了。
有一种叫PROJ的牌子的轮子不知道听说过没有。反正开源界现在最好的有关坐标转换的轮子差不多就是它了。对于C/C++人来说,很幸运,因为这个库可以直接调用。如果是在linux平台下的java人,也很幸运,因为可以用jni来间接调用。如果你是在windows平台下的java人,可能是不幸的。我和我老大折腾了两天都没有编译成功,查了无数的帖子,也没有发现有解决方案(如果你成功了,一定要告诉我哦~~~)。对于业余的python人来说,前两天我还以为我们是不幸的,但是今天我觉得上帝还是偏爱python的,让我找到了这个,还有这个。呜啦~~~
对于pyproj,我也是“不太熟悉的”:-),那我们也就看看osr里面坐标转换是怎么玩的吧。 先翻译官方教程,我翻译不对的地方请回头看英文原版
1.1. 官方的教程
1.1.1. 简介
OGRSpatialReference和OGRCoordinateTransformation类提供了用来描绘坐标系统(投影和基准面)以及坐标系统相互之间转换的服务。这些服务在OpenGIS坐标转换说明文档里有模型,并且有对应的WKT描述。 一些OpenGIS坐标系统和服务的背景资料可以在COM的简单要素(Simple Features for COM)和空间参考系统抽象模型文档(可以从opengis.org获取)中找到。GeoTiff投影变换列表也可以辅助地用来理解WKT中的投影公式。EPSG测地学网页也是一个有用的资源。
1.1.2. 定义一个地理坐标系
地理坐标系被封装进了OGRSpatialReference类。有几种办法来初始化OGRSpatialReference对象以形成一个合法的坐标系统。有两种主要的坐标系统。一种是地理坐标(用经纬度表示)。另一种是投影坐标(如UTM,用米等单位量度来定位)。
一个地理坐标包含基准面(它包含了由长半轴描述的椭球体和反扁率),本初子午线(一般是格林威治),和一个角度单位(一般是度)。下面的代码就初始化了一个地理坐标系。它提供了围绕用户定义名字的地理坐标系的所有信息。
1 OGRSpatialReference oSRS;
2 oSRS.SetGeogCS( "My geographic coordinate system",
3 "WGS_1984",
4 "My WGS84 Spheroid",
5 SRS_WGS84_SEMIMAJOR, SRS_WGS84_INVFLATTENING,
6 "Greenwich", 0.0,
7 "degree", SRS_UA_DEGREE_CONV );
在这些值中,"My geographic coordinate system", "My WGS84 Spheroid", "Greenwich" 和 "degree" 不是关键字。但是被用于显示给用户看。但是基准面名"WGS_1984" 经常被用于定义基准面。而且哪些值可以用哪些不行都是有规矩的。不能乱写。
OGRSpatialReference 已经支持一些标准的坐标系统。比如"NAD27", "NAD83", "WGS72" and "WGS84"。要建造它们只要用一个函数SetWellKnownGeogCS().
oSRS.SetWellKnownGeogCS( "WGS84" );
如果EPSG数据库存在的话,所有EPSG中的地理坐标系都可以用GCS编码来表示。
oSRS.SetWellKnownGeogCS( "EPSG:4326" );
在各种坐标系统表达方式中,WKT是个纽带,通过它,各种表达方式可以互换。一个OGRSpatialReference可以用一个wkt来初始化,或者转换出wkt表达。
1 char *pszWKT = NULL;
2 oSRS.SetWellKnownGeogCS( "WGS84" );
3 oSRS.exportToWkt( &pszWKT );
4 printf( "%s\n", pszWKT );
比如这样一个wkt:
1 GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,
2 AUTHORITY["EPSG",7030]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG",6326]],
3 PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],UNIT["DMSH",0.0174532925199433,
4 AUTHORITY["EPSG",9108]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG",
5 4326]]
或者让它更好看一点吧:
1 GEOGCS["WGS 84",
2 DATUM["WGS_1984",
3 SPHEROID["WGS 84",6378137,298.257223563,
4 AUTHORITY["EPSG",7030]],
5 TOWGS84[0,0,0,0,0,0,0],
6 AUTHORITY["EPSG",6326]],
7 PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],
8 UNIT["DMSH",0.0174532925199433,AUTHORITY["EPSG",9108]],
9 AXIS["Lat",NORTH],
10 AXIS["Long",EAST],
11 AUTHORITY["EPSG",4326]]
OGRSpatialReference::importFromWkt()方法可以用来把wkt坐标系统设置到OGRSpatialReference中。