Diferença em geometrias de shapefiles
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
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ê :)