本文介绍了如何在Python中使用osr模块定义和操作地理坐标系（经纬度）和投影坐标系（米制），包括WGS84、UTM投影及EPSG编码的应用。通过实例展示了如何创建、转换和验证不同的空间参考系，并探讨了栅格和矢量数据的坐标获取方法。
OSR模块和空间参考系
定义地理坐标系
地理坐标系封装在osr.SpatialReference类中。
主要有两种坐标系：经纬度的地理坐标系和米的投影坐标系。
# -*- coding: utf-8 -*-
# 学习时间： 2022/11/26 17:05
__author__ = 'He XK'
from osgeo import ogr, osr, gdal
# 定义SR空间参考对象
sr = osr.SpatialReference()
# 给SR对象设定地理坐标系为WGS84
sr.SetWellKnownGeogCS('WGS84')
# 导出成wkt格式
wkt = sr.ExportToWkt()
print(wkt)
输出结果：
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"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
注意：
以上声明坐标系统采用的参数是'WGS84'。也可以使用标准的EPSG声明方式，比如用'EPSG:4326'进行等价的代替。
可以通过sr.ImportFromWkt(wkt)方法将自定义的wkt格式写入SR中。
在哪里查看SetWellKnownGeoCS的const参数？
定义投影系统
定义在WGS84基准面下的北半球UTM17带。
sr.SetProjCS(char const * name="unnamed")设置一个定义用户名字的投影坐标系统并确定系统被投影过，应该是参数是用户随意输入，只表示用于定义坐标系的名称；
sr.SetWellKnownGeogCS(...)设定一个坐标系统，EPSG中已有的各种参数；
sr.SetUTM(int zone, int north=1)设置投影转换参数细节，包含两个参数，zone表示投影带号，north表示是否为北半球，默认是True。
sr.SetProjCS('UTM 17 (WGS84) in northern hemisphere.')
sr.SetWellKnownGeogCS('WGS84')
sr.SetUTM(17 , True)
wkt = sr.ExportToWkt()
pprint(wkt) # pretty print
关于地理坐标系和投影坐标系
sr = osr.SpatialReference()
sr.SetWellKnownGeogCS('EPSG:4214')
wkt = sr.ExportToWkt()
pprint(wkt)
print(sr.IsGeographic()) # 1
sr2 = osr.SpatialReference()
sr2.SetProjCS('4214 ZONE 18 > 4573 metre')
sr2.SetWellKnownGeogCS('EPSG:4214')
sr2.SetUTM(18, 1)
wkt2 = sr2.ExportToWkt()
pprint(wkt2)
print(sr2.IsProjected()) # 1
地理坐标系：Degree
投影坐标系：Metre
投影坐标系是以地理坐标系为基准面。
以EPSG：4214为例，表示Beijing 1954地理坐标系。
SR只设定了地理基准面，所以sr.IsGeographic()的判断为True；
SR2依据此地理基准面进行投影，同时设定18带，经过sr2.IsProjected()的判断为True
可以 漂亮的输出所有空间参考系信息。
print(sr2.ExportToPrettyWkt())
可以输出空间投影坐标系的长半轴SemiMajor、短半轴SemiMinor和扁率InvFlattening
print(sr2.GetSemiMajor()) # 6378245.0
print(sr2.GetSemiMinor()) # 6356863.018773047
print(sr2.GetInvFlattening()) # 298.3
可以获取投影坐标系中的各种参数。
名称可以用ExportToPrettyWkt输出后获取查看。
print(sr2.GetAttrValue('PROJCS')) # 4214 ZONE 18 > 4573 metre
print(sr2.GetAttrValue('GEOGCS')) # Beijing 1954
print(sr2.GetAttrValue('DATUM')) # Beijing_1954
print(sr2.GetAttrValue('PRIMEM')) # Greenwich
print(sr2.GetAttrValue('PROJECTION')) # Transverse_Mercator
栅格数据的空间参考
数据集的坐标系统由OpenGIS WKT字符串定义，它包含了：
-
一个总体的坐标系名；
-
一个地理图形坐标系统名；
-
一个基准面定义；
-
一个椭球体的名字。长半轴(semi-major axis)和反扁率(inverse flattening)；
-
本初子午线(prime meridian)名和其与格林威治子午线的偏移值；
-
投影的方法类型（如横轴莫卡托）；
-
投影的参数列表（如中央经线等）；
-
一个单位的名称和其到米和弧度单位的转换参数；
-
轴线的名称和顺序；
-
在预定义的权威坐标系中的编码（如EPSG）。
tif1 = gdal.Open('I:\\MasterLife\\资源与环境遥感\\实验作业一\\caijian2.tif')
print(tif1.GetProjection())
print(tif1.GetProjectionRef())
两个获得的结果一致，均输出栅格影像的投影坐标系。
地理仿式变换参数。
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 */
坐标转换
常规的WKT格式。
__author__ = 'He XK'
from osgeo import gdal, osr
from osgeo.gdalconst import *
ds = gdal.Open('I:\\MasterLife\\资源与环境遥感\\实验作业一\\caijian2.tif')
sr = ds.GetProjection()
print(sr)
# 声明一个osr的对象
osr_obj = osr.SpatialReference()
# 将获取的空间参考系导入到osr对象中
# 此时 osr对象也就具备了投影的所有信息
osr_obj.ImportFromWkt(sr)
# 常规的wkt格式
print(osr_obj)
接上代码。
osr_obj.MorphToESRI()将wkt格式转化成ESRI的格式。
输出结果会有一点变化。
osr_obj.MorphToESRI()
print(osr_obj)
自动识别EPSG。
# 对于WGS84
print(osr_obj.AutoIdentifyEPSG()) # WGS84的排序是0
# 对于Beijing 1954
print(sr2.AutoIdentifyEPSG()) # 7
判断两个SR对象是否相同。
0是不相同。
print(sr2.IsSame(sr)) # 0
导出成米单位的坐标系统。
因为osr_obj是遥感影像已有坐标系，所以具备经纬度信息，自定义不确定到具体带的没有经纬度信息。
print(osr_obj.ExportToMICoordSys())
# Earth Projection 8, 104, "m", 117, 0, 0.9996, 500000, 0
导出成PCI的格式，我也不知道PCI是什么格式。
print(osr_obj.ExportToPCI())
# ['UTM 50 E012', 'METRE', (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)]
还有其他的一些导出格式。
ExportToWkt()
ExportToPrettyWkt()
ExportToProj4()
ExportToPCI()
ExportToUSGS()
ExportToXML()
矢量数据和栅格数据获取空间参考系的方法
矢量数据：
layer.GetSpatialRef()
矢量数据获取的空间坐标系是OSR对象。
栅格数据：
projection = dataset.GetProjection() # 或者
projection = dataset.GetProjectionRef()
geotransform = dataset.GetGeoTransform()
注意，以上栅格数据获取的投影坐标系信息，都不是OSR对象。
地理坐标系和投影坐标系之间的坐标转换
相当于之前使用过的重投影。
我们可以创建一个新的坐标系统，对某一个SHP文件进行重投影。
下面我们进行一个对shp文件重投影的实例。
# -*- coding: utf-8 -*-
# 学习时间： 2022/11/29 16:47
__author__ = 'He XK'
import os
from osgeo import osr, ogr, gdal
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
gdal.SetConfigOption("SHAPE_ENCODING", "UTF8")
driver = ogr.GetDriverByName('ESRI SHAPEFILE')
origin_path = 'shp/province_1.shp'
output_path = 'shp/demo32_reproject.shp'
dataset = driver.Open(origin_path)
layer = dataset.GetLayer(0)
origin_ref = layer.GetSpatialRef() # WGS84
# 定义新的投影
sr = osr.SpatialReference()
# sr.SetWellKnownGeogCS('EPSG:4610')
# sr.SetProjCS('xian 80')
# sr.SetUTM(19, True)
sr.ImportFromEPSG(3857)
transform = osr.CoordinateTransformation(origin_ref, sr)
if os.path.exists(output_path):
driver.DeleteDataSource(output_path)
else:
print('THIS SHP IS NOT EXIST!')
output_ds = driver.CreateDataSource(output_path)
output_layer = output_ds.CreateLayer('repro', srs=sr, geom_type=ogr.wkbPolygon)
output_feature_defn = output_layer.GetLayerDefn()
# 创建字段
output_field_def = ogr.FieldDefn('Name', ogr.OFTString)
output_field_def.SetWidth(10)
output_layer.CreateField(output_field_def)
layer.ResetReading()
i = 0
while i < layer.GetFeatureCount():
feature = layer.GetFeature(i)
i += 1
geom = feature.GetGeometryRef()
output_geom = geom.Transform(transform)
output_feature = ogr.Feature(output_feature_defn)
output_feature.SetGeometry(geom)
output_feature.SetField('Name', feature.GetField('Name'))
output_layer.CreateFeature(output_feature)
output_ds.Release()
dataset.Release()
put_feature)
output_ds.Release()
dataset.Release()
其中只能对geometry进行重投影。
那么我们就要获取到geometry然后投影，之后创建新要素挂载到新图层。
从QGIS上并不能直接看到不同投影的变化，可能是因为QGIS自己进行了重投影，但是查看坐标系可以发现投影坐标系。
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
