2

我有两个包含多边形的形状文件。我试图从中找到增量。尝试通过以下代码执行此操作，但未按我预期的方式工作。以下是两个形状文件蓝色一个是缓冲区形状文件，我需要删除与蓝色缓冲区相交的缓冲区。即需要获得与Qgis 差异 函数相同的几何差异

在此处输入图像描述

import fiona
from shapely.geometry import shape, mapping, Polygon

green = fiona.open(
    "/home/gulve/manual_geo_ingestion/probe-data/images/r/shape_out/dissolved.shp")
blue = fiona.open(
    "/home/gulve/manual_geo_ingestion/probe-data/images/g/shape/shape.shp")

print([not shape(i['geometry']).difference(shape(j['geometry'])).is_empty for i, j in zip(list(blue), list(green))])

schema = {'geometry': 'Polygon',
          'properties': {}}
crs = {'init': u'epsg:3857'}

with fiona.open(
        '/home/gulve/manual_geo_ingestion/probe-data/images/r/shape_out/diff.shp', 'w',
        driver='ESRI Shapefile', crs=crs, schema=schema
) as write_shape:
    for geom in [shape(i['geometry']).difference(shape(j['geometry'])) for i, j in zip(list(blue), list(green))]:
        if not geom.empty:
            write_shape.write({'geometry': mapping((shape(geom))), 'properties': {}})

预期输出： 在此处输入图像描述

4

2 回答 2

2

将shapefile 导入 PostgreSQL后，只需执行以下查询：

CREATE TABLE not_intersects AS 
SELECT * FROM shape 
WHERE id NOT IN (SELECT DISTINCT shape.id
         FROM buffer,shape 
         WHERE ST_Intersects(buffer.geom,shape.geom));

此查询将创建第三个表（称为此处not_intersects），其中包含两个表（形状文件）之间不相交的多边形。

在此处输入图像描述

  • 结果以黄色表示。
于 2018-07-05T13:11:32.720 回答
0

我可以用匀称的函数来解决这个问题

这是我的代码..

import fiona
from shapely.geometry import shape, mapping, Polygon
from shapely.ops import unary_union

buffered_shape = fiona.open(
    "dissolved.shp", 'r', encoding='UTF-8')
color_shape = fiona.open(
    "shape.shp", 'r', encoding='UTF-8')

print([not shape(i['geometry']).difference(shape(j['geometry'])).is_empty for i, j in
       zip(list(color_shape), list(buffered_shape))])

outmulti = []
for pol in color_shape:
    green = shape(pol['geometry'])
    for pol2 in buffered_shape:
        red = shape(pol2['geometry'])
        if red.intersects(green):
            # If they intersect, create a new polygon that is
            # essentially pol minus the intersection
            intersect = green.intersection(red)
            nonoverlap = green.symmetric_difference(intersect)
            outmulti.append(nonoverlap)
        else:
            outmulti.append(green)

finalpol = unary_union(outmulti)

schema = {'geometry': 'MultiPolygon',
          'properties': {}}
crs = {'init': u'epsg:4326'}

with fiona.open(
        'shape_out/diff.shp', 'w',
        driver='ESRI Shapefile', crs=crs, schema=schema
) as write_shape:
    for geom in finalpol:
        # if not geom.empty:
        write_shape.write({'geometry': mapping(Polygon(shape(geom))), 'properties': {}})
于 2018-07-09T04:16:32.943 回答