GDAL python教程基础篇（8）osr 模块简介与用法

在学习osr模块之前我们先来了解一下大地水准面和地球椭球体以及地理坐标系与投影坐标系
大地水准面与椭球体
地球是一个表面很复杂的球体， 人们以假想的平均静止的海水面形成的“大地体”为参照， 推求出近似的椭球体，并通过了理论和实践证明， 该椭球体以地球短轴为轴的椭圆而旋转的椭球面， 而且这个椭球面也可以用数学公式来进行表达。
为此，人们从数学角度出发，选择了一个形状和大小与大地球体均极为相似的旋转椭球面来描述地球表层，这也被称为地球椭球面。具体的定义这里就不叙述了

我国使用的地球椭球体参数列表：

a) 海福特椭球(1910) 我国52年以前基准椭球
a = 6378388m
b = 6356911.9461279m
alpha = 0.33670033670

b) 克拉索夫斯基椭球(1940 Krassovsky) 　北京54坐标系基准椭球
a = 6378245m
b = 6356863.018773m
alpha = 0.33523298692

c) 1975年I.U.G.G推荐椭球(国际大地测量协会1975) 　西安80坐标系基准椭球
a = 6378140m
b = 6356755.2881575m
alpha = 0.0033528131778

d) WGS-84椭球(GPS全球定位系统椭球、17届国际大地测量协会) WGS-84 GPS 基准椭球
a = 6378137m
b = 6356752.3142451m
alpha = 0.00335281006247

sr.SpatialReference 和 osr.CoordinateTransformation 类提供了用来描绘坐标系统（投影和基准面）以及坐标系统相互转换的服务。 这些服务在OpenGIS坐标转换文档里有模型， 并且有对应的WKT描述。
一些OpenGIS坐标系统和服务的背景资料可以在COM的简单要素（Simple Features for COM）和空间参考系统抽象模型文档 (可以从http://opengis.org 获取) 中找到。 GeoTiff投影变换列表也可以辅助地用来理解WKT中的投影。 EPSG测地学网页也是一个有用的资源。

1.定义地理坐标系

地理坐标系被封装进了osr.SpatialReference类。 有几种办法来初始化osr.SpatialReference对象以形成一个合法的坐标系统。 主要有两种坐标系统：一种是地理坐标(用经纬度表示)； 另一种是投影坐标(如 UTM ，用米等单位量度来定位)。

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

OGRSpatialReference oSRS;
oSRS.SetGeogCS( "My geographic coordinate system",
  "WGS_1984",
  "My WGS84 Spheroid",
  SRS_WGS84_SEMIMAJOR, SRS_WGS84_INVFLATTENING,
  "Greenwich", 0.0,
  "degree", SRS_UA_DEGREE_CONV );

OGRSpatialReference 支持一些标准的坐标系统。 比如"NAD27"， “NAD83”， “WGS72” 和 “WGS84”。 要建造它们只要用一个函数SetWellKnownGeogCS()。
oSRS.SetWellKnownGeogCS("WGS84")
如果EPSG数据库存在的话，所有EPSG中的地理坐标系都可以用GCS编码来表示。
oSRS.SetWellKnownGeogCS("EPSG:4326")

在各种坐标系统表达方式中，WKT是个纽带， 通过它，可以互换各种表达方式。 一个OGRSpatialReference可以用一个WKT来初始化，或者转换出WKT表达。

>>> from osgeo import osr
>>> sr = osr.SpatialReference()
>>> sr.SetWellKnownGeogCS('WGS84')
>>> wkt = sr.ExportToWkt()

2.定义投影坐标系

一个投影坐标系统（如UTM等）是基于一个地理坐标系统进行定义的， 它便于在线性单位和角度位置之间进行转换。 接下来的代码定义了一个在WGS84基准面下的 地理坐标系统为基础的UTM17带投影坐标系统。

>>> from osgeo import osr
>>> sr = osr.SpatialReference()
>>> sr.SetProjCS( 'UTM 17 (WGS84) in northern hemisphere.' )
>>> sr.SetWellKnownGeogCS( 'WGS84' )
>>> sr.SetUTM( 17, True )
>>> wkt = sr.ExportToWkt()

调用SetProjCS()设置一个定义用户名字的投影坐标系统并确定系统被投影过。 SetWellKnownGeogCS()分配一个地理坐标系统， SetUTM()设置投影转换参数细节。 上面的步骤直到这时才能创建一个合法的定义。 但是这个对象需要时将会重新自动编排内在表达，以保持合法性。

投影表达坐标系的方式包括：OpenGIS的WKT（Well-Know Text：知名文本）； PROJ.4 的表达方式；EPSG权威编码方式等等

1.3查询坐标系统

osr.SpatialReference建立起来，可以查询多种坐标系统信息。 用IsProjected() 和IsGeographic() 可以看出他是投影坐标还是地理坐标。

>>> from osgeo import osr
>>> sr = osr.SpatialReference()
>>> sr.SetProjCS( 'UTM 17 (WGS84) in northern hemisphere.' )
>>> sr.SetWellKnownGeogCS( 'WGS84' )
>>> sr.SetUTM( 17, True )
>>> sr.IsGeographic()
0

GetSemiMajor(), GetSemiMinor() 和GetInvFlattening()可以获取椭球体信息。 GetAttrValue() 可以用来获取PROJCS, GEOGCS, DATUM, SPHEROID, 和PROJECTION 的名字表达字符串。 GetProjParm()可以用来获取投影参数。 GetLinearUnits()可以用来获取线性单位，并转换成米。
下面的代码来自于osr的测试程序。 示范了GetProjParm()、GetAttrValue() 的用法。 已经定义的投影参数（如SRS_PP_CENTRAL_MERIDIAN） 应该在GetProjParm()获取投影参数的时候使用。 ogrspatialreference.cpp中设置不同投影的方法的代码 可以用来查找哪个投影用哪个参数。

>>> srs = osr.SpatialReference()
>>> from osgeo import gdal
>>> srs.ImportFromUSGS(8, 0,
>>>   (0.0, 0.0,
>>>   gdal.DecToPackedDMS(47.0), gdal.DecToPackedDMS(62.0),
>>>   gdal.DecToPackedDMS(45.0), gdal.DecToPackedDMS(54.5),
>>>   0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
>>>   15)
>>> srs.GetProjParm(osr.SRS_PP_STANDARD_PARALLEL_1)