OGR-空间参考

空间参考系统

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)