[Python] GDAL/OGR操作矢量数据（shp、GeoJSON）

逃离我推掉我的手 2022-03-16 02:52 2750阅读 0赞

> GDAL项目旨于地理数据抽象模型对地理数据文件进行读写管理；而其项目下有两大类模块：GDAL和OGR  
> OGR提供操作矢量数据的API，GDAL模块提供栅格数据的API

【相关链接】

1.[GDAL/OGR例子合集][GDAL_OGR]  
2.[GDAL官网API][GDAL_API]  
3.[网上的学习笔记，教程来自于犹他州立大学][Link 1]  
4.[中国OSGeo教程][OSGeo]

### 文章目录 ###

*  OGR中的数据结构
 *   *  driver解析数据的驱动
     *  datasource管理数据源
     *  layer管理图层
     *  feature管理要素类
     *   *   *   *  从Layer中获得feature
                 *  更新要素
                 *  删除要素
                 *  关闭要素
                 *  创建要素
     *  geometry管理要素的几何属性
     *   *   *   *  Geometry的几何属性
                 *  从Feature获取Geometry
                 *  常用方法
                 *  创建Geometry
     *  field管理要素的属性字段
     *   *   *   *  创建字段
                 *  更改属性定义
                 *  其他
                 *  经验
     *  wkt
 *  查询
 *   *  属性查询
     *  空间查询
     *  SQL查询
 *  其他
 *   *  设置编码
 *  示例
 *   *  打开shp
     *  删除数据
     *  遍历所有属性值
     *  创建
     *   *   *   *  创建点状数据集
                 *  创建线状数据集
                 *  创建多边形数据集
                 *  创建属性
                 *  创建点Geometry
                 *  创建线Geometry
                 *  创建多边形Geometry
     *  拷贝
     *   *   *   *  datasource层次的拷贝
                 *  layer层次的拷贝
                 *  feature层次的拷贝

# OGR中的数据结构 #

![OGR中的数据结构][OGR]

## driver解析数据的驱动 ##

# 【方法一】从已有数据源中获取驱动变量
    ds = ogr.open(r'D:\....\...shp')
    driver = ds.GetDriver()
    # 【方法二】通过名称获取
    json_driver = ogr.GetDrverByName('GeoJSON')
        # 获得驱动名的两种方法
        # 1. OGR网站上有介绍，通过GDAL/ORG自带的ogrinfo
        # 2. 代码中提供的print_drivers函数来获取驱动程序的名字

## datasource管理数据源 ##

**（1）创建数据源**  
有了驱动程序之后，提供数据源名称来使用它创建一个空的数据源。新创建的数据源会自动打开等待写入

# 创建一个功能齐全的SpatialLite数据源，而不是使用SQLite
    fn = "filename"
    
    driver = ogr.GetDriverByName("SQLite")
    # 创建新的数据源时，不能覆盖现有的数据源。如果你的代码可能会覆盖现有数据，那么在创建新数据之前需要删除旧数据
    if os.path.exists(fn):
    	driver.DeleteDataSource(fn)
    ds= driver.CreateDataSource(fn, 
    	options=[ #创建选项：传递一个字符串列表（参数在OGR网站上有文档介绍）
    		'SPATIALITE=yes'
    	]
    )
    if ds is None:
    	#【注意】如果数据源没有创建成功，那么CreateDataSource函数会返回为空。如果为空对象，之后使用会报AttributeError错误
    	sys.exit('Could not create {0}'.formaat(json_fn ) )

**（2）删除数据源**

datasource.Destory() #关闭
    del datasource # 删除datasource
    	""" 删除ds变量强制关闭文件，并将所有的编辑成果写入磁盘中 【注意】删除图层变量并不会触发这个操作，必须关闭数据源才行 【备注】如果想保持数据源打开的话，可以通过图层对象或者数据源对象调用ds.SyncToDisk() 【警告】为了使你的编辑写入到磁盘中，必须关闭文件或者调用SyncToDisk函数。如果没有这么做，并且在交互环境中还打开数据源，那么你会很失望地发现创建了一个空的数据集 """

## layer管理图层 ##

layer = datasource.GetLayer(0) #根据下标来获取图层（对于shp而言只有一个图层）
    layer.GetGeomType() #图层几何类型
    n = layer.GetFeatureCount() #要素数量
    extent = layer.GetExtent() #上下左右边界
    readedNum = GetFeaturesRead() #已经读取多少条Feature
    layer.ResetReading() #重置内部feature指针，指向FID=0的要素
    layer.GetSpatialRef() #使用wkt文本表示空间参考系统的实例
    layer.schema #获取FieldDefn对象列表（属性名）

