python gdal使用wkt格式面裁剪shp

ylfmsn

已于 2024-03-25 14:35:47 修改

阅读量121

点赞数 2

分类专栏： python gdal 文章标签： python 开发语言

于 2024-03-25 14:34:59 首次发布

本文链接：https://blog.csdn.net/ylfmsn/article/details/137013227

版权

python 同时被 2 个专栏收录

12 篇文章 3 订阅

订阅专栏

gdal

2 篇文章 0 订阅

订阅专栏

利用wkt对shp矢量图层进行裁剪

from osgeo import ogr
from osgeo.ogr import (
	Driver,
	DataSource,
    Layer,
    FeatureDefn,
    Feature,
    Geometry
)

wkt = 'MULTIPOLYGON (((10959600.1600258 4491987.15336097,11663148.6944404 4780944.58713839,12448359.1123138 4598775.77019176,12517457.6290866 4203029.71958356,12555147.7291446 3832410.4023473,12203373.4619373 3738185.15220249,11732247.2112132 4133931.20281069,11361627.893977 3782156.9356034,10815121.4431371 3637678.21871469,10815121.4431371 3637678.21871469,10959600.1600258 4491987.15336097)))'
shp_file = '/data/test/province.shp'
output_file = '/data/test/clip.shp'

# wkt创建geometry，要求wkt数据坐标系要与shp坐标系一致
clip_polygon: Geometry = ogr.CreateGeometryFromWkt(wkt)

# 打开shp数据
driver: Driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource: DataSource = driver.Open(shp_file, 0)
layer: Layer = dataSource.GetLayer()
layerDefinition: FeatureDefn = layer.GetLayerDefn()

# 创建目标图层
target_ds: DataSource = driver.CreateDataSource(output_file)
target_layer: Layer = target_ds.CreateLayer("clip", layer.GetSpatialRef(), ogr.wkbUnknown, options=["ENCODING=GBK"])
for i in range(layerDefinition.GetFieldCount()):
    target_layer.CreateField(layerDefinition.GetFieldDefn(i))

# 先过滤，筛选出相交部分
layer.SetSpatialFilter(clip_polygon)

# 对相交部分求相切
for feature in layer:
    feature: Feature = feature
    geom: Geometry = feature.GetGeometryRef()
    intersection = geom.Intersection(clip_polygon)
    feature.SetGeometry(intersection)
    target_layer.CreateFeature(feature.Clone())
    feature.Destroy()
# 关闭数据源
dataSource.Destroy()
target_ds.Destroy()