GDAL库学习笔记(五):坐标系之间的转化
1. 坐标转化 坐标转化总是一个令人头疼的问题。不知道你有没有写过坐标转换的程序代码,如果写过就会有体会。那一长串的公式和成群出现的大括号,真是太可怕了!只有写过这种东西,你才会对我们生活的这个家园肃然起敬。才会发现原来地球居然如此复杂而显...
- 作者:lilin来源:啄木鸟CPUG站|2007年02月08日
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中。