矢量裁剪矢量

　　近日在博客中看到一篇使用ogr裁剪矢量数据的文章，觉得挺好，就做个笔记来学习。文章链接https://www.cnblogs.com/ljwgis/p/14071604.html

核心代码

from osgeo import ogr, osr
import os

def createShpByLayer(shp, layer, fileType):
    '''
    根据layer创建shpfile
    '''
    driver = ogr.GetDriverByName("ESRI Shapefile")
    ds = driver.CreateDataSource(shp)
    pt_layer = ds.CopyLayer(layer, 'layername')
    ds.Destroy()
    
def splitShp(shpfile, outPath, splitField):
    '''
    按属性字段分割shpfile
    '''
    ds = ogr.Open(shpfile)
    lyr = ds.GetLayer(0)
    for feat in lyr:
        cityName = feat.GetField(splitField)  # 以字段名为文件名称
        outShp = os.path.join(outPath, str(cityName)+'.shp')
        geom = feat.GetGeometryRef()
        driver = ogr.GetDriverByName("ESRI Shapefile")
        outDs = driver.CreateDataSource(outShp)
        outLyr = outDs.CreateLayer("layername", lyr.GetSpatialRef(), ogr.wkbMultiPolygon)
        outLyr.CreateFields(lyr.schema)  # 创建字段属性
        outFeat = ogr.Feature(lyr.GetLayerDefn())
        for i in range(feat.GetFieldCount()):
            val = feat.GetField(i)
            outFeat.SetField(i, val)
        outFeat.SetGeometry(geom)
        outLyr.CreateFeature(outFeat)
        outDs = None
        
def splitShpByShp(covershp, splitshp, outPath):
    '''
    @desc: 将一个shp按属性分割后，再用来分割另一个shp（空间查询）
    '''
    ds = ogr.Open(covershp)
    lyr = ds.GetLayer(0)
    splitDs = ogr.Open(splitshp)
    splitLyr = splitDs.GetLayer(0)
    for feat in lyr:
        cityName = feat.GetField("NAME")  # 以字段名为文件名称
        outShp = os.path.join(outPath, str(cityName)+'.shp')
        geom = feat.GetGeometryRef()
        driver = ogr.GetDriverByName("memory")
        outDs = driver.CreateDataSource("temp")
        outLyr = outDs.CreateLayer("layername", lyr.GetSpatialRef(), ogr.wkbMultiPolygon)
        outLyr.CreateFields(lyr.schema)  # 创建字段属性
        outFeat = ogr.Feature(lyr.GetLayerDefn())
        for i in range(feat.GetFieldCount()):
            val = feat.GetField(i)
            # val = val.encode('utf8') # 将unicode编码转化为中文，还有点问题
            outFeat.SetField(i, val)
        outFeat.SetGeometry(geom)
        outLyr.CreateFeature(outFeat)
        outDs = None
        
        splitLyr.SetSpatialFilter(geom) # 按空间查询
        createShpByLayer(outShp, splitLyr, "ESRI Shapefile")
        print("按空间查询的要素数量："+ str(splitLyr.GetFeatureCount()))