GDAL python教程基础篇（5）OSR添加投影

空间参考（Spatial Reference），实则是给数字附有实际含义。添加了空间参考，才能依据这些数字在地图上定位。同一个Layer需要有相同的空间参考。

1.获取投影

投影使用的是Spatial Reference对象。目前Projections较为多样化，GDAL支持WKT, PROJ.4, ESPG, USGS, ESRI.prj等。我们可以从layer和Geometry中读取Projections：

from osgeo import ogr

inshp = r'D:\pythonlianxi\gdallianxi\Line.shp'
ds = ogr.Open(inshp)
layer = ds.GetLayer(0)
spatialRef = layer.GetSpatialRef()
print(spatialRef)

输出结果

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"]]

Shapefile的投影信息一般存储在.prj文件中，如果没有这个文件，上述函数生成结果为None。

2.创建投影

建立一个新的Projection的方法如下：首先导入osr库，使用osr.SpatialReference()创建SpatialReference对象，再向SpatialReference对象导入投影信息：

ImportFromWkt()
ImportFromEPSG()
ImportFromProj4()
ImportFromESRI(<proj_lines>)
ImportFromPCI(, , )
ImportFromUSGS(<proj_code>, )
ImportFromXML()

使用下列语句导出Projection字符串：

ExportToWkt()
ExportToPrettyWkt()
ExportToProj4()
ExportToPCI()
ExportToUSGS()
ExportToXML()

3.根据已有的投影创建新投影

创建空间参考最简单的办法是用ImportFromWkt导入Wkt中。Wkt获得最简便的方式是安装一个PostgreSQL，再安装PostGIS。在Windows系统中，有安装PostGIS的选项，而Linux系统需要自己编译。PostGIS的spatial_reference表中包含全套的空间参考，几乎所有的投影坐标信息的Wkt表达全在里面。按需拷贝，开箱即用。

假设已有Wkt，创建一个空间参考：
from osgeo import osr

wkt = spatialRef.ExportToWkt()
spatial = osr.SpatialReference()
spatial.ImportFromWkt(wkt)

最后用spatial替代CreateLayer的第二个参数None，完成建立空间参考的矢量数据投影

4.对矢量数据进行投影转换

可对矢量数据要素的几何特征进行投影转换，首先先将两个Projection初始化，再创建一个CoordinateTransformation对象：

targetSR = osr.SpatialReference()
targetSR.ImportFromEPSG(4326) #Geo WGS84
coordTrans = osr.CoordinateTransformation(spatial, targetSR)
feature = layer.GetFeature(0)
geom = feature.GetGeometryRef()
geom.ExportToWkt()
geom.Transform(coordTrans)

另外还需注意：要在适当的时候编辑Geometry，投影变换后最好不要再修改Geometry。一个数据源（DataSource）里面的所有Geometry都得作投影变换，可以用循环语句。

5.将投影写入文件

若单独处理投影信息，将投影信息单独写入文件时，需将你的投影写入.prj文件，用MorphToESRI()转成字符串，然后新建文本文件，写入文件即可：

targetSR.MorphToESRI()
file = open('/tmp/test.prj', 'w')
file.write(targetSR.ExportToWkt())
file.close()