目录

上一个主题

4.4. osr 模块简介与用法

下一个主题

5. 矢量数据的空间分析:使用Shapely

关注公众号


常见问题

  1. Windows下的安装说明
  2. Jupyter免费在线实验环境
  3. 勘误与补充


>>> from helper import info; info()
页面更新时间: 2020-03-28 11:10:08
操作系统/OS: Linux-4.19.0-8-amd64-x86_64-with-debian-10.3 ;Python: 3.7.3

4.5. 坐标转换

OGRCoordinateTransformation 类可用在不同坐标系统间转换坐标。 新的转换对象由OGRCoordinateTransformation()创建。 创建后 OGRCoordinateTransformation::Transform() 方法就可以用来转换不同坐标系的点。

导致在转换时会出错的两个原因:

第一,OGRCreateCoordinateTransformation() 可能失败。 一般是因为程序认为两个指定坐标系统不能建立转换关系, 这可能是由于使用的投影Proj内部暂时不支持, 这样两个不一样的基准面就没有已知关系。 或者一个投影定义得不适当。 如果OGRCreateCoordinateTransformation() 失败了将返回空。

OGRCoordinateTransformation::Transform() 方法本身也可能失败。 这可能是因为一个有问题积累后形成的问题。 或者是“因为一个或者多个点因数值上未定义而省略操作”起因的后果。 Transform()函数如果成功,则返回为TRUE, 如果某个点转换失败,剩下的点就会以一个错误的状态保持。

坐标转换服务其实可以处理3D点的。 并且服务会根据椭球体和基准面上的高程差异调节高程。 将来,也有可能会在不同矢量基准面上进行的点间的转换的应用。 如果没有Z值,则假设点都在大地水准面上。

下面的例子显示了如何方便创建一个地理坐标系统 到一个基于同一个地理坐标系统的投影坐标系统。并在两个坐标系统之间转换。

4.5.1. 从栅格数据中得到

>>> from osgeo import gdal
>>> from osgeo import osr
>>> # import osr
>>> dataset = gdal.Open("/gdata/geotiff_file.tif")

从数据集中获取空间参考并且建立一个SpatialReference对象:

>>> sr = dataset.GetProjectionRef()
>>> osrobj = osr.SpatialReference()
>>> osrobj.ImportFromWkt(sr)
0

4.5.2. 输出格式

>>> osrobj.ExportToWkt()
'PROJCS["Albers_Beijing54",GEOGCS["GCS_Krasovsky_1940",DATUM["Not_specified_based_on_Krassowsky_1940_ellipsoid",SPHEROID["Krassowsky 1940",6378245,298.2999999999998,AUTHORITY["EPSG","7024"]],AUTHORITY["EPSG","6024"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",47],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",105],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]'
>>> osrobj.MorphToESRI() # 重点,转成Esri的wkt格式
0
>>> osrobj.ExportToWkt()
'PROJCS["Albers_Beijing54",GEOGCS["GCS_Krasovsky_1940",DATUM["D_Krasovsky_1940",SPHEROID["Krasovsky_1940",6378245,298.3]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Albers"],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",47],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",105],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]'

4.5.3. 其他操作

>>> osrobj.IsGeographic()
0
>>> osrobj.IsProjected()
1

再获取一个进行比较:

>>> dataset2 = gdal.Open("/gdata/lu75c.tif")
>>> sr2 = dataset2.GetProjectionRef()
>>> osrobj2 = osr.SpatialReference()
>>> osrobj2.ImportFromWkt(sr2)
>>> osrobj2.IsSame(osrobj)
0

创建一个地理坐标系,然后和已有的坐标系统进行比较:

>>> osrobj3 = osr.SpatialReference()
>>> osrobj3.SetWellKnownGeogCS("WGS84")
>>> osrobj3.IsSame(osrobj2)
>>> osrobj3.IsSame(osrobj)
>>> osrobj3.ExportToWkt()
'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"]]'
>>> osrobj3.ExportToMICoordSys()
'Earth Projection 1, 104'
>>> osrobj3.ExportToPCI()
['LONG/LAT    D000',
 'DEGREE',
 (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)]
>>> osrobj3.ExportToUSGS()
[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),
 12]
