osr简介

osr.SpatialReferenceosr.CoordinateTransformation 类提供了用来描绘坐标系统(投影和基准面)以及坐标系统相互转换的服务。 这些服务在OpenGIS坐标转换文档里有模型, 并且有对应的WKT描述。

一些OpenGIS坐标系统和服务的背景资料可以在COM的简单要素(Simple Features for COM)和空间参考系统抽象模型文档 (可以从http://opengis.org 获取) 中找到。 GeoTiff投影变换列表也可以辅助地用来理解WKT中的投影。 EPSG测地学网页也是一个有用的资源。

Python 绑定

C的接口在ogr_srs_api.h中定义, Python绑定在osr.py模块中。 一些C++的方法在C和Python绑定中没有定义。

服务的实现是建立在 PROJ.4 库之上, 该库是最先由USGS的Gerald Evenden 编写的。

Pyton绑定的形式为:

class osr.SpatialReference
    def __init__(self,obj=None):
    def ImportFromWkt( self, wkt ):
    def ExportToWkt(self):
    def ImportFromEPSG(self,code):
    def IsGeographic(self):
    def IsProjected(self):
    def GetAttrValue(self, name, child = 0):
    def SetAttrValue(self, name, value):
    def SetWellKnownGeogCS(self, name):
    def SetProjCS(self, name = "unnamed" ):
    def IsSameGeogCS(self, other):
    def IsSame(self, other):
    def SetLinearUnits(self, units_name, to_meters ):
    def SetUTM(self, zone, is_north = 1):
class CoordinateTransformation:
    def __init__(self,source,target):
    def TransformPoint(self, x, y, z = 0):
    def TransformPoints(self, points):    

定义地理坐标系

地理坐标系被封装进了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 );

在这些值中,"My geographic coordinate system""My WGS84 Spheroid""Greenwich""degree" 不是关键字。 但是被显示给用户看。 但是基准面名“WGS_1984” 经常被用做定义基准面。 而且哪些值可以用哪些不行都是有规矩的,且不能乱写。

OGRSpatialReference 支持一些标准的坐标系统。 比如"NAD27""NAD83""WGS72""WGS84"。 要建造它们只要用一个函数SetWellKnownGeogCS()

oSRS.SetWellKnownGeogCS("WGS84");

如果EPSG数据库存在的话,所有EPSG中的地理坐标系都可以用GCS编码来表示。

oSRS.SetWellKnownGeogCS("EPSG:4326");

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

In [21]:
from osgeo import osr
sr = osr.SpatialReference()
sr.SetWellKnownGeogCS('WGS84')
wkt = sr.ExportToWkt()
wkt
Out[21]:
'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"]]'

比如这样一个WKT,进行了格式排版,以便让它更美观一些。

GEOGCS["WGS 84",
  DATUM["WGS_1984",
    SPHEROID["WGS 84",6378137,298.257223563,
        AUTHORITY["EPSG",7030]],
    TOWGS84[0,0,0,0,0,0,0],
    AUTHORITY["EPSG",6326]],
  PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],
  UNIT["DMSH",0.0174532925199433,AUTHORITY["EPSG",9108]],
  AXIS["Lat",NORTH],
  AXIS["Long",EAST],
  AUTHORITY["EPSG",4326]]

ImportFromWkt()函数可以把 WKT坐标系统设置到OGRSpatialReference中。

定义一个投影系统

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

In [22]:
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()
wkt
Out[22]:
'PROJCS["UTM 17 (WGS84) in northern hemisphere.",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["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-81],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]'
In [23]:
%pprint
wkt
Pretty printing has been turned ON
Out[23]:
'PROJCS["UTM 17 (WGS84) in northern hemisphere.",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["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-81],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]'
In [24]:
from pprint import pprint
pprint(wkt)
('PROJCS["UTM 17 (WGS84) in northern hemisphere.",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["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-81],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]')

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

注意:定义对象时的步骤和顺序!!

上面的定义会定义一个和下面一样的WKT。注意UTM117被扩展成横轴莫卡托定义UTM带的参数。

美化输出 :

