osgeo：OSR

osgeo：OSR

from osgeo import osr
# 初始化osr.SpatialReference对象
srs = osr.SpatialReference()

wkt格式

一、定义地理坐标系

1.1 标准系统定义

srs.SetWellKnownGeoCS(‘WGS84’)

• 一些标准的坐标系统：“NAD27”， “NAD83”， “WGS72” 和 “WGS84”

1.2 自定义

srs.SetGeogCS(GeoName, DatumName, EllipsoidName, SemiMajor, InvFlattening, PMName=‘Greenwich’, PMOffset=0.0, ConvertToRadians=0.0174532925199433)

• GeoName/DatumName/EllipsoidName/PMName:系统名/基本面名/椭球名/子午线名
• SemiMajor：osr.SRS_WGS84_SEMIMAJOR，椭球体
• InvFlattening：osr.SRS_WGS84_INVFLATTENING, 反扁率
• PMOffset：0.0，本初子午线和其与格林威治子午线(Greenwich)的偏移值
• ConvertToRadians：osr.SRS_UA_DEGREE_CONV，角度单位，一般为度

# 标准系统定义
srs = osr.SpatialReference()
srs.SetWellKnownGeogCS('WGS84')
wkt = srs.ExportToPrettyWkt()
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"]]

# 自定义
srs = osr.SpatialReference()
srs.SetGeogCS("My geographic coordinate system", 'WGS_1984', "My WGS84 Spheroid",
              osr.SRS_WGS84_SEMIMAJOR, osr.SRS_WGS84_INVFLATTENING, "Greenwich",
              0.0, osr.SRS_UA_DEGREE_CONV)
print(srs.ExportToPrettyWkt())

GEOGCS["My geographic coordinate system",
    DATUM["WGS_1984",
        SPHEROID["My WGS84 Spheroid",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["0.0174532925199433",0.0174532925199433],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST]]

二、定义投影坐标系

1. srs.SetProjCS(name):设置投影坐标系名称
2. srs.SetWellKnownGeogCS:分配地理坐标系统
3. srs.SetUTM(int zone, int north=1):设置投影转换参数细节
   • SetUTM:(Transverse Mercator)
   • SetLCC:(Lambert Conformal Conic)
   • SetMercator:(Merator)

srs = osr.SpatialReference()
srs.SetProjCS('UTM 17 (WGS84) in northern hemisphere.')
srs.SetWellKnownGeogCS('WGS84')
srs.SetUTM(17, True)
wkt = srs.ExportToPrettyWkt()
print(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["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]

三、查询坐标系统

• IsProjected：是否是投影坐标系
• IsGeographic：是否是地理坐标系
• IsSame(srs2)：参考坐标是否相同
• GetSemiMajior/GetSemiMinor/GetInvFlattening：获取椭球体信息
• GetAttrValue(‘PROJCS’)：获取PROJCS, GEOGCS, DATUM, SPHEROID和PROJECTION的名字表达字符串
• GetProjParm(osr.SRS_PP_**)：获取投影参数, 例如：osr.SRS_PP_FALSE_EASTING
• GetLinearUnits：获取线性单位，并转换成米

srs.IsGeographic(), srs.IsProjected()

(0, 1)

srs.GetAttrValue('PROJCS'), srs.GetProjParm(osr.SRS_PP_FALSE_EASTING), srs.GetLinearUnits()

(‘UTM 17 (WGS84) in northern hemisphere.’, 500000.0, 1.0)

四、投影

• 创建SpatialReference对象：osr.SpatialReference()
• 向SpatialReference对象导入投影信息
  • ImportFromWkt(<wkt>)
  • ImportFromEPSG(<epsg>)
  • ImportFromProj4(<proj4>)
  • ImportFromESRI(<proj_lines>)
• 导出Projection字符串
  • ExportToWkt()
  • ExportToProj4()
  • ExportToPROJJSON() == javascript风格
  • ExportToUSGS
• 转换格式
  • MorphToESRI：转成Esri的wkt格式

4.1 从参考文件创建新投影

from osgeo import osr, ogr
ds = ogr.Open('./osgeo_data/Test.shp')
inLayer = ds.GetLayer(0)
spatialRef = inLayer.GetSpatialRef()
wkt = spatialRef.ExportToWkt()
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
srs.ExportToWkt()

'PROJCS["WGS 84 / UTM zone 50N",GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",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",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32650"]]'

4.2 坐标转换

• CoordinateTransformation(srsSRS, dstSRS)

sourceSRS = osr.SpatialReference()
sourceSRS.ImportFromEPSG(2927)
targetSRS = osr.SpatialReference()
targetSRS.ImportFromEPSG(4326)
# 建立当前坐标系与目标坐标系的坐标转换
coordTrans = osr.CoordinateTransformation(sourceSRS, targetSRS)
point = ogr.CreateGeometryFromWkt("POINT (1120351.57 741921.42)")
point.Transform(coordTrans)
point.ExportToWkt()

'POINT (47.3488013802885 -122.598135130878)'

# 转换shapefile实例
from osgeo import ogr, osr
import os

driver = ogr.GetDriverByName('ESRI Shapefile')
ds_input = driver.Open('./osgeo_data/Test.shp')
inLayer = ds_input.GetLayer(0)
inputSRS = inLayer.GetSpatialRef()
outSRS = osr.SpatialReference()
outSRS.ImportFromEPSG(4326)
coordTrans = osr.CoordinateTransformation(inputSRS, outSRS)
outpath = './osgeo_data/Test_new.shp'
if os.path.exists(outpath):
    driver.DeleteDataSource(outpath)
ds_out = driver.CreateDataSource(outpath)
outLayer = ds_out.CreateLayer('transform', srs=outSRS, geom_type=ogr.wkbMultiPolygon)
inLayerdefn = inLayer.GetLayerDefn()
for i in range(0, inLayerdefn.GetFieldCount()):
    fieldDefn = inLayerdefn.GetFieldDefn(i)
    outLayer.CreateField(fieldDefn)
outLayerDefn = outLayer.GetLayerDefn()
for j in range(inLayer.GetFeatureCount()):
    inFeature = inLayer.GetFeature(j)
    geom = inFeature.GetGeometryRef()
    geom.Transform(coordTrans)
    outFeature = ogr.Feature(outLayerDefn)
    outFeature.SetGeometry(geom)
    for k in range(outLayerDefn.GetFieldCount()):
        outFeature.SetField(outLayerDefn.GetFieldDefn(k).GetNameRef(), inFeature.GetField(k))
    outLayer.CreateFeature(outFeature)
    outFeature = None
ds_input = None
ds_out = None

4.3 将投影写入文件

targetSRS = osr.SpatialReference()
targetSRS.ImportFromEPSG(4326)
targetSRS.MorphToESRI()
file = open('./osgeo_data/write.prj', 'w')
file.write(targetSRS.ExportToWkt())
file.close()