## feature管理要素类 ##

##### 从Layer中获得feature #####

# 1. 【获得图层中的要素】
    feature = layer.GetFeature(0) 
    # 2. 【获取图层中所有要素】
    feat = layer.GetNextFeature()
    while feat:
    	feat = layer.GetNextFeature()
    layer.ResetReading()

##### 更新要素 #####

> 要素更新操作也大同小异，只不过你操作的是图层中的现有要素，而不是一个空白的对象

# 【例子】为图层的每一个要素添加一个唯一的ID
    layer.CreateField( # 为图层添加一个ID字段
    	ogr.FieldDefn('ID', ogr.OFTInteger)
    ) 
    n = 1
    for feat in layer:
    	feat.SetField('ID', n) #设置字段的值
    	layer.SetFeature(feat) #将数值传递给SetFeature函数，来更新图层中的要素信息（这里不是CreateFeature函数）
    	n += 1

##### 删除要素 #####

> 删除要素更容易，需要知道的就是待删除要素的FID

#【例子】删除字段City_Name为"Seattle"的要素
    # 找到目标的FID进行删除
    for feat in layer:
    	if feat.GetField('City_Name')=='Seattle':
    		layer.DeleteFeature(feat.GetFID() ) #删除操作
    
    """ 【但是】 1. 某些特定的数据格式使用此方式并不能完全删除要素。 2. 你可能看不出，有时要素只是被标记为删除了，而不是被完全抛出，它仍然还潜伏着。 3. 正因如此，你不会看到其他要素被分配到刚才删除的FID 4. 这意味着，如果已经删除了很多要素，在文件中可能存在大量不必要的已用空间 【所以】删除这些要素可以回收对应的空间 【如何回收对应空间】如果你有一些关系数据库的经验，应该熟悉这一点。它类似于在微软的Access数据库上运行压缩和修复或在PostgreSQL数据库中使用重组功能（VACUUM） 【具体做法】打开数据源，然后执行一条SQL语句来压缩数据库（不同的数据格式回收的方法不同） """
    ds.ExecuteSQL('REPACK' + layer.GetName() ) #shapefile格式回收方式
    ds.ExecuteSQL('RECOMPUTE EXTENT ON' + layer.GetName() ) #确保空间范围更新
    	# 当现有要素发生改变或被删除后，它不会更新元数据中的空间范围

##### 关闭要素 #####

feature.Destory() #关闭

##### 创建要素 #####

> 往现有图层中添加新要素和往全新的图层中添加要素的操作一样，具体请看示例

【步骤】  
1.基于图层字段创建一个空要素，填充它  
2.然后把它插入到该图层

## geometry管理要素的几何属性 ##

##### Geometry的几何属性 #####

【org.Geometry中的几何类型】ogr.wkbPoint、org.wkbLineString、org.wkbPolygon等，请看wkt

# 获得Geometry的几何属性
    layer.GetGeomType() == ogr.wkbPolygon
    layer.GetGeometryName()

##### 从Feature获取Geometry #####

# 两种方法
    geometry = feature.geometry()
    geometry = feat.GetGeometryRef()

##### 常用方法 #####

geom.GetGeometryName() #要素类型
    geom.GetPointCount() #要素点个数
    geom.GetX(0);geom.GetY(0);geom.GetZ(0) #获得第一个坐标点的X、Y、Z
    print(geom) #打印出所有点
    geom.ExportToWkt() #导出WKT

##### 创建Geometry #####

请看示例

## field管理要素的属性字段 ##

【字段类型常量】

> 有一些不同的属性字段类型，但并不是每一种字段类型都支持所有的数据格式。这时候用于描述各种数据格式的在线文档就派上用场了。