PROJCS["UTM 17 (WGS84) in northern hemisphere.",
  GEOGCS["WGS 84",
    DATUM["WGS_1984",
      SPHEROID["WGS 84",6378137,298.257223563,
        AUTHORITY["EPSG",7030]],
      TOWGS84[0,0,0,0,0,0,0],
      AUTHORITY["EPSG",6326]],
    PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],
    UNIT["DMSH",0.0174532925199433,AUTHORITY["EPSG",9108]],
    AXIS["Lat",NORTH],
    AXIS["Long",EAST],
    AUTHORITY["EPSG",4326]],
  PROJECTION["Transverse_Mercator"],
  PARAMETER["latitude_of_origin",0],
  PARAMETER["central_meridian",-81],
  PARAMETER["scale_factor",0.9996],
  PARAMETER["false_easting",500000],
  PARAMETER["false_northing",0]]

不同的投影会有不同的函数来设定参数:

SetTM()  
(Transverse Mercator)

SetLCC()  
(Lambert Conformal Conic)

SetMercator()  
(Merator)

投影的表示方法

投影表达坐标系的方式包括:OpenGIS的WKT(Well-Know Text:知名文本); PROJ.4 的表达方式;EPSG权威编码方式等等。

常用投影定义的表达方式

可以进行地图投影定义的常用格式:

  • NAD27/NAD83/WGS84/WGS72: 这些常见的地理坐标系统可以通过名字来直接使用。
  • EPSG:n坐标系统(投影或者地理坐标)可以通过EPSG码来选择。 例如EPSG 27700是英国国家网格。 一系列的EPSG坐标系统可以在GDAL数据文件gcs.csv和pcs.csv中找到。
  • PROJ.4定义: 一个PROJ.4定义字符串可以用坐标系统定义。 例如“+proj=utm +zone=11 +datum=WGS84”。 注意在命令行中要保持PROJ.4字符串在一起作为一个单独的参数(一般用双引号引起来)。
  • OpenGIS 知名文本: OpenGIS标准定义了一个文本格式来描述坐标系统作为简单要素规范的一个部分。 这个格式是GDAL中使用的坐标系统的内部工作格式。 包含WKT坐标系统描述的文件名可以被用来作为坐标系统参数, 或者坐标系统元素本身也可以被用来作为命令行参数(虽然WKT中释放所有引号是很有争议的)。
  • ESRI知名文本: ESRI 在他们的ArcGIS产品(ArcGIS .prj文件)中使用了一种经过精简OGC WKT的格式, 而且这个格式被用在一个和WKT相似的风格的文件中。 但是文件名要被加以ESRI::前缀。比如“ESRI::NAD 1927 StatePlane Wyoming West FIPS 4904.prj”。
In [ ]:
 

查询坐标系统

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

In [25]:
from osgeo import osr
sr = osr.SpatialReference()
sr.SetProjCS( 'UTM 17 (WGS84) in northern hemisphere.' )
sr.SetWellKnownGeogCS( 'WGS84' )
sr.SetUTM( 17, True )
sr.IsGeographic()
Out[25]:
0
In [26]:
sr.ExportToProj4()
Out[26]:
'+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs '
In [27]:
sr.AutoIdentifyEPSG()
Out[27]:
0
In [28]:
print(sr.ExportToPrettyWkt())
PROJCS["UTM 17 (WGS84) in northern hemisphere.",
    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["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-81],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1],
    AUTHORITY["EPSG","32617"]]
In [ ]:
 
In [29]:
sr.IsProjected()
Out[29]:
1

GetSemiMajor(), GetSemiMinor()GetInvFlattening()可以获取椭球体信息。 GetAttrValue() 可以用来获取PROJCS, GEOGCS, DATUM, SPHEROID, 和PROJECTION 的名字表达字符串。 GetProjParm()可以用来获取投影参数。 GetLinearUnits()可以用来获取线性单位,并转换成米。

下面的代码来自于osr的测试程序。 示范了GetProjParm()GetAttrValue() 的用法。 已经定义的投影参数(如SRS_PP_CENTRAL_MERIDIAN) 应该在GetProjParm()获取投影参数的时候使用。 ogrspatialreference.cpp中设置不同投影的方法的代码 可以用来查找哪个投影用哪个参数。

In [30]:
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)
Out[30]:
47.0
In [31]:
srs.GetProjParm(osr.SRS_PP_STANDARD_PARALLEL_2)
Out[31]:
62.0
In [32]:
srs.GetProjParm(osr.SRS_PP_LATITUDE_OF_CENTER)
Out[32]:
54.5
In [33]:
srs.GetProjParm(osr.SRS_PP_LONGITUDE_OF_CENTER)
Out[33]:
45.0
In [34]:
srs.GetProjParm(osr.SRS_PP_FALSE_EASTING)
Out[34]:
0.0
In [35]:
srs.GetProjParm(osr.SRS_PP_FALSE_NORTHING)
Out[35]:
0.0
In [ ]:
 

