使用ogr裁剪矢量数据
使用ogr裁剪矢量数据
由来:
近期有个需求,内容是这样的:我们有两个矢量数据,现在要求以一个矢量文件为底板,按字段对另一个矢量文件进行分割,生成若干小的shpfile文件
分析:
经过分析之后,将步骤拆解如下:
-
首先确保两个shpfile投影坐标系统一
如果出现不统一的情况,那么用Arcgis的工具Project进行投影转换。(
Data Management Tools--Projections and Transformations--Project
),链接如下https://www.sohu.com/a/292487012_488161。 -
其次编写按属性字段分割shpfile,生成若干小shpfile的代码
-
接着编写根据layer创建shpfile的代码(https://www.cnblogs.com/bobird/articles/3079523.html)-->使用ogr中拷贝方法创建新的shpfile
-
最后进行循环,对每次生成的shpfile和原先的shpfile进行空间查询,得到layer后生成需要的shpfile
核心代码:
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()))
不足之处
- 转出来的shpfile属性字段中,中文部分乱码,目前未能解决。