Python地理数据处理十一：空间参照系统(SRS)

1. 空间参照

空间参考系统由三个部分组成：坐标系统，基准面和投影，所有这些都会影响一组坐标在地球上的对应位置。使用基准来表示地球的曲率，而投影则将坐标从三维地球仪转换为二维地图。

地球椭球具有多种模型，这些模型被称为基准，每一个空间参考系统是基于其中之一。广泛使用的一种全球基准，世界大地测量系统于1984年进行了修订，该基准简称为WGS84。是一种用于具有全球覆盖范围的数据，包括全球定位系统（GPS）。大多数基准被设计为在更局部的区域中模拟地球的曲率。为一个区域设计的基准将无法在其他区域使用，如，欧洲不应该使用1983年的北美基准面（NAD83）。
中断地图时从三维转换到二维的一种方法，下面列出一些投影图像：

2. OSR空间参考

2.1 空间参考对象

osgeo包中包含一个叫作OSR（OGR Spatial Reference）空间参考的模块，用于处理SRS。我们可以利用它来进行查看数据，并且转换数据。
可以使用GetSpatialRef()函数来得到SRS，如果图层不具有SRS信息，则返回None。

import os
from osgeo import ogr
import ospybook as pb
from ospybook.vectorplotter import VectorPlotterdata_dir = r'E:\gis with python\osgeopy data'from osgeo import osr# 查看SRS
ds = ogr.Open(os.path.join(data_dir, 'US', 'states_48.shp'))
srs = ds.GetLayer().GetSpatialRef()# 知文本Well Known Text (WKT)
print(srs)

WKT:

这不是一个投影的SRS，因为它不具有PROTCS条目，只有一个GEOGCS对象。
用WKT表示的NAD83基准面UTM12带的空间参考系统：

WKT不是表示SRS的唯一字符串，可以用更加简洁的PROJ.4字符串：

>>> print(srs.ExportToProj4())
+proj=longlat +datum=NAD83 +no_defs

XML格式：

>>> print(srs.ExportToXML())

显示结果：

查询SRS是地理投影还是已投影，并不需要输出什么东西，可以使用 IsGeographic 和 IsProjected 函数提供这些信息。可以使用GetAttrValue()函数获得SRS对应的关键字，比如，

1.获得投影的名称：

>>> print(utm_sr.GetAttrValue('PROJCS'))
NAD83 / UTM zone 12N

2. 查看UTM SRS的AUTHORITY:

>>> # Get the authority.
print(utm_sr.GetAttrValue('AUTHORITY'))
EPSG

3.返回AUTHORITY的第二个数值：

>>> print(utm_sr.GetAttrValue('AUTHORITY', 1))
26912

26912是SRS的最后一个，而不是第一个，因为其最少嵌套，所有是函数返回的第一个。

4.获取感兴趣的数值：

>>> print(utm_sr.GetAuthorityCode('DATUM'))
6269

用GetAuthorityCode（）函数传递感兴趣的SRID的键值。

5.使用GetProjParm（）的PARAMETER键值获取对应的数值：

>>> print(utm_sr.GetProjParm(osr.SRS_PP_FALSE_EASTING))
500000.0

2.2 创建控件参考对象

因为不能总是可以从一个图层中或几何对象中获取合适的空间参考对象，所有我们需要根据要求自定义。这里有两种比较简短的表示方式：
1.使用标准的 EPSG（European Petroleum Survey Group，欧洲石油调查组织）代码；
2.使用 PROJ.4 字符串。

在前面的数据中，NAD83 UTM 12N 对应的EPSG代码是：26912，可以引入osr模块，将其传递给ImportFromEPSG函数，并创建一个空的Spatial Reference对象：

from osgeo import osr# 从EPSG代码创建一个UTM SRS
sr = osr.SpatialReference()
print(sr.ImportFromEPSG(26912))
print(sr.GetAttrValue('PROJCS'))

0
NAD83 / UTM zone 12N

ImportFromEPSG（）函数返回的0表示SRS被成功导入。

