|
|
|
|

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,用米等单位量度来定位)。

一个地理坐标包含基准面(它包含了由长半轴描述的椭球体和反扁率),本初子午线(一般是格林威治),和一个角度单位(一般是度)。下面的代码就初始化了一个地理坐标系。它提供了围绕用户定义名字的地理坐标系的所有信息。

Toggle line numbers
						
								   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表达。

Toggle line numbers
						
								   1 
								char
								*
								pszWKT
								=
								NULL
								;
						
						
								   2 
								oSRS
								.
								SetWellKnownGeogCS
								(
								"WGS84"
								)
								;
						
						
								   3 
								oSRS
								.
								exportToWkt
								(
								&
								pszWKT
								)
								;
						
						
								   4 
								printf
								(
								"%s\n"
								,
								pszWKT
								)
								;
								
								
						
				

比如这样一个wkt:

Toggle line numbers
						
								   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
								]
								]
								
								
						
				

或者让它更好看一点吧:

Toggle line numbers
						
								   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中。

上一篇:GDAL库学习笔记(四):关于空间参考

下一篇:GDAL库学习笔记(六): 把dem地形转化成vrm…