目录
空间参考系统
OGC (Open Geospatial Consortium) 空间参考信息是一种用于定义地理数据空间坐标系统的标准。它允许不同系统之间进行地理数据的交换和比较,确保数据的一致性和准确性。以下是一些关于OGC空间参考信息的关键点:
-
坐标参考系统 (CRS):CRS 是定义地球表面上点的坐标位置的系统。它可以是地理的(基于经纬度)或投影的(基于某种数学投影)。
-
EPSG代码:EPSG(欧洲石油调查组)代码是用于唯一标识CRS的代码。每个EPSG代码都对应一个特定的坐标系统。
-
WKT:WKT(Well-Known Text)是一种文本表示法,用于描述CRS。它包含了坐标系统的所有必要信息,如单位、轴顺序、基准面等。
-
PROJ:PROJ是一个用于处理地理投影的库,它支持多种坐标转换和CRS定义。
-
GDAL/OGR:GDAL(Geospatial Data Abstraction Library)和OGR(OGR Simple Features Library)是两个流行的开源库,用于读写栅格和矢量地理空间数据格式。OGR是GDAL的一部分,专注于矢量数据。
-
空间数据框架:空间数据框架是一组规则,定义了空间数据如何组织和存储。例如,UTM(Universal Transverse Mercator)是一种常用的投影框架。
-
转换和变换:在不同的CRS之间转换地理数据时,可能需要进行坐标转换或变换。这涉及到将数据从一个系统映射到另一个系统,可能还需要考虑地球曲率等因素。
-
元数据:空间参考信息通常作为地理数据的元数据存储,确保数据的地理上下文和准确性。
1.1 获取空间参考
from osgeo import osr
ds = ogr.Open(os.path.join(data_dir, 'US', 'states_48.shp'))
# 获取WKT格式的空间参考文本
srs = ds.GetLayer().GetSpatialRef()
# Well Known Text (WKT)
print(srs)
# PROJ.4,将WKT格式文本转换为Proj4
print(srs.ExportToProj4())
# XML 转换为XML格式
print(srs.ExportToXML())
# Look at a UTM SRS.
utm_sr = osr.SpatialReference()
# EPSG UTM 数值:26912
# 查询网址:https://epsg.io/
# 介绍网址:https://blog.csdn.net/qq_28419035/article/details/141067699
utm_sr.ImportFromEPSG(26912)
print(utm_sr) # WKT
print(utm_sr.ExportToProj4()) # PROJ.4
print(utm_sr.ExportToXML()) # XML
# Get the projection name.
print(utm_sr.GetAttrValue('PROJCS'))
# 获取授权信息[字典形式]
print(utm_sr.GetAttrValue('AUTHORITY'))
print(utm_sr.GetAttrValue('AUTHORITY', 1))
# Get the datum code. 地理参考基准
print(utm_sr.GetAuthorityCode('DATUM'))
# Get the false easting.
print(utm_sr.GetProjParm(osr.SRS_PP_FALSE_EASTIN
1.2 创建空间参考对象
# Create a UTM SRS from an EPSG code.
sr = osr.SpatialReference()
sr.ImportFromEPSG(26912)
print(sr.GetAttrValue('PROJCS'))
# Create a UTM SRS from a PROJ.4 string.
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'))
# Create a unprojected SRS from a WKT string.
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)
# Create an Albers SRS using parameters.
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)
1.3 为创建的适量添加空间参考
# Make sure that the output folder exists in your data directory before
# trying this example.
out_fn = os.path.join('testdata.shp')
# Create an empty shapefile that uses a UTM SRS. If you run this it will
# create the shapefile with a .prj file containing the SRS info.
sr = osr.SpatialReference()
sr.ImportFromEPSG(26912)
ds = ogr.GetDriverByName('ESRI Shapefile').CreateDataSource(out_fn)
lyr = ds.CreateLayer('counties', sr, ogr.wkbPolygon)
1.4 重投影
# Get the world landmasses and plot them.
world = pb.get_shp_geom(os.path.join('.shp'))
# Create a point for the Eiffel Tower.
tower = ogr.Geometry(wkt='POINT (2.294694 48.858093)')
# 为坐标点赋予空间参考:WGS84
tower.AssignSpatialReference(osr.SpatialReference(osr.SRS_WKT_WGS84))
# 重投影为 Web Mercator
web_mercator_sr = osr.SpatialReference()
web_mercator_sr.ImportFromEPSG(3857)
# 对几何对象进行重投影(对点坐标进行操作):会报错
world.TransformTo(web_mercator_sr)
# 设置参数,重新进行投影
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)
# 根据坐标转换参数进行重投影
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)
# Create an unprojected NAD27 SRS and add datum shift info.
# 创建一个NAD27投影坐标系统,并添加偏移信息
sr = osr.SpatialReference()
sr.SetWellKnownGeogCS('NAD27')
sr.SetTOWGS84(-8, 160, 176)
1.5Pyproj 库
# Transform lat/lon to UTM.
import pyproj
utm_proj = pyproj.Proj('+proj=utm +zone=31 +ellps=WGS84')
x, y = utm_proj(2.294694, 48.858093)
print(x, y)
# 将坐标转为经纬度:投影转为椭球
# Go back to lat/lon.
x1, y1 = utm_proj(x, y, inverse=True)
print(x1, y1)
# Convert UTM WGS84 coordinates to UTM NAD27.
wgs84 = pyproj.Proj('+proj=utm +zone=18 +datum=WGS84')
nad27 = pyproj.Proj('+proj=utm +zone=18 +datum=NAD27')
# NAD27 投影基准 :580744.32, 4504695.26
x, y = pyproj.transform(wgs84, nad27, 580744.32, 4504695.26)
print(x, y)
# 计算大远距离
#r Los Angeles经纬度 and Berlin经纬度.
la_lat, la_lon = 34.0500, -118.2500
berlin_lat, berlin_lon = 52.5167, 13.3833
# 创建地理基准
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))
# 根据起始坐标、方位、举例来计算目标的经纬度与方位
x, y, bearing = geod.fwd(berlin_lon, berlin_lat, back, dist)
print('{}, {}\n{}'.format(x, y, bearing))
coords = geod.npts(la_lon, la_lat, berlin_lon, berlin_lat, 100)
# Only print the first 3.
for i in range(3):
print(coords[i])
1.6矢量文件重投影
# Script to reproject a shapefile.
from osgeo import ogr, osr
# 将shp文件投影到新的坐标系下
# Create an output 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''')
# Don't forget to change your directory here.
ds = ogr.Open(r'D:\osgeopy-data\US', 1)
# Get the input layer.
in_lyr = ds.GetLayer('us_volcanos')
# Create the empty output layer.
out_lyr = ds.CreateLayer('us_volcanos_aea', sr,
ogr.wkbPoint)
out_lyr.CreateFields(in_lyr.schema)
# Loop through the features in the input layer.
out_feat = ogr.Feature(out_lyr.GetLayerDefn())
for in_feat in in_lyr:
# Clone the geometry, project it, and add it to the feature.
geom = in_feat.geometry().Clone()
geom.TransformTo(sr)
out_feat.SetGeometry(geom)
# Copy attributes.
for i in range(in_feat.GetFieldCount()):
out_feat.SetField(i, in_feat.GetField(i))
# Insert the feature
out_lyr.CreateFeature(out_feat)
总结
空间参考系统的坐标基准为参考椭球,如WGS84。数据投影的操作对象是坐标点,因此投影处理的基本单元是geometry;投影变换有两种方式1:Geo.TransformTo();2:计算两个投影之间的转换参数:
# 计算转换参数
# source->target
ct = osr.CoordinateTransformation(web_mercator_sr, peters_sr)
# 调用转换方法进行转换
world.Transform(ct)