Differenza nelle geometrie degli shapefile

Question

Ho la geometria della regione del contorno come parte di uno shapefile (che è solo una caratteristica e ci sono altre caratteristiche oltre a questa). I confini di altre regioni del poligono che sono colorati sono presenti in un altro shapefile.

C'è un modo per scoprire la geometria delle regioni non colorate (complessivamente come un poligono o più poligoni) dell'immagine?

Ho provato con il codice qui sotto ma mi dava la geometria dell'intera regione.

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}})

Answer 1

Questo è completamente non testato ma dovrebbe funzionare.

L'idea è di creare un poligono dal contorno in modo da poter effettivamente calcolare le differenze di area. Per le regioni ho creato una singola geometria in modo che il codice non debba ricalcolare le differenze ancora e ancora, dovrebbe essere più veloce.

Ho anche modificato la gestione dei file in best practice con i withcontesti. E ho aggiunto alcune affermazioni per cogliere malintesi.

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)

Potresti perdere un SR per l'output, ma dipende da te :)