geo-l2-nyc-park-symdiff

Rerun

import geopandas as gpd
import pandas as pd
from shapely.geometry import MultiPolygon, Polygon, Point
from shapely.ops import unary_union

def solve():
    # 1. Load data
    gdf_official = gpd.read_file("nyc_parks.gpkg", layer="parks_official")
    gdf_osm = gpd.read_file("nyc_parks.gpkg", layer="parks_osm")

    # Ensure CRS is consistent (though they should be from same GPKG)
    if gdf_official.crs != gdf_osm.crs:
        gdf_osm = gdf_osm.to_crs(gdf_official.crs)

    # 2. Compute Differences
    # "Claimed by exactly one side" means:
    # (Official - OSM) -> Official only
    # (OSM - Official) -> OSM only
    
    # Dissolve both to ensure we are dealing with unified geometries for difference ops
    union_official = gdf_official.unary_union
    union_osm = gdf_osm.unary_union

    diff_official = union_official.difference(union_osm)
    diff_osm = union_osm.difference(union_official)

    # 3. Handle MultiPolygon/Polygon fragmentation
    # We need to identify "adjacent disagreement patches merged together".
    # A "shifted park's gain and loss" should be a single feature.
    # Actually, if a park shifts, the area it left is (Official - OSM) 
    # and the area it gained is (OSM - Official).
    # These two areas are adjacent.
    
    # Let's build a single geometry representing all disagreements.
    # All disagreement = (Official XOR OSM)
    # But we need to group them.
    
    # Approach: 
    # Create a single geometry of all disagreements.
    # Then, find connected components of these disagreement areas.
    # Wait, the prompt says "every patch claimed by exactly one side, with adjacent 
    # disagreement patches merged together".
    
    # If a park shifts, the "loss" (Official \ OSM) and "gain" (OSM \ Official) 
    # share a boundary. 
    # We can merge all disagreement polygons into one set of polygons, 
    # then group them by adjacency.
    
    # Let's collect all "disagreement" polygons.
    all_disagreements = []
    
    # Convert diffs to a list of Polygons
    def get_polys(geom):
        if geom.is_empty:
            return []
        if isinstance(geom, Polygon):
            return [geom]
        if isinstance(geom, MultiPolygon):
            return list(geom.geoms)
        return []

    polys_official_only = get_polys(diff_official)
    polys_osm_only = get_polys(diff_osm)

    # We need to group them. If a poly from official_only touches a poly from osm_only, 
    # they belong to the same "disagreement group".
    
    # Combine all disjoint pieces into a list of geometries with their intended source.
    # However, the "source" attribute for a merged group must be one of:
    # 'parks_official', 'parks_osm', or 'both'.
    # Wait, the prompt says: "A source attribute whose value is exactly one of 
    # 'parks_official', 'parks_osm', or 'both'".
    # But it also says: "Every patch claimed by exactly one side".
    # If a group contains patches from BOTH sides (like a shifted park), 
    # the source should be 'both'.
    
    all_pieces = []
    for p in polys_official_only:
        all_pieces.append({'geometry': p, 'source': 'parks_official'})
    for p in polys_osm_only:
        all_pieces.append({'geometry': p, 'source': 'parks_osm'})

    if not all_pieces:
        print("No disagreements found.")
        return

    # Grouping by adjacency:
    # We can use a Disjoint Set Union (DSU) or a graph approach.
    # Two pieces are "adjacent" if they intersect or touch (share a boundary).
    # In Shapely, .touches() or .intersects()? 
    # If they share a boundary, they are adjacent. 
    # Since they are from a difference operation, they won't overlap, they will touch.
    
    num_pieces = len(all_pieces)
    adj = [[] for _ in range(num_pieces)]
    for i in range(num_pieces):
        for j in range(i + 1, num_pieces):
            if all_pieces[i]['geometry'].touches(all_pieces[j]['geometry']):
                adj[i].append(j)
                adj[j].append(i)

    # Find connected components
    visited = [False] * num_pieces
    groups = []
    for i in range(num_pieces):
        if not visited[i]:
            component = []
            stack = [i]
            visited[i] = True
            while stack:
                u = stack.pop()
                component.append(u)
                for v in adj[u]:
                    if not visited[v]:
                        visited[v] = True
                        stack.append(v)
            groups.append(component)

    # Create the merged features
    merged_features = []
    for group_indices in groups:
        group_geoms = [all_pieces[idx]['geometry'] for idx in group_indices]
        merged_geom = unary_union(group_geoms)
        
        # Determine source
        sources = set(all_pieces[idx]['source'] for idx in group_indices)
        if len(sources) > 1:
            final_source = 'both'
        else:
            final_source = list(sources)[0]
        
        # Label anchor: point strictly inside
        # representative_point() is guaranteed to be inside.
        anchor = merged_geom.representative_point()
        
        merged_features.append({
            'geometry': merged_geom,
            'source': final_source,
            'anchor': anchor
        })

    # Create GeoDataFrames
    # For parks_disagreement.geojson
    df_disag = gpd.GeoDataFrame([
        {'geometry': f['geometry'], 'source': f['source']} for f in merged_features
    ], crs=gdf_official.crs)
    df_disag.to_file("parks_disagreement.geojson", driver="GeoJSON")

    # For park_label_anchors.geojson
    df_anchors = gpd.GeoDataFrame([
        {'geometry': f['anchor'], 'source': f['source']} for f in merged_features
    ], crs=gdf_official.crs)
    df_anchors.to_file("park_label_anchors.geojson", driver="GeoJSON")
    
    print("Successfully created output files.")

if __name__ == "__main__":
    solve()