Diferença em geometrias de shapefiles

Aug 21 2020

Eu tenho a geometria da região de contorno como parte de um shapefile (que é apenas um recurso e há outros recursos além deste). Os limites de outras regiões poligonais que são coloridas estão presentes em outro shapefile.

Existe uma maneira de descobrir a geometria de regiões não coloridas (como um polígono ou vários polígonos) da imagem?

Tentei com o código abaixo, mas ele estava me dando a geometria de toda a região.

from shapely.geometry import shape, mapping, Polygon
import matplotlib.pyplot as plt
import geopandas as gpd
import fiona

schema = {'geometry': 'Polygon','properties': {'test': 'float'}}

outline_shape = fiona.open(shapefile1)
region_shape = fiona.open(shapefile2)

for feature in outline_shape:
    if feature['properties']['Name'] == 'Required field':
        schema = {'geometry': 'Polygon', 'properties': {'test': 'int'}}
        with fiona.open('diff.shp', 'w', 'ESRI Shapefile', schema) as e:
            for geom in [shape(feature['geometry']).difference(shape(j['geometry'])) for j in region_shape]:
                if not geom.is_empty:
                    e.write({'geometry': mapping(geom), 'properties': {'test': 1}})

Respostas

1 bugmenot123 Aug 22 2020 at 23:45

Isso não foi testado, mas deve funcionar.

A ideia é fazer um polígono a partir do contorno para que você possa realmente calcular as diferenças de área. Para as regiões, fiz uma única geometria para que o código não tenha que recalcular as diferenças repetidamente, isso deve ser mais rápido.

Eu também mudei seu tratamento de arquivo para melhores práticas com withcontextos. E acrescentei algumas afirmações para pegar mal-entendidos.

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


# load all regions and make a single geometry of them
with fiona.open(regions_file) as regions:
    regions_geometries = [shape(f["geometry"]) for f in regions]
    regions_geometry = unary_union(regions_geometries)

# load the wanted outline
# via https://gis.stackexchange.com/questions/91676/select-by-attributes-within-the-fiona-python-module
with fiona.open(outlines_file) as outlines:
    filtered = filter(lambda f: f["properties"]["Name"] == "Required field", outlines)
    assert len(filtered) == 1, "more than 1 match?"
    outline = filtered[0]

outline_geometry = shape(outline["geometry"])

# assuming you actually have an outLINE (as LineString)
outline_geometry = Polygon(outline_geometry)

# compute the difference
difference_geometry = outline_geometry.difference(regions_geometry)
assert not difference_geometry.is_empty, "no remaining difference geometry?"

# ready to write to output
schema = {"geometry": "Polygon", "properties": {"test": "int"}}

difference_feature = {
    "geometry": mapping(difference_geometry),
    "properties": {"test": 1},
}

with fiona.open("diff.shp", "w", "ESRI Shapefile", schema) as output:
    output.write(difference_feature)

Pode estar faltando um CRS para a saída, mas isso é com você :)