Python空间数据处理—-矢量数据在不同空间参考下的转换

1. 流程图如下：
   

2. 需要用到的函数：
   （1）spatialref=osr.SpatialReference()
   该函数是SpatialReference类的一个构造函数，它构造一个 SpatialReference对象，另外，它接受wkt格式的字符串来构建SpatialReference对象，如果没有接受，可以通过SpatialReference对象调用成员函数ImportFromWkt()、ImportFromEPSG()来构建。
   （2）spatialref.ImportFromEPSG(4326)
   该函数从EPSG中导入一个空间参考。
   （3）trans=osr.CoordinateTransformation(spatialref_source,spatialref_target)
   该函数是CoordinateTransformation（坐标系统转换）类的构造函数，用来创建坐标系统转换对象，后续可以通过该函数创建的坐标系统转换对象来对数据的空间参考进行转换。
   （4）driver=ogr.GetDriverByName(“ESRI Shapefile”)
   注册驱动，可以理解为：这一个函数的作用是告诉程序需要在磁盘中创建多少个格式的文件，如：对于 .shp文件，这一函数则告诉程序需要在文件中创建 .shp、.shx、.prj、.dbf等文件，创建完之后，可以将数据写入相应的文件中。该函数返回一个Driver类的对象。
   （5）datasource=driver.CreateDataSource(fn_path)
   创建数据源，矢量的组织形式如下图：
   从上图可以看出，若要将矢量数据输出，需要先创建一个DataSource，然后再在DataSource中添加Layer（图层），并且创建要素，将空间数据（Geometry）以及属性添加到要素中，最后在图层中添加要素。该函数返回一个DataSource类的对象。
   （6）layer_out=datasource.CreateLayer(“layer_out”,srs=spatialref_target,geom_type=ogr.wkbPolygon)
   该函数在DataSource中添加一个矢量图层，图层的名字为Layer_out，关键字参数srs=spatialref_target表示图层的空间参考为spatialref_targe，而geom_type=ogr.wkbPolygon表示图层中要素的集合类型为wkbPolygon。该函数返回一个Layer类的对象。
   （7）datasource_in=ogr.Open(fn_in_path,False)
   该函数打开一个矢量图层，第一个形式参数是数据的路径，第二个参数表示以只读的方式打开（默认也是只读方式）。该函数返回Datasource类的对象。如果返回None，则表示打开失败。
   （8）layer_in=datasource_in.GetLayer(0)
   该函数是DataSource类的一个成员函数，用来获取第一个图层，返回的是一个Layer类的对象。若返回None，表示获取第一个图层失败。另外，矢量图层的计数是从0开始的，而栅格数据则从1开始，其实shp文件只有一个图层。
   （9）layer_in_defn=layer_in.GetLayerDefn()
   该函数是Layer类的一个成员函数，作用是获取图层的属性表信息（或者说，获取图层的属性表表头信息），返回一个FeatureDefn类的对象。
   
   （10）field_count=layer_in_defn.GetFieldCount()
   该函数是FeatureDefn类的成员函数，可以获得属性列的列数。返回值是int型的整数，为属性列列数。
   （11）field_defn=layer_in_defn.GetFieldDefn(i)
   该函数是FeatureDefn类的成员函数，用来获取第i个字段信息（如：第i列的名称、数据类型等）。该函数返回FieldDefn类的对象。
   
   在该函数的返回对象中，可以获取上图中的Name、Precision、Type等信息。
   （12）layer_out.CreateField(field_defn)
   该函数是Layer类的一个对象，用来在图层中创建一个字段，其形式参数是FieldDefn类对象，若创建失败，则返回None。
   （13）feature_in=layer_in.GetNextFeature()
   该函数是Layer类的一个对象，作用是获取图层中的下一个要素，返回值为Feature类对象。
   （14）feature_out=ogr.Feature(layer_out_defn)
   该函数是Feature类的构造函数，构造一个Feature对象（构造一个要素），其形式参数是FeatureDefn类的对象（输出图层的属性表表头信息），返回值是Feature类对象。构造完Feature之后，需要创建Geometry以及字段信息并将二者塞入创建的Feature中。
   （15）geomref_in=feature_in.GetGeometryRef()
   该函数是Feature类的成员函数，作用是获取要素的空间坐标信息，返回值为Geometry类对象。
   （16）geomref_in.Transform(trans)
   该函数是Geometry类的成员函数，其形参为CoordinateTransform类的对象，用来将矢量数据从一种空间参考转为另一种空间参考。
   （17）feature_out.SetGeometry(geomref_in)
   该函数是Feature类的成员函数，其形参为Geometry类对象，用来将构成Feature（要素）的空间信息塞入Feature（要素）中。
   （18）feature_out.SetField(name,field_value)
   该函数是Feature类的成员函数，用来将构成Feature（要素）的属性信息塞入Feature（要素）中（将列名为name的值塞入Feature（要素）中）。
   （19）layer_out.CreateFeature(feature_out)
   该函数是Layer类的成员函数，形参为Feature类对象，用来将构建好的Feature（要素）写入图层中。

3. 代码如下：

#coding=utf-8
from osgeo import osr,ogr
#定义原始空间参考信息
spatialref_source=osr.SpatialReference()
spatialref_source.ImportFromEPSG(4326)
#定义目标空间参考信息
spatialref_target=osr.SpatialReference()
spatialref_target.ImportFromEPSG(3857)
#创建坐标转换对象
trans=osr.CoordinateTransformation(spatialref_source,spatialref_target)

#创建输出矢量数据
fn_out=r"F:\envi\text\shp\shp\shape_out.shp"
driver=ogr.GetDriverByName("ESRI Shapefile")
datasource_out=driver.CreateDataSource(fn_out)
layer_out=datasource_out.CreateLayer("layer_out",srs=spatialref_target,geom_type=ogr.wkbPolygon)


#读取矢量
fn_in=r"F:\envi\text\shp\shp\shape.shp"
datasource_in=ogr.Open(fn_in,False)
layer_in=datasource_in.GetLayer(0)
layer_in_defn=layer_in.GetLayerDefn()
field_count=layer_in_defn.GetFieldCount()


#给输出图层定义一个属性表,并创建字段
for i in range(field_count):
    field_defn=layer_in_defn.GetFieldDefn(i)
    layer_out.CreateField(field_defn)  #创建字段


layer_out_defn=layer_out.GetLayerDefn()



#从输入图层中取出要素，并进行坐标转换，以及创建要素以写入输出图层中
feature_in=layer_in.GetNextFeature()
i=0
while feature_in:
    feature_out=ogr.Feature(layer_out_defn)
    geomref_in=feature_in.GetGeometryRef()
    geomref_in.Transform(trans)
    feature_out.SetGeometry(geomref_in)
    for j in range(field_count):
        field_defn=layer_in_defn.GetFieldDefn(j)
        name=field_defn.name
        field_value=feature_in.GetField(name)
        feature_out.SetField(name,field_value)
    layer_out.CreateFeature(feature_out)
    feature_in=layer_in.GetNextFeature()
    
print "Finished!"
datasource_in=None
datasource_out=None

参考：http://pcjericks.github.io/py-gdalogr-cookbook/projection.html#reproject-a-layer

转载自：https://blog.csdn.net/qq_36181130/article/details/83543544