<table> 
 <thead> 
  <tr> 
   <th>字段数据类型</th> 
   <th>OGR常量</th> 
  </tr> 
 </thead> 
 <tbody> 
  <tr> 
   <td>Integer</td> 
   <td>OFTInteger</td> 
  </tr> 
  <tr> 
   <td>List of integers</td> 
   <td>OFTIntegerList</td> 
  </tr> 
  <tr> 
   <td>Floating point number</td> 
   <td>OFTReal</td> 
  </tr> 
  <tr> 
   <td>List of floating point numbers</td> 
   <td>OFTRealList</td> 
  </tr> 
  <tr> 
   <td>String</td> 
   <td>OFTString</td> 
  </tr> 
  <tr> 
   <td>List of strings</td> 
   <td>OFTStringList</td> 
  </tr> 
  <tr> 
   <td>Date</td> 
   <td>OFTDate</td> 
  </tr> 
  <tr> 
   <td>Time of day</td> 
   <td>OFTTime</td> 
  </tr> 
  <tr> 
   <td>Date and time</td> 
   <td>OFTDateTime</td> 
  </tr> 
 </tbody> 
</table>

##### 创建字段 #####

在Layer中新建属性字段

# 创建并添加第一个属性字段
    coord_fld = ogr.FieldDefn('X', ogr.OFTReal)
    # 设置限制
    coord_fld.SetWidth(8) 
    coord_fld.SetPrecision(3)
    out_lyr.CreateField(coord_fld)
    # 重用FieldDefn对象来创建第二个字段
    coor_fld.SetName('Y')
    out_lyr.CreateField(coord_fld)

##### 更改属性定义 #####

【函数】`AlterFieldDefn(iField, field_def, nFlags)`  
【作用】用新字段的规则替换现有的字段  
【参数说明】

1.  iField是想改变的属性字段对应的索引值
2.  field\_def是新属性字段的定义对象
3.  nFlags是一个整数，是下表中的一个或多个常数的总和

【用于指明字段定义的哪个属性可以更改的标记，如果想使用多个，一起添加即可】

<table> 
 <thead> 
  <tr> 
   <th>需要更改的字段属性</th> 
   <th>OGR常量</th> 
  </tr> 
 </thead> 
 <tbody> 
  <tr> 
   <td>Field name only</td> 
   <td>ALTER_NAME_FLAG</td> 
  </tr> 
  <tr> 
   <td>Field type only</td> 
   <td>ALTER_TYPE_FLAG</td> 
  </tr> 
  <tr> 
   <td>Field width and/or precision only</td> 
   <td>ALTER_WIDTH_PRECISION_FLAG</td> 
  </tr> 
  <tr> 
   <td>All of the above</td> 
   <td>ALTER_ALL_FLAG</td> 
  </tr> 
 </tbody> 
</table>

# 【例子一】 更改原字段的名字
    index = layer.GetLayerDefn().GetFieldIndex('NAME') #获取字段的索引值
    fld_defn = ogr.FieldDefn('City_Name', ogr.OFTString) #创建新属性的字段定义
    layer.AlterFieldDefn(index, fld_defn, ogr.ALTER_NAME_FLAG) 
    
    # 【例子二】 更改字段的多个属性，例如名称和精度
    lyr_defn = layer.GetLayerDefn()
    index = lyr_defn.GetFieldIndex('X')
    width = lyr_defn.GetFieldDefn(index).GetWidth()
    fld_defn = ogr.FieldDefn('X_coord', ogr.OFTReal)
    fld_defn.SetWidth(width)
    fld_defn.SetPrecision(4)
    flag = ogr.ALTER_NAME_FLAG + ogr.ALTER_WIDTH_PRECISION_FLAG
    layer.AlterFieldDefn(index, fld_defn, flag)
    	"""【注意】在此例子中创建新的字段定义时使用的是原始字段宽度 如果没有为原始数据设置足够大的字段宽度，结果将不正确 为了解决这个问题，需要使用与原始数据一样的字段宽度 为了使更改的精度生效，所有记录必须重写 为数据设置一个比它自身更高的精度并不会提高精度，因为数据不能凭空产生，但如果设置的精度不够高，精度就会降低 """

##### 其他 #####

feature.GetFieldCount() # 属性总数
    feature.keys() # 属性名
    feat.GetField('AREA') # 获取字段
    feature.GetFieldAsString("FID") # 获取字段成string

##### 经验 #####

1.没能为GeoJSON文件设置一个精度  
2.想在shapefile中设置字段精度，必须设置字段宽度  
3.设置字段的默认值：`FieldDefn.SetDefault("我是默认值")`

## wkt ##

> WKB（二进制Well-KnownBinary），用于不同软件程序间进行几何要素类型转换的一种二进制表示标准。  
> 因为它是二进制格式，所以人们无法直接阅读获取其表示的内容，但是熟知文本格式（Well-Known Text，WKT）可以阅读。

