创建和比较

从栅格数据中得到

In [2]:
from osgeo import gdal
from osgeo import osr
# import osr
dataset = gdal.Open("/gdata/")
# 从数据集中获取空间参考并且建立一个SpatialReference对象
sr = dataset.GetProjectionRef()
osrobj = osr.SpatialReference()
osrobj.ImportFromWkt(sr)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-2-1bde19fb4b37> in <module>()
      4 dataset = gdal.Open("/gdata/K52E015007.tif")
      5 # 从数据集中获取空间参考并且建立一个SpatialReference对象
----> 6 sr = dataset.GetProjectionRef()
      7 osrobj = osr.SpatialReference()
      8 osrobj.ImportFromWkt(sr)

AttributeError: 'NoneType' object has no attribute 'GetProjectionRef'

输出格式

In [ ]:
osrobj.ExportToWkt()
In [ ]:
osrobj.MorphToESRI() # 重点,转成Esri的wkt格式
osrobj.ExportToWkt()

其他操作

In [ ]:
osrobj.IsGeographic()
In [ ]:
osrobj.IsProjected()

再获取一个进行比较:

In [ ]:
dataset2 = gdal.Open("/gdata/lu75c.tif")
sr2 = dataset2.GetProjectionRef()
osrobj2 = osr.SpatialReference()
osrobj2.ImportFromWkt(sr2)
osrobj2.IsSame(osrobj)

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

In [ ]:
osrobj3 = osr.SpatialReference()
osrobj3.SetWellKnownGeogCS("WGS84")
osrobj3.IsSame(osrobj2)
osrobj3.IsSame(osrobj)
osrobj3.ExportToWkt()
osrobj3.IsGeographic()
In [ ]:
# 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()

从矢量数据中获取

In [ ]:
from osgeo import ogr, osr
driver = ogr.GetDriverByName('ESRI Shapefile')
dataset = driver.Open('/gdata/world_borders.shp')
layer = dataset.GetLayer()
spatialRef = layer.GetSpatialRef()

spatialRef.ExportToWkt()
In [ ]:
spatialRef.ExportToPrettyWkt()
In [ ]:
spatialRef.ExportToPCI()
In [ ]:
spatialRef.ExportToUSGS()
In [ ]:
spatialRef.ExportToXML()

写入到文件

In [ ]:
from osgeo import ogr, osr

spatialRef = osr.SpatialReference()
spatialRef.ImportFromEPSG(26912)

spatialRef.MorphToESRI()
file = open('yourshpfile.prj', 'w')
file.write(spatialRef.ExportToWkt())
file.close()

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

spot的图像投影坐标范围:

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

用CoordinateTransformation转换的结果:

In [ ]:
import osr
import gdal
dataset = gdal.Open("/gdata/K52E015007.tif")
# 从数据集中获取空间参考并且建立一个SpatialReference对象
sr = dataset.GetProjectionRef()
osrobj = osr.SpatialReference()

osrobj3 = osr.SpatialReference()
osrobj3.SetWellKnownGeogCS("WGS84")
osrobj3.IsSame(osrobj)
osrobj3.ExportToWkt()
In [ ]:
osrobj3.IsGeographic()
In [ ]:
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 [ ]:
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()

转换Shapefile的实例

In [ ]:
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/world_borders.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 [ ]:
import shapefile
uu = shapefile.Reader('/tmp/xx_ab.shp')
In [ ]:
 

原型

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

原型:

OGRSpatialReference oSourceSRS, oTargetSRS;
OGRCoordinateTransformation *poCT;
double  x, y;
oSourceSRS.importFromEPSG( atoi(papszArgv[i+1]) );
oTargetSRS.importFromEPSG( atoi(papszArgv[i+2]) );
poCT = OGRCreateCoordinateTransformation( &oSourceSRS,
                      &oTargetSRS );
x = atof( papszArgv[i+3] );
y = atof( papszArgv[i+4] );
if( poCT == NULL || !poCT->Transform( 1, &x, &y ) )
  printf( "Transformation failed.\n" );
else
  printf( "(%f,%f) -> (%f,%f)\n",
    atof( papszArgv[i+3] ),
    atof( papszArgv[i+4] ),
    x, y );

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

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

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

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

实例

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

OGRSpatialReference  oUTM, *poLatLong;
OGRCoordinateTransformation *poTransform;
oUTM.SetProjCS("UTM 17 / WGS84");
oUTM.SetWellKnownGeogCS( "WGS84" );
oUTM.SetUTM( 17 );
poLatLong = oUTM.CloneGeogCS();
poTransform = OGRCreateCoordinateTransformation( &oUTM, poLatLong );
if( poTransform == NULL )
{
  ...
}
...
if( !poTransform->Transform( nPoints, x, y, z ) )
...
In [ ]: