シェープファイルのジオメトリの違い

Question

アウトライン領域のジオメトリを1つのシェープファイルの一部として持っています（これは1つの機能であり、これとは別に他の機能があります）。色付きの他のポリゴン領域の境界は、別のシェープファイルに存在します。

画像の色のない領域（1つのポリゴンまたは複数のポリゴンとして）のジオメトリを見つける方法はありますか？

以下のコードで試してみましたが、領域全体のジオメトリが得られました。

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

これは完全にテストされていませんが、うまくいくはずです。

アイデアは、実際に面積の違いを計算できるように、アウトラインからポリゴンを作成することです。コードが違いを何度も再計算する必要がないように、単一のジオメトリを作成した領域の場合、これはより高速になるはずです。

また、ファイル処理をwithコンテキストを使用したベストプラクティスに変更しました。そして、誤解を招くためにいくつかのアサーションを追加しました。

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)

出力のCRSが欠落している可能性がありますが、それはあなた次第です:)