<table> 
 <thead> 
  <tr> 
   <th>几何要素类型</th> 
   <th>对应OGR常量</th> 
  </tr> 
 </thead> 
 <tbody> 
  <tr> 
   <td>Point</td> 
   <td>wkbPoint</td> 
  </tr> 
  <tr> 
   <td>Multipoint</td> 
   <td>wkbMultiPoint</td> 
  </tr> 
  <tr> 
   <td>Line</td> 
   <td>wkbLineString</td> 
  </tr> 
  <tr> 
   <td>Multiline</td> 
   <td>wkbMultiLineString</td> 
  </tr> 
  <tr> 
   <td>Polygon</td> 
   <td>wkbPolygon</td> 
  </tr> 
  <tr> 
   <td>Multipolygon</td> 
   <td>wkbMultiPolygon</td> 
  </tr> 
  <tr> 
   <td>Unknown geometry type</td> 
   <td>wkbUnknown（图层如果有多种几何类型就返回wkbUnknown）</td> 
  </tr> 
  <tr> 
   <td>No geometry</td> 
   <td>wkbNone</td> 
  </tr> 
 </tbody> 
</table>

# 查询 #

## 属性查询 ##

from osgeo import ogr
    import os
    shpfile = r'C:\tmp\data.shp'
    ds = ogr.Open(shpfile)
    layer = ds.GetLayer(0) #得到图层
    lyr_count = layer.GetFeatureCount()
    print(lyr_count) #原先的要素总数
    layer.SetAttributeFilter("AREA > 800000") #属性查询条件
    lyr_count = layer.GetFeatureCount() #查询过后，layer被更新
    print(lyr_count) #满足条件的layer
    # 将该layer保存成shp
    driver = ogr.GetDriverByName("ESRI Shapefile")
    extfile = r'C:\tmp\data.shp'
    if os.access( extfile, os.F_OK ):
    	driver.DeleteDataSource( extfile )
    newds = driver.CreateDataSource(extfile)
    layernew = newds.CreateLayer('rect',None,ogr.wkbPolygon)
    # 遍历，复制
    feat = layer.GetNextFeature()
    while feat is not None:
    	layernew.CreateFeature(feat)
    	feat = layer.GetNextFeature()
    newds.Destroy()

## 空间查询 ##

## SQL查询 ##

from osgeo import ogr
    driver = ogr.GetDriverByName("ESRI Shapefile")
    world_shp = r'C:\tmp\test.shp'
    world_ds = ogr.Open(world_shp)
    world_layer = world_ds.GetLayer()
    world_layer_name = world_layer.GetName()
    result = world_ds.ExecuteSQL("select * from %s where prov_id = '22' order by BNDRY_ID desc" % (world_layer_name)) # ) # ExecuteSQL是基于数据集进行的，而不是图层
    resultFeat = result.GetNextFeature ()
    out_shp = r'C:\tmp\test\test_sql_result.shp'
    create_shp_by_layer(out_shp, result) #保存结果
    # 对查询结果进行遍历
    while resultFeat :
        print resultFeat.GetField('BNDRY_ID')
        resultFeat = result.GetNextFeature ()
    #执行下一条SQL语句之前一定要先释放
    world_ds.ReleaseResultSet(result)

# 其他 #

## 设置编码 ##

# 全局设置
    gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES") #支持文件名称及路径内的中文
    gdal.SetConfigOption("SHAPE_ENCODING", "GBK")        #支持属性字段中的中文
    # 为图层指定编码
    out_layer = ds.CreateLayer( "CommercialHousing",
        srs = in_layer.GetSpatialRef(),
        geom_type= ogr.wkbPoint,
        options  = [ 
            "ENCODING=UTF-8"
        ]
    )

# 示例 #

## 打开shp ##

from osgeo import ogr
    inshp_path = 'D:\mycode\GISandPython\0data\park_point_shp\xiamen_20181128_park.shp'
    driver = ogr.GetDriverByName('ESRI Shapefile') #查找一个特定的驱动程序
    datasource = driver.Open(inshp_path, 0) #0只读，1可写
    dir( datasource) #使用Python的内省函数dir()查看所有方法

## 删除数据 ##

# 删除矢量数据文件
    if os.path.exists(out_shp)
    	driver.DeleteDataSource(out_shp)
    ds2 = driver.CreateDataSource(out_shp)

## 遍历所有属性值 ##