>>> osrobj3.ExportToProj4()
'+proj=longlat +datum=WGS84 +no_defs '
>>> osrobj3.IsGeographic()
1

proj.4 字符串转成Esri wkt字符串的方法 April 12, 2016barry.z proj.4 字符串转成Esri wkt字符串的方法 使用gdal/ogr库进行转换

>>> import os
>>> import sys
>>> import string
>>> import osgeo.osr
>>>
>>> srs = osgeo.osr.SpatialReference()
>>> srs.ImportFromProj4(sys.argv[1])
>>> srs.MorphToESRI() # 重点,转成Esri的wkt格式
>>> srs.ExportToWkt()
''

4.5.4. 从矢量数据中获取

>>> from osgeo import ogr, osr
>>> driver = ogr.GetDriverByName('ESRI Shapefile')
>>> dataset = driver.Open('/gdata/region_popu.shp')
>>> layer = dataset.GetLayer()
>>> spatialRef = layer.GetSpatialRef()
>>>
>>> spatialRef.ExportToWkt()
'PROJCS["WGS 84 / Pseudo-Mercator",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["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]'
>>> spatialRef.ExportToPrettyWkt()
'PROJCS["WGS 84 / Pseudo-Mercator",n    GEOGCS["WGS 84",n        DATUM["WGS_1984",n            SPHEROID["WGS 84",6378137,298.257223563,n                AUTHORITY["EPSG","7030"]],n            AUTHORITY["EPSG","6326"]],n        PRIMEM["Greenwich",0,n            AUTHORITY["EPSG","8901"]],n        UNIT["degree",0.0174532925199433,n            AUTHORITY["EPSG","9122"]],n        AUTHORITY["EPSG","4326"]],n    PROJECTION["Mercator_1SP"],n    PARAMETER["central_meridian",0],n    PARAMETER["scale_factor",1],n    PARAMETER["false_easting",0],n    PARAMETER["false_northing",0],n    UNIT["metre",1,n        AUTHORITY["EPSG","9001"]],n    AXIS["X",EAST],n    AXIS["Y",NORTH],n    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],n    AUTHORITY["EPSG","3857"]]'
>>> spatialRef.ExportToPCI()
['MER         D000',
 'METRE',
 (0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  1.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0,
  0.0)]
>>> spatialRef.ExportToUSGS()
[5,
 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),
 12]
