python 使用GDAL 面积相交查询

最新推荐文章于 2023-09-07 15:06:06 发布

ZGFshm

最新推荐文章于 2023-09-07 15:06:06 发布

阅读量497

点赞数

文章标签： python json

本文链接：https://blog.csdn.net/ZGFshm/article/details/130722263

版权

1、先试用GDAL导入shp 文件,

2、有A,B两个SHP范围 A shp相当于市的范围,B shp 相当于县的范围

,以B的shp的范围选取Ashp的范围选取出A的geometry

3，GDAL使用 Intersects进行相交查询

from osgeo import ogr
import os
from server.app.Common import CACHE_DIR
from osgeo import gdal
import json

#根据自己的python换将引入GDAL
os.environ['PROJ_LIB'] = r'C:\Users\Administrator\AppData\Local\Programs\Python\Python39\Lib\site-packages\osgeo\data\proj'
# 开启异常
gdal.UseExceptions()
shp_path = os.path.join(CACHE_DIR, 'tmp')
def duibiShp(multipolygons_json):
 

    # 打开 A 和 B 的 Shapefile 文件
    driver = ogr.GetDriverByName('ESRI Shapefile')
    A = driver.Open(os.getcwd()+'/lyshp/ly.shp', 0)
    B = driver.Open(multipolygons_json+'/out.shp', 0)



    # 获取 A 和 B 的 Layer
    A_layer = A.GetLayer()
    B_layer = B.GetLayer()

    # 创建一个空的输出图层
    driver = ogr.GetDriverByName('GeoJSON')
    fileOut = driver.CreateDataSource(os.getcwd()+'/geojsonOut/result.json')
    outLayer = fileOut.CreateLayer('result', B_layer.GetSpatialRef(), ogr.wkbMultiPolygon)
    for field in B_layer.schema:
        outLayer.CreateField(field)

    # 创建一个空间过滤器，用于查询 A 和 B 中相交大于 50% 的数据
    extent = B_layer.GetExtent()
    wktaa = '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])
    # spatial_filter = ogr.Geometry(ogr.wkbPolygon)
    spatial_filter=ogr.CreateGeometryFromWkt(wktaa)

    # 遍历 B 中的每个 Feature，查询与 A 中相交大于 50% 的 Feature，并打印查询结果
    for B_feature in B_layer:
        B_geometry = B_feature.GetGeometryRef()
        spatial_filter.Intersects(B_geometry)
        A_layer.SetSpatialFilter(spatial_filter)
        for A_feature in A_layer:
            A_geometry = A_feature.GetGeometryRef()
            intersection = A_geometry.Intersection(B_geometry)
            intersection_area = intersection.Area()
            # A_area = A_geometry.Area()
            A_area = B_geometry.Area()

            if intersection_area / A_area > 0.1:
                # print('B Feature ID: {}, A Feature ID: {}'.format(B_feature.GetFID(), A_feature.GetFID()))
                # print(intersection_area, '----------', A_area, intersection_area / A_area)
                # 添加 B 的几何到输出图层
                # 添加 B 的几何到输出图层
                outFeature = ogr.Feature(outLayer.GetLayerDefn())
                outFeature.SetGeometry(B_geometry)
                for i in range(B_feature.GetFieldCount()):
                    outFeature.SetField(outLayer.GetLayerDefn().GetFieldDefn(i).GetNameRef(),
                                        B_feature.GetField(i))
                # 添加输出图层
                outLayer.CreateFeature(outFeature)
    # result ={}
    # 读取输出的 GeoJSON 文件为 Python 对象
    # with open(os.getcwd()+'/geojsonOut/result.json') as f:
    #     result = json.load(f)
    # 保存输出文件
    outLayer = None
    fileOut.FlushCache()
    fileOut = None
    try:
        with open(os.getcwd()+'/geojsonOut/result.json') as f:
            data = json.load(f)
    except Exception as e:
        print('Error:', e)
    finally:
        f.close()
    # 删除输出的 GeoJSON 文件
    driver.DeleteDataSource(os.getcwd()+'/geojsonOut/result.json')
    print(data)

    # 关闭 Shapefile 文件
    A = None
    B = None
    return data