#【遍历所有属性值】
    for i in range(feature.GetFieldCount() ):
    	print( feature.GetField(i) )
    ### 【查看表的结构，各个字段的名称等信息】在layer附加信息中看
    layerdef = layer.GetLayerDefn()
    for i in range(layerdef.GetFieldCount() ):
    	defn = layerdef.GetFieldDefn(i)
    	print(defn.GetName(), defn.GetWidth(), defn.GetType(), defn.GetPrecision() )

## 创建 ##

##### 创建点状数据集 #####

from osgeo import ogr
    import os,math
    driver = ogr.GetDriverByName("ESRI Shapefile")
    extfile = 'point_demo.shp'
    point_coors = [[300,450], [750, 700], [1200, 450], [750, 200], [750, 450]]
    print point_coors
    driver = ogr.GetDriverByName("ESRI Shapefile")
    if os.access( extfile, os.F_OK ):
        driver.DeleteDataSource( extfile )
    newds  = driver.CreateDataSource(extfile)
    layernew = newds.CreateLayer('point',None,ogr.wkbPoint)
    for aa in point_coors:
        wkt = 'POINT (' + str(aa[0]) + ' ' + str(aa[1]) + ')'
        geom = ogr.CreateGeometryFromWkt(wkt)
        feat = ogr.Feature(layernew.GetLayerDefn())
        feat.SetGeometry(geom)
        layernew.CreateFeature(feat)
    newds.Destroy()

##### 创建线状数据集 #####

from osgeo import ogr
    import os,math
    driver = ogr.GetDriverByName("ESRI Shapefile")
    extfile = 'line_demo.shp'
    point_coors = [300,450, 750, 700, 1200, 450, 750, 200]
    print point_coors
    driver = ogr.GetDriverByName("ESRI Shapefile")
    if os.access( extfile, os.F_OK ):
        driver.DeleteDataSource( extfile )
    newds  = driver.CreateDataSource(extfile)
    layernew = newds.CreateLayer('point',None,ogr.wkbLineString)
    wkt = 'LINESTRING (%f %f,%f %f,%f %f,%f %f,%f %f)' % (point_coors[0],point_coors[1],
        point_coors[2],point_coors[3], point_coors[4],point_coors[5],
        point_coors[6],point_coors[7], point_coors[0],point_coors[1])
    print(wkt)
    geom = ogr.CreateGeometryFromWkt(wkt)
    feat = ogr.Feature(layernew.GetLayerDefn())
    feat.SetGeometry(geom)
    layernew.CreateFeature(feat)
    newds.Destroy()

##### 创建多边形数据集 #####

from osgeo import ogr
    import os,math
    driver = ogr.GetDriverByName("ESRI Shapefile")
    extfile = 'rect_demo.shp'
    if os.access( extfile, os.F_OK ):
     driver.DeleteDataSource( extfile )
    extent = [400, 1100, 300, 600]
    newds = driver.CreateDataSource(extfile)
    layernew = newds.CreateLayer('rect',None,ogr.wkbPolygon)
    width = math.fabs(extent[1]-extent[0])
    height = math.fabs(extent[3]-extent[2])
    tw = width/2
    th = width/2
    extnew = extent[0]+tw
    wkt = 'POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f))' % (extent[0],extent[3],
        extent[1],extent[3], extent[1],extent[2],
        extent[0],extent[2], extent[0],extent[3])
    geom = ogr.CreateGeometryFromWkt(wkt)
    feat = ogr.Feature(layernew.GetLayerDefn())
    feat.SetGeometry(geom)
    layernew.CreateFeature(feat)
    newds.Destroy()

##### 创建属性 #####

from osgeo import ogr
    import os,math
    driver = ogr.GetDriverByName("ESRI Shapefile") #shp驱动器
    extfile = 'rect_field_demo.shp' #输出的shp路径
    if os.access( extfile, os.F_OK ): #文件是否已存在
        driver.DeleteDataSource( extfile )
    
    extent = [400, 1100, 300, 600] #范围
    newds = driver.CreateDataSource(extfile) #按照驱动器创建数据源
    layernew = newds.CreateLayer('rect',None,ogr.wkbPolygon) #创建图层
    # 创建字段
    fieldcnstr = ogr.FieldDefn("fd",ogr.OFTString) #创建字段(字段名，类型)
    fieldcnstr.SetWidth(32) #设置宽度
    layernew.CreateField(fieldcnstr) #将字段设置到layer
    fieldf = ogr.FieldDefn("f",ogr.OFTReal) 
    layernew.CreateField(fieldf)
    # wkt格式的polygon
    wkt = 'POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f))' % (extent[0],extent[3],
        extent[1],extent[3], extent[1],extent[2],
        extent[0],extent[2], extent[0],extent[3])
    geom = ogr.CreateGeometryFromWkt(wkt) #根据wkt创建geometry
    feat = ogr.Feature( layernew.GetLayerDefn() ) #从layer中得到feature的要求，创建一个Feature对象
    feat.SetField('fd',"这里是字段的值") #设置字段值
    feat.SetGeometry(geom) #设置要素的geometry属性
    layernew.CreateFeature(feat) #在图层中创建该feature
    newds.Destroy() #删除（包括保存）