1. 使用PROJ.4字符串：

# 从project .4字符串创建一个UTM SRS.
sr = osr.SpatialReference()
sr.ImportFromProj4('''+proj=utm +zone=12 +ellps=GRS80+towgs84=0,0,0,0,0,0,0 +units=m +no_defs ''')
print(sr.GetAttrValue('PROJCS'))

0
unknown

2.从WKT字符串来创建空间参考对象，而不需要使用导入程序函数：

# 从WKT字符串创建一个未投影的SRS
wkt = '''GEOGCS["GCS_North_American_1983",DATUM["North_American_Datum_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]'''
sr = osr.SpatialReference(wkt)
print(sr)

GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433],AXIS["Longitude",EAST],AXIS["Latitude",NORTH]]

3.自定义构建一个空间参考对象，从UTM分支创建阿尔伯斯（Albers）等面积投影（美国地质调查局用于美国48个州）：

# 使用参数创建Albers sr
# 创建一个空的SRS
# 设置名称
# 指定一个基准
# 最后为Albers投影提供所需参数
sr = osr.SpatialReference()
sr.SetProjCS('USGS Albers')
sr.SetWellKnownGeogCS('NAD83')
sr.SetACEA(29.5, 45.5, 23, -96, 0, 0)
#sr.Fixup()  # 为缺少的参数添加默认值，并重新排序项目，但这样做似乎有问题
sr.Validate()
print(sr)