栅格数据的空间参考

空间信息包括两个部分: 一个是坐标系统,另一个是栅格坐标和地理坐标之间的转换。

GDAL文档的解释

数据集的坐标系统由OpenGIS WKT字符串定义,它包含了:

  1. 一个总体的坐标系名
  2. 一个地理图形坐标系统名
  3. 一个基准面定义
  4. 一个椭球体的名字。长半轴(semi-major axis)和反扁率(inverse flattening)
  5. 本初子午线(prime meridian)名和其与格林威治子午线的偏移值
  6. 投影的方法类型(如横轴莫卡托)
  7. 投影的参数列表(如中央经线等)
  8. 一个单位的名称和其到米和弧度单位的转换参数
  9. 轴线的名称和顺序
  10. 在预定义的权威坐标系中的编码(如EPSG)

更多信息请参考OpenGIS WKT坐标系统定义,以及osr教程文档和OGRSpatialReference类的描述文档。

在GDAL中,返回坐标系统的函数是GDALDataset::GetProjectionRef()。 它返回的坐标系统描述了地理参考坐标,暗含着仿射地理参考转换, 地理参考转换是由GDALDataset::GetGeoTransform()来返回。

GCPs地理参考坐标描述的坐标系统是由GDALDataset::GetGCPProjection()返回的。

注意:返回的坐标系统字符串“”表示未知的地理参考坐标系统。

坐标系统

坐标系统在GDAL中涉及的坐标系统概念。 根据GDAL数据模型文档 [1]的解释, GDAL数据集类中包括的坐标系统解释是这样的:

仿射地理变换

把数据在坐标系统中定位,需要了解数据的坐标。 在这一点的处理上,栅格数据与矢量数据采用了完全不同的数据模型和处理方式。

矢量数据直接把坐标信息存储到数据本身。 每一个点 [2],都具有其相应的地理坐标。

栅格数据则只存储了左上角像元的坐标, 其他各个像元的坐标,则依靠像元大小, 以及在X方向与Y方向和原点(左上角像元)的偏移来计算。

使用仿射变换进行空间定位

In [36]:
from osgeo import gdal
dataset = gdal.Open("/gdata/geotiff_file.tif")
dataset.GetGCPs()
dataset.GetGeoTransform()
Out[36]:
(1868454.913, 30.0, 0.0, 5353126.266, 0.0, -30.0)

只有通过这几个数据, 我们才可以知道自己手中这张图该如何铺展到绘制平面上。 上面说过了,栅格坐标到地理坐标之间的转换可以通过两种方式来进行: 一种是仿射坐标转换,一种是GCPs。 一般来说,仿射坐标转换比较常见。 仿射坐标转换就是通过几个值来进行栅格框架到地理框架的映射。 上面的例子中出现了6个值。下面这是六个值可以用的: 

Xgeo = GT(0) + XpixelGT(1) + YlineGT(2)

Ygeo = GT(3) + XpixelGT(4) + YlineGT(5)

通过公式来解释,使数组顺序对应起来, 理解起来也很简单。下面是几个参数的意思:

adfGeoTransform[0] /* top left x */
adfGeoTransform[1] /* w-e pixel resolution */
adfGeoTransform[2] /* rotation, 0 if image is "north up" */
adfGeoTransform[3] /* top left y */
adfGeoTransform[4] /* rotation, 0 if image is "north up" */
adfGeoTransform[5] /* n-s pixel resolution */

看起来和ESRI为GeoTiff定义的tfw文件格式很像。

有了上面的公式, 可以把栅格的每个点代入公式求出其在地球上的实际位置 (当然,求出的数值的单位是在坐标系统中定义的那个,而非经纬度坐标, 要求经纬度坐标,还要通过一步转换)。 当然我们也可以 只求出整个图像外框上下左右四个边界就可以定下整个图像的位置, 因为图像都是矩形的。 图像中某点的地理位置可以通过外框四点线性计算。

现在看看我们这张图的边界地理范围:

In [37]:
dataset.GetGeoTransform()
dataset.RasterXSize
dataset.RasterYSize
Out[37]:
900

上边框和左边框都不算,则结果是4928000和590000。 求右下角的坐标,代入公式:

右边框 = Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2) 
       = 590000.0 + 950*20.0+ 700*0.0 = 609000.0
下边框 = Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5) 
       = 4928000.0 + 950*0.0+ 700*(-20.0) = 4914000.0