##### 创建点Geometry #####

from osgeo import ogr
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(10,20) #向Point添加多个点，不会报错，但最终只会用最后添加的点

##### 创建线Geometry #####

from osgeo import ogr
    line = ogr.Geometry(ogr.wkbLineString) # 创建
    line.AddPoint(10, 20) # 添加点
    line.SetPoint(1, 50, 50) # 修改点坐标
    line.GetPointCount() # 得到所有点数目
    line.GetX(0); line.GetY(0) # 读取0号点的坐标

##### 创建多边形Geometry #####

from osgeo ipmort ogr
    # 创建环
    ring = ogr.Geometry(ogr.wkbLinearRing) #创建
    ring.AddPoint(10,10) #添加点
    ring.CloseRings() #闭合 或 添加的最后一个点和第一个相同
    # 将环添加到多边形中
    polygon = ogr.Geometry(ogr.wkbPolygon)
    polygon.AddGeometry(outring)
    polygon.AddGeometry(inring)
    ring = polygon.GetGeometryRef(0) #得到多边形内的几何对象

## 拷贝 ##

##### datasource层次的拷贝 #####

from osgeo import ogr
    import os,math
    inshp = "C:\tmp\test.shp"
    ds = ogr.Open(inshp)
    driver = ogr.GetDriverByName("ESRI Shapefile")
    outputfile = 'C:\tmp\test_copy.shp'
    if os.access( outputfile, os.F_OK) : #已经有该文件
    	driver.DeleteDataSource(outputfile)  #删除
    pt_cp = driver.CopyDataSource(ds, outputfile) #根据数据源ds创建数据，返回指针
    pt_cp.Release() #释放指针，将数据写入到磁盘

##### layer层次的拷贝 #####

from osgeo import ogr
    import os,math
    inshp = "C:\tmp\test.shp"
    ds = ogr.Open(inshp)
    driver = ogr.GetDriverByName("ESRI Shapefile")
    outputfile = "C:\tmp\test_copy2.shp"
    if os.access( outputfile, os.F_OK):
    	driver.DeleteDataSourse( outputfile )
    layer = ds.GetLayer()
    newds = driver.CreateDataSource(outputfile) #首先得有数据源
    pt_layer  = newds.CopyLayer(layer, 'abcd') #复制图层，返回指针
    newds.Destroy() #对newds进行Destroy()操作，才能将数据写入磁盘

##### feature层次的拷贝 #####

from osgeo import ogr
    import os,math
    inshp = 'C:\tmp\test.shp'
    ds = ogr.Open(inshp)
    driver = ogr.GetDriverByName("ESRI Shapefile")
    outputfile = 'c:\tmp\test.shp'
    if os.access( outputfile, os.F_OK ):
        driver.DeleteDataSource( outputfile )
    newds = driver.CreateDataSource(outputfile) #按outputfile的数据源创建一个newds
    layernew = newds.CreateLayer('worldcopy',None,ogr.wkbLineString) #创建一个图层
    layer = ds.GetLayer() #得到原数据源的图层
    extent = layer.GetExtent() #图层范围
    # 遍历旧图层的所有feature，进行复制
    feature = layer.GetNextFeature()
    while feature is not None: 
        layernew.CreateFeature(feature)
        feature = layer.GetNextFeature()
    newds.Destroy() # 包括了将数据flush到磁盘

[GDAL_OGR]: http://pcjericks.github.io/py-gdalogr-cookbook/#
[GDAL_API]: https://gdal.org/python/index.html
[Link 1]: https://www.osgeo.cn/python_gdal_utah_tutorial/
[OSGeo]: https://www.osgeo.cn/pygis/ogr-intro.html
[OGR]: /images/20220316/72d88d010fab4a9fa6a58521b8d61ec4.png