PROJCS["USGS Albers",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",23],PARAMETER["longitude_of_center",-96],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

SetACEA（stdp1，stdp2，clat，clong，fe，fn）介绍：
参数依次是：标准平行线1，标准平行线2，中央纬线，中央经线，东移假定值和北移假定值。

注意： sr.Validate() 返回0，表示一切正常；若返回5，表示SRS已损坏。

2.3 分配SRS

最好将SRS信息附加到数据集上，可以为图层和单个几何对象分配SRS。由于数据中每个图层都可以具有不同的SRS，所以不能向数据源分配SRS。
CreateLayer（） 函数的一个参数是空间参考对象，但osr不能确定数据自身使用的SRS，所以这个参数的默认值为None。如有SRS，需要在 新建图层时就提供，因为不能向现有的图层添加SRS。

# 确保之前的数据目录中存在输出文件夹
out_fn = os.path.join(data_dir, 'output', 'testdata.shp')# 创建一个使用UTM SRS的空shapefile文件
# 如果运行，它将使用包含SRS信息的.prj文件创建shapefile文件
sr = osr.SpatialReference()
sr.ImportFromEPSG(26912)
ds = ogr.GetDriverByName('ESRI Shapefile').CreateDataSource(out_fn)
# 新建图层时，SRS是第二个参数sr
lyr = ds.CreateLayer('counties', sr, ogr.wkbPolygon)

向图层分配SRS时，并不会将所有的数据转换到该坐标系中，它所做的就是提供相关信息。
若使用的单个几何对象，而不是图层，可能需要为几何对象分配SRS，使用 AssignSpatialReference（） 方法来操作：

geom.AssignSpatialReference(sr)

不强制使用该SRS，只是提供关于SRS的信息。

2.4 重投影

当我得到一个新的数据集时，但是它与我通常使用的SRS不同，就必须进行重投影。
重投影的两种方法：假设这个几何对象已经有一个SRS分配给它，而另一个没有。
1.首先，获取覆盖全球的复合多边形：

# ne_110m_land_1p.shp包含一个表示世界陆地的复合多边形
data_dir = r'E:\Googlechrome'world = pb.get_shp_geom(os.path.join(data_dir, 'global', 'ne_110m_land_1p.shp'))
vp = VectorPlotter(True)
vp.plot(world)
vp.draw()# 为埃菲尔铁塔创建一个点
tower = ogr.Geometry(wkt='POINT (2.294694 48.858093)')  # 经纬度坐标
tower.AssignSpatialReference(osr.SpatialReference(osr.SRS_WKT_WGS84))# 将世界多边形重新投影到Web墨卡托
# 这是一个错误
web_mercator_sr = osr.SpatialReference()
web_mercator_sr.ImportFromEPSG(3857)
world.TransformTo(web_mercator_sr)

使用内置函数跳过某些点（如tower、北极、南极等），通过导入 gdal 模块并使用其中的 SetConfigOption 函数来修复错误：

# 设置一个配置变量，然后再试一次投影
from osgeo import gdal
gdal.SetConfigOption('OGR_ENABLE_PARTIAL_REPROJECTION', 'TRUE')
web_mercator_sr = osr.SpatialReference()
web_mercator_sr.ImportFromEPSG(3857)
world.TransformTo(web_mercator_sr)
tower.TransformTo(web_mercator_sr)
print(tower)
vp.clear()
vp.plot(world)
vp.draw()

使用Web墨卡托投影绘制显示的全球陆地图：

假设全球几何对象没有空间参考，使用 CoordinateTransformation 方法从Web墨卡托投影转换为Gall-Peters投影，使用Transform函数，需要一个CoordinateTransformation对象：

from osgeo import gdal
gdal.SetConfigOption('OGR_ENABLE_PARTIAL_REPROJECTION', 'TRUE')
web_mercator_sr = osr.SpatialReference()
web_mercator_sr.ImportFromEPSG(3857)
world.TransformTo(web_mercator_sr)# 在Web Mercator和Gall-Peters之间创建一个坐标转换
peters_sr = osr.SpatialReference()
peters_sr.ImportFromProj4("""+proj=cea +lon_0=0 +x_0=0 +y_0=0+lat_ts=45 +ellps=WGS84 +datum=WGS84+units=m +no_defs""")
ct = osr.CoordinateTransformation(web_mercator_sr, peters_sr)
world.Transform(ct)
vp.clear()
vp.plot(world, 'c')
vp.draw()

使用Gall-Peters投影绘制显示的全球陆地：

2.5 更改基准

比如将NAD27基准数据转换为需要的NAD83基准，可以使用TransformTo和Transform函数在基准面之间进行转换。
基准之间进行转换的数学公式并不总是存在，所以很多GIS软件使用称为网格移位文件的数据文件来帮助转换。这些文件包含将坐标从一个基准转换为另一个基准所需要的信息。osr模块可以在系统中找到这些文件，但必须保证原始SRS和目标SRS都包含基准信息。
如果没有用于基准转换的适当栅格移位文件，则可以为源和目标SRS设置 towgs84 参数。这些参数描述了从特定基准到WGS84基准的近似转换。如果不知道所需的参数，可以从这里查看：EPSG Geodetic Parameter Dataset

# 创建一个未投影的NAD27 SRS，并添加数据移位信息
sr = osr.SpatialReference()
sr.SetWellKnownGeogCS('NAD27')
sr.SetTOWGS84(-8, 160, 176)
print(sr)

2.6 重投影整个图层

创建新图层后，需要循环遍历原始图层中的每个要素，获取几何对象，并对其将进行转换，然后将包含已转换几何对象的要素插入到新图层中。可能需要保留图层的所有属性字段，因此需要在创建新字段时从原始图层复制字段定义。

投影一个点图层：

# 重新投影shapefile的脚本
from osgeo import ogr, osr# 创建一个输出SRS
sr = osr.SpatialReference()
sr.ImportFromProj4('''+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23+lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80+datum=NAD83 +units=m +no_defs''')# 不要忘记在这里更改目录
ds = ogr.Open(r'E:\gis with python\osgeopy data\US', 1)# 获取输入层
in_lyr = ds.GetLayer('us_volcanos')# 创建空的输出层
out_lyr = ds.CreateLayer('us_volcanos_aea', sr,ogr.wkbPoint)# 复制字段定义
out_lyr.CreateFields(in_lyr.schema)# 循环遍历输入层中的属性字段
out_feat = ogr.Feature(out_lyr.GetLayerDefn())
for in_feat in in_lyr:# 转换几何对象，并复制属性信息geom = in_feat.geometry().Clone()geom.TransformTo(sr)out_feat.SetGeometry(geom)# 复制属性for i in range(in_feat.GetFieldCount()):out_feat.SetField(i, in_feat.GetField(i))# 插入属性out_lyr.CreateFeature(out_feat)

3. 使用pyproj空间参考

PROJ.4制图投影库适用于在SRS之间转换数据的C语言库。可以被各种开源库使用，包括OSR。pyproj模块提供了一个用于PROJ.4的python包装器。
与OSR处理几何对象不同，它是处理坐标值列表，可以作为python列表、元组、数组、NumPy数组或标量提供。如果包含坐标集的文本文件，则pyproj模块中包含的函数将是它们转换为其他坐标系统的理想方式。

3.1 在不同SRS中转换坐标

两种方法：使用Proj类在地理坐标和投影坐标之间进行转换；使用模块级Transform函数在两个SRS之间进行转换。
将埃菲尔铁塔坐标从经纬度转换为UTM Zone 31N：
先用PROJ.4字符串为UTM SRS初始化Proj对象，然后使用它来转换坐标：

# 需要下载pyproj模块
import pyproj
utm_proj = pyproj.Proj('+proj=utm +zone=31 +ellps=WGS84')
x, y = utm_proj(2.294694, 48.858093)
print(x, y)# 从投影坐标转换到地理坐标，需要将inverse设置为Tue
x1, y1 = utm_proj(x, y, inverse=True)   # 传入UTM坐标
print(x1, y1)  # 获取了原始的经纬度坐标值，但是存在四舍五入

448265.9146797117 5411920.651461348
2.2946940000000002 48.858093000000004

使用transform函数，进行基准转换：

# 将UTM WGS84坐标转换为UTM NAD27
# 三种方法初始化Proj对象
from functools import partial
from pyproj import Proj, transformwgs84 = pyproj.Proj('+proj=utm +zone=18 +datum=WGS84')
nad27 = pyproj.Proj('+proj=utm +zone=18 +datum=NAD27')
x, y = pyproj.transform(wgs84, nad27, 580744.32, 4504695.26) # WGS84基准
print(x, y) # NAD27基准

580710.2007512685 4504483.45735296

两个基准在东西方向上相差30米左右，南北方向上相差200米左右

由于更新，出现了以下的修改：

3.2 计算大圆距离

地球上两个点之间的距离叫做大圆距离。使用 pyproj 获得两组经纬度坐标之间的距离，和它们之间大圆线的起始和结束方位。
首先用要使用的椭圆体实例化一个Geod类对象，然后传递十进制的开始和结束坐标给其inv函数，来获得向前方向、向后方向和距离：

import pyproj# 设置洛杉矶和柏林的经纬度坐标
la_lat, la_lon = 34.0500, -118.2500
berlin_lat, berlin_lon = 52.5167, 13.3833# 创建一个WGS84 Geod
geod = pyproj.Geod(ellps='WGS84')# 找到洛杉矶和柏林之间的方位和距离
forward, back, dist = geod.inv(la_lon, la_lat, berlin_lon, berlin_lat)
print('forward: {}\nback: {}\ndist: {}'.format(forward, back, dist))# 如果从柏林出发，找到你的最终坐标,然后从后方出发,这些坐标应该与LA匹配
x, y, bearing = geod.fwd(berlin_lon, berlin_lat, back, dist)
print('{}, {}\n{}'.format(x, y, bearing))# 找到一份从洛杉矶到柏林的等距坐标表
# npts函数来生成用于绘制两地之间的大圆路径
coords = geod.npts(la_lon, la_lat, berlin_lon, berlin_lat, 100)# 只打印前3个
for i in range(3):print(coords[i])

forward: 27.23284045673669
back: -38.49148498662069
dist: 9331934.878166698
-118.24999999999996, 34.05
27.23284045673668
(-117.78803196383676, 34.789725145004155)
(-117.31774994946879, 35.52757560403803)
(-116.83878951054419, 36.2634683783333)