基于gdal的矢量交集（python）

基于gdal的矢量交集（python）

文章目录

​​代码​​​​数据结果可视化​​​​参考​​

代码

from pathlib import Pathimport ogrimport gdalimport osimport globimport shutildef get_outside(inShp, outputShp, field, fieldValue): driver = ogr.GetDriverByName("ESRI Shapefile") dataSource = driver.Open(inShp, 1) layer = dataSource.GetLayer() # 新建DataSource，Layer out_ds = driver.CreateDataSource(outputShp) out_lyr = out_ds.CreateLayer(outputShp, layer.GetSpatialRef(), ogr.wkbPolygon) def_feature = out_lyr.GetLayerDefn() for feature in layer: geom = feature.GetGeometryRef() value = feature.GetField(field) if value == fieldValue: out_feature = ogr.Feature(def_feature) out_feature.SetGeometry(geom) out_lyr.CreateFeature(out_feature) out_feature = None out_ds.FlushCache() del dataSource, out_dsdef intersection(borderShp, roadShp, fname): driver = ogr.GetDriverByName("ESRI Shapefile") dataSource = driver.Open(borderShp, 1) layer = dataSource.GetLayer() RdataSource = driver.Open(roadShp, 1) Rlayer = RdataSource.GetLayer() # 新建DataSource，Layer out_ds = driver.CreateDataSource(fname) out_lyr = out_ds.CreateLayer(fname, layer.GetSpatialRef(), ogr.wkbPolygon) def_feature = out_lyr.GetLayerDefn() # 遍历原始的Shapefile文件给每个Geometry做Buffer操作 # current_union = layer[0].Clone() print('the length of layer:', len(layer)) if len(layer) == 0: return for feature in layer: geometry = feature.GetGeometryRef() for Rfeature in Rlayer: Rgeometry = Rfeature.GetGeometryRef() inter = Rgeometry.Intersection(geometry).Clone() out_feature = ogr.Feature(def_feature) out_feature.SetGeometry(inter) out_lyr.ResetReading() out_lyr.CreateFeature(out_feature) del dataSource, RdataSource, out_dsdef mkdir(path): if not os.path.exists(path): os.mkdir(path)if __name__ == '__main__': f = open('config_order.txt') data = f.readlines() inShp = data[0].replace('\n', '') roadshp = data[1].replace('\n', '') outshp = data[2].replace('\n', '') field = 'gj' fieldValue = 0 mkdir('temp') tempShp = './temp/temp.shp' get_outside(inShp, tempShp, field, fieldValue) intersection(tempShp, roadshp, outshp) if os.path.exists('temp'): shutil.rmtree('temp')

数据结果可视化

参考

​​https://osgeo-/python_gdal_utah_tutorial/ch04.html#id1​​

版权声明：本文内容由网络用户投稿，版权归原作者所有，本站不拥有其著作权，亦不承担相应法律责任。如果您发现本站中有涉嫌抄袭或描述失实的内容，请联系我们jiasou666@gmail.com 处理，核实后本网站将在24小时内删除侵权内容。