>>> spatialRef.ExportToXML()
'<gml:ProjectedCRS gml:id="ogrcrs1">n  <gml:srsName>WGS 84 / Pseudo-Mercator</gml:srsName>n  <gml:srsID>n    <gml:name codeSpace="urn:ogc:def:crs:EPSG::">3857</gml:name>n  </gml:srsID>n  <gml:baseCRS>n    <gml:GeographicCRS gml:id="ogrcrs2">n      <gml:srsName>WGS 84</gml:srsName>n      <gml:srsID>n        <gml:name codeSpace="urn:ogc:def:crs:EPSG::">4326</gml:name>n      </gml:srsID>n      <gml:usesEllipsoidalCS>n        <gml:EllipsoidalCS gml:id="ogrcrs3">n          <gml:csName>ellipsoidal</gml:csName>n          <gml:csID>n            <gml:name codeSpace="urn:ogc:def:cs:EPSG::">6402</gml:name>n          </gml:csID>n          <gml:usesAxis>n            <gml:CoordinateSystemAxis gml:id="ogrcrs4" gml:uom="urn:ogc:def:uom:EPSG::9102">n              <gml:name>Geodetic latitude</gml:name>n              <gml:axisID>n                <gml:name codeSpace="urn:ogc:def:axis:EPSG::">9901</gml:name>n              </gml:axisID>n              <gml:axisAbbrev>Lat</gml:axisAbbrev>n              <gml:axisDirection>north</gml:axisDirection>n            </gml:CoordinateSystemAxis>n          </gml:usesAxis>n          <gml:usesAxis>n            <gml:CoordinateSystemAxis gml:id="ogrcrs5" gml:uom="urn:ogc:def:uom:EPSG::9102">n              <gml:name>Geodetic longitude</gml:name>n              <gml:axisID>n                <gml:name codeSpace="urn:ogc:def:axis:EPSG::">9902</gml:name>n              </gml:axisID>n              <gml:axisAbbrev>Lon</gml:axisAbbrev>n              <gml:axisDirection>east</gml:axisDirection>n            </gml:CoordinateSystemAxis>n          </gml:usesAxis>n        </gml:EllipsoidalCS>n      </gml:usesEllipsoidalCS>n      <gml:usesGeodeticDatum>n        <gml:GeodeticDatum gml:id="ogrcrs6">n          <gml:datumName>WGS_1984</gml:datumName>n          <gml:datumID>n            <gml:name codeSpace="urn:ogc:def:datum:EPSG::">6326</gml:name>n          </gml:datumID>n          <gml:usesPrimeMeridian>n            <gml:PrimeMeridian gml:id="ogrcrs7">n              <gml:meridianName>Greenwich</gml:meridianName>n              <gml:meridianID>n                <gml:name codeSpace="urn:ogc:def:meridian:EPSG::">8901</gml:name>n              </gml:meridianID>n              <gml:greenwichLongitude>n                <gml:angle uom="urn:ogc:def:uom:EPSG::9102">0</gml:angle>n              </gml:greenwichLongitude>n            </gml:PrimeMeridian>n          </gml:usesPrimeMeridian>n          <gml:usesEllipsoid>n            <gml:Ellipsoid gml:id="ogrcrs8">n              <gml:ellipsoidName>WGS 84</gml:ellipsoidName>n              <gml:ellipsoidID>n                <gml:name codeSpace="urn:ogc:def:ellipsoid:EPSG::">7030</gml:name>n              </gml:ellipsoidID>n              <gml:semiMajorAxis uom="urn:ogc:def:uom:EPSG::9001">6378137</gml:semiMajorAxis>n              <gml:secondDefiningParameter>n                <gml:inverseFlattening uom="urn:ogc:def:uom:EPSG::9201">298.257223563</gml:inverseFlattening>n              </gml:secondDefiningParameter>n            </gml:Ellipsoid>n          </gml:usesEllipsoid>n        </gml:GeodeticDatum>n      </gml:usesGeodeticDatum>n    </gml:GeographicCRS>n  </gml:baseCRS>n  <gml:definedByConversion>n    <gml:Conversion gml:id="ogrcrs9">n      <gml:coordinateOperationName>Mercator_1SP</gml:coordinateOperationName>n    </gml:Conversion>n  </gml:definedByConversion>n  <gml:usesCartesianCS>n    <gml:CartesianCS gml:id="ogrcrs10">n      <gml:csName>Cartesian</gml:csName>n      <gml:csID>n        <gml:name codeSpace="urn:ogc:def:cs:EPSG::">4400</gml:name>n      </gml:csID>n      <gml:usesAxis>n        <gml:CoordinateSystemAxis gml:id="ogrcrs11" gml:uom="urn:ogc:def:uom:EPSG::9001">n          <gml:name>Easting</gml:name>n          <gml:axisID>n            <gml:name codeSpace="urn:ogc:def:axis:EPSG::">9906</gml:name>n          </gml:axisID>n          <gml:axisAbbrev>E</gml:axisAbbrev>n          <gml:axisDirection>east</gml:axisDirection>n        </gml:CoordinateSystemAxis>n      </gml:usesAxis>n      <gml:usesAxis>n        <gml:CoordinateSystemAxis gml:id="ogrcrs12" gml:uom="urn:ogc:def:uom:EPSG::9001">n          <gml:name>Northing</gml:name>n          <gml:axisID>n            <gml:name codeSpace="urn:ogc:def:axis:EPSG::">9907</gml:name>n          </gml:axisID>n          <gml:axisAbbrev>N</gml:axisAbbrev>n          <gml:axisDirection>north</gml:axisDirection>n        </gml:CoordinateSystemAxis>n      </gml:usesAxis>n    </gml:CartesianCS>n  </gml:usesCartesianCS>n</gml:ProjectedCRS>n'

4.5.5. 写入到文件

>>> from osgeo import ogr, osr
>>>
>>> spatialRef = osr.SpatialReference()
>>> spatialRef.ImportFromEPSG(26912)
>>>
>>> spatialRef.MorphToESRI()
>>> file = open('yourshpfile.prj', 'w')
>>> file.write(spatialRef.ExportToWkt())
>>> file.close()

4.5.6. 地理坐标系和投影坐标系之间的坐标转换

spot的图像投影坐标范围:

In projected or local coordinates
Left: 590000.000000
Right: 609000.000000
Top: 4928000.000000
Bottom: 4914000.000000
help(ct.TransformPoint)

用CoordinateTransformation转换的结果:

>>> import osr
>>> import gdal
>>> dataset = gdal.Open("/gdata/geotiff_file.tif")

从数据集中获取空间参考并且建立一个SpatialReference对象

>>> sr = dataset.GetProjectionRef()
>>> osrobj = osr.SpatialReference()
>>>
>>> osrobj3 = osr.SpatialReference()
>>> osrobj3.SetWellKnownGeogCS("WGS84")
>>> osrobj3.IsSame(osrobj)
>>> osrobj3.ExportToWkt()
'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"]]'
>>> osrobj3.IsGeographic()
1
>>> ct = osr.CoordinateTransformation(osrobj,osrobj3)
>>>
>>> # ct.TransformPoint([590000.0,4928000.0,0])
>>> # ct.TransformPoint(609000,4928000)
>>> # ct.TransformPoint(609000,4914000)
>>> # ct.TransformPoint(590000,4914000)

下面的代码,展示了从 EPSG 2927 到 EPSG 4326 的投影转换:

>>> from osgeo import ogr
>>> from osgeo import osr
>>>
>>>
>>> source = osr.SpatialReference()
>>> source.ImportFromEPSG(2927)
>>>
>>> target = osr.SpatialReference()
>>> target.ImportFromEPSG(4326)
>>>
>>> transform = osr.CoordinateTransformation(source, target)
>>>
>>> point = ogr.CreateGeometryFromWkt("POINT (1120351.57 741921.42)")
>>> point.Transform(transform)
>>>
>>> point.ExportToWkt()
'POINT (-122.598135130878 47.3488013793486)'

4.5.7. 转换Shapefile的实例

>>> from osgeo import ogr, osr
>>> import os
>>>
>>> driver = ogr.GetDriverByName('ESRI Shapefile')
>>>
>>> # input SpatialReference
>>> inSpatialRef = osr.SpatialReference()
>>> inSpatialRef.ImportFromEPSG(4326)
>>>
>>> # output SpatialReference
>>> outSpatialRef = osr.SpatialReference()
>>> outSpatialRef.ImportFromEPSG(3857)
>>>
>>> # create the CoordinateTransformation
>>> coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)
>>>
>>> # get the input layer
>>> inDataSet = driver.Open('/gdata/region_popu.shp')
>>> inLayer = inDataSet.GetLayer()
>>>
>>> # create the output layer
>>> outputShapefile = '/tmp/xx_ab.shp'
>>> if os.path.exists(outputShapefile):
>>>     driver.DeleteDataSource(outputShapefile)
>>> outDataSet = driver.CreateDataSource(outputShapefile)
>>> outLayer = outDataSet.CreateLayer("basemap_4326", geom_type=ogr.wkbMultiPolygon)
>>>
>>> # add fields
>>> inLayerDefn = inLayer.GetLayerDefn()
>>> for i in range(0, inLayerDefn.GetFieldCount()):
>>>     fieldDefn = inLayerDefn.GetFieldDefn(i)
>>>     outLayer.CreateField(fieldDefn)
>>>
>>> # get the output layer's feature definition
>>> outLayerDefn = outLayer.GetLayerDefn()
>>>
>>> # loop through the input features
>>> inFeature = inLayer.GetNextFeature()
>>> while inFeature:
>>>     # get the input geometry
>>>     geom = inFeature.GetGeometryRef()
>>>     # reproject the geometry
>>>     geom.Transform(coordTrans)
>>>     # create a new feature
>>>     outFeature = ogr.Feature(outLayerDefn)
>>>     # set the geometry and attribute
>>>     outFeature.SetGeometry(geom)
>>>     for i in range(0, outLayerDefn.GetFieldCount()):
>>>         outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
>>>     # add the feature to the shapefile
>>>     outLayer.CreateFeature(outFeature)
>>>     # dereference the features and get the next input feature
>>>     outFeature = None
>>>     inFeature = inLayer.GetNextFeature()
>>>
>>> # Save and close the shapefiles
>>> inDataSet = None
>>> outDataSet = None