有了坐标,就可以把图像和地理坐标联系起来了。 有了地理坐标,就可以把图像在照屏幕上预先定好的坐标系画出来了。

需要注意的是,对于GeoTiff数据源来说, 坐标系统和仿射变换参数有两种存储方式, 一种是直接存储在.tif文件内部, 这种存储方式是按照tif内部的的键值方式来存储的。 官方说明pdf可以在这里下载。 也可以把图像存储到tif文件中, 而把坐标系统和仿射参数什么的单独提取出来, 创建两个文件,一个是prj文件,存放WKT坐标系字符串,另一个是tfw文件,存放仿射转换参数。 还有一些软件在处理时将二者都保存起来:一方面将投影的信息存储到数据文件本身, 另一方面使用单独的文本文件保存投影信息。 这两种方式都有好处:即是将投影信息存储到数据文本本身进行数据管理, 在进行数据的拷贝到转移时不会因为落下文件而造成投影信息的丢失; 在进行处理时使用单独文本文件保存投影信息会比较方便, 只需要使用简单的文本处理程序,而不必使用复杂的地理信息软件。

GDAL存储栅格数据坐标

本节来看一下在GDAL中如何实现存储栅格数据坐标的。

GDAL数据集有两种模式是描述栅格位置(用点/线坐标系)以及地理参考坐标系之间的关系: 一种是使用仿射转换, 另一种则是GCPs(多控制点定位方式)

一个仿射变换,是一个用于栅格数据的公式,公式如下:

这些操作中一个以上的可以立即应用;这允许您执行复杂的转换(如旋转)。

地面控制点(GCPs),是将栅格影像中的点与实际的地面坐标联系起来的点。 在GDAL中,也可以保存在遥感影像内部。

注意: GDAL 并不使用 GCPs 来进行空间定位,这个是留给应用时使用的, 并且一般要使用较为复杂的数学函数来执行这个变换。

让我们看看坐标系统究竟是什么样子。

In [38]:
sr = dataset.GetProjectionRef()
sr
Out[38]:
'PROJCS["Albers_Beijing54",GEOGCS["GCS_Krasovsky_1940",DATUM["Not_specified_based_on_Krassowsky_1940_ellipsoid",SPHEROID["Krassowsky 1940",6378245,298.2999999999998,AUTHORITY["EPSG","7024"]],AUTHORITY["EPSG","6024"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",47],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",105],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]'

GDAL 的坐标系统的标示方法是用OpenGIS 的 WKT 字符串表示, 当然也混合了 EPSG 权威编码。

使用GCPs表示空间参考

使用GCPs进行空间定位

对于GCPs, 需要配准一张栅格图(把一个没有空间信息的图片变成有空间信息的栅格数据集), 只要在一个已知的坐标系统中定位几个控制点,然后输入这几个点的地理位置。 软件通过读取这几个点在图像中的位置,以及输入的这几个点对应的地理位置, 建立一个拟合方程,用方程来描述其他点在这个坐标系统中的地理位置。 当然,在不同投影不同坐标系不同软件中拟合方程处理可能不同。 所以拟合出来的坐标可能有一些差异。

数据集可以由一系列控制点来定义空间参考坐标系。 所有的GCPs共用一个地理参考坐标系(由GDALDataset::GetGCPProjection()返回, 每个GCP(由GDAL_GCP定义)包含下面内容:

typedef struct
{
char *pszId;
char *pszInfo;
double dfGCPPixel;
double dfGCPLine;
double dfGCPX;
double dfGCPY;
double dfGCPZ;
} GDAL_GCP;

pszid是本GCP在数据集中一系列GCP点中唯一的标示字符串, 它常常是数字,但不一定总是。 有可能在GCP状态中包含机器可分析信息,虽然现在还不行。

(dfGCPPixel,dfGCPLine)位置是栅格中的GCP位置, (dfGCPX,dfGCPY,dfGCPZ)位置是联合的地理参考位置(Z通常是0)。

GDAL数据模型没有实现由GCPs产生坐标系的变化的机制,而是把它留给实际应用。 但是1到5阶多项式是通常使用的方法。

通常一个数据集会包含仿射地理变换和GCPS中的一个或两个都没有。 两个都有的情况是很少见,而且无法用权威坐标系定义。


In [39]:
import helper; helper.info()
页面更新时间: 2019-02-28 18:41:39
操作系统/OS: Linux-4.9.0-8-amd64-x86_64-with-debian-9.8
Python: 3.5.3
In [ ]: