导致在转换时会出错的两个原因：
第一，OGRCreateCoordinateTransformation() 可能失败。
一般是因为程序认为两个指定坐标系统不能建立转换关系，
这可能是由于使用的投影Proj内部暂时不支持，
这样两个不一样的基准面就没有已知关系。 或者一个投影定义得不适当。
如果OGRCreateCoordinateTransformation() 失败了将返回空。
OGRCoordinateTransformation::Transform() 方法本身也可能失败。
这可能是因为一个有问题积累后形成的问题。
或者是“因为一个或者多个点因数值上未定义而省略操作”起因的后果。
Transform()函数如果成功，则返回为TRUE，
如果某个点转换失败，剩下的点就会以一个错误的状态保持。
坐标转换服务其实可以处理3D点的。 并且服务会根据椭球体和基准面上的高程差异调节高程。 将来，也有可能会在不同矢量基准面上进行的点间的转换的应用。 如果没有Z值，则假设点都在大地水准面上。
下面的例子显示了如何方便创建一个地理坐标系统 到一个基于同一个地理坐标系统的投影坐标系统。并在两个坐标系统之间转换。
In [1]:
from osgeo import gdal
from osgeo import osr
# import osr
dataset = gdal.Open("/gdata/geotiff_file.tif")
从数据集中获取空间参考并且建立一个SpatialReference对象：
In [2]:
sr = dataset.GetProjectionRef()
osrobj = osr.SpatialReference()
osrobj.ImportFromWkt(sr)
Out[2]:
In [3]:
osrobj.ExportToWkt()
Out[3]:
In [4]:
osrobj.MorphToESRI() # 重点，转成Esri的wkt格式
osrobj.ExportToWkt()
Out[4]:
In [5]:
osrobj.IsGeographic()
Out[5]:
In [6]:
osrobj.IsProjected()
Out[6]:
再获取一个进行比较：
In [7]:
dataset2 = gdal.Open("/gdata/lu75c.tif")
sr2 = dataset2.GetProjectionRef()
osrobj2 = osr.SpatialReference()
osrobj2.ImportFromWkt(sr2)
osrobj2.IsSame(osrobj)
Out[7]:
创建一个地理坐标系，然后和已有的坐标系统进行比较：
In [8]:
osrobj3 = osr.SpatialReference()
osrobj3.SetWellKnownGeogCS("WGS84")
osrobj3.IsSame(osrobj2)
osrobj3.IsSame(osrobj)
osrobj3.ExportToWkt()
osrobj3.IsGeographic()
Out[8]:
In [9]:
# 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()
Out[9]:
In [10]:
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()
Out[10]:
In [11]:
spatialRef.ExportToPrettyWkt()
Out[11]:
In [12]:
spatialRef.ExportToPCI()
Out[12]:
In [13]:
spatialRef.ExportToUSGS()
Out[13]:
In [14]:
spatialRef.ExportToXML()
Out[14]:
In [15]:
from osgeo import ogr, osr
spatialRef = osr.SpatialReference()
spatialRef.ImportFromEPSG(26912)
spatialRef.MorphToESRI()
file = open('yourshpfile.prj', 'w')
file.write(spatialRef.ExportToWkt())
file.close()
help(ct.TransformPoint)用CoordinateTransformation转换的结果：
In [16]:
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()
Out[16]:
In [17]:
osrobj3.IsGeographic()
Out[17]:
In [18]:
ct = osr.CoordinateTransformation(osrobj,osrobj3)
# ct.TransformPoint([590000.0,4928000.0,0])
# ct.TransformPoint(609000,4928000)
# ct.TransformPoint(609000,4914000)
# ct.TransformPoint(590000,4914000)
ArcGIS里面的地理坐标范围的结果：
In decimal degrees
West: -103.870326
East: -103.628957
North: 44.501659
South: 44.373036用proj计算的和用ArcGIS计算的有些差距(分别有一个很准，一个差了0.01度)。但是我觉得这种误差应该是被允许的。
In [19]:
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()
Out[19]:
In [20]:
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
In [21]:
import shapefile
uu = shapefile.Reader('/tmp/xx_ab.shp')
In [22]:
import helper; helper.info()