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 get_representative_point(geom):
    """Get a point strictly inside the geometry."""
    return geom.representative_point()

def solve():
    # Load layers
    official = gpd.read_file("nyc_parks.gpkg", layer="parks_official")
    osm = gpd.read_file("nyc_parks.gpkg", layer="parks_osm")

    # Ensure CRS matches
    if official.crs != osm.crs:
        osm = osm.to_crs(official.crs)

    # Dissolve layers to get single geometry per side
    # (In case official/osm have multiple polygons that should be treated as one unioned area)
    official_union = official.geometry.unary_union
    osm_union = osm.geometry.unary_union

    # We want areas that exist in ONLY one.
    # diff_official = official - osm
    # diff_osm = osm - official
    diff_official = official_union.difference(osm_union)
    diff_osm = osm_union.difference(official_union)

    # Now we have two sets of geometries.
    # The user wants "adjacent disagreement patches merged together".
    # This implies that if a diff_official polygon is adjacent to a diff_osm polygon,
    # they should be merged.
    # But the prompt says: "One MultiPolygon per merged group with a source attribute..."
    # AND "exactly one of 'parks_official', 'parks_osm', or 'both'".
    # Wait, if it's "claimed by exactly one side", it CANNOT be 'both'.
    # Re-reading: "I need every patch claimed by exactly one side... source attribute... exactly one of parks_official, parks_osm, or both".
    # If it is claimed by "exactly one side", it's either official or osm.
    # Why is 'both' mentioned? 
    # "those exact strings, matching the input GPKG layer names".
    # "exactly one of parks_official, parks_osm, or both (those exact strings...)"
    # If the requirement is "claimed by exactly one side", then 'both' shouldn't appear in the disagreement file.
    # UNLESS 'both' refers to something else? 
    # "I need every patch claimed by exactly one side... One MultiPolygon per merged group with a source attribute whose value is exactly one of parks_official, parks_osm, or both"
    # If the patch is claimed by both, it is NOT a disagreement.
    # Maybe "disagreement" means (Official UNION OSM) - (Official INTERSECTION OSM)? No, that's symmetric difference.
    # Symmetric difference is (Official - OSM) U (OSM - Official).
    # This is exactly "claimed by exactly one side".
    # So 'both' might be a distractor or for a different task, or I'm misinterpreting.
    # Let's look at the requirement again.
    # "I need every patch claimed by exactly one side, with adjacent disagreement patches merged together"
    # If they are adjacent, and one is Official-only and one is OSM-only, they merge into one MultiPolygon.
    # This merged MultiPolygon would have a 'source'... but which one?
    # "exactly one of parks_official, parks_osm, or both"
    # If a merged group contains both Official-only and OSM-only parts, what is its source?
    # If the source can be 'both', then the "merged group" is likely the union of all disagreement parts that touch.
    # Let's check: If Official-only is adjacent to OSM-only, they form a single group.
    # This group's source would be 'both'.
    
    # Let's refine the logic:
    # 1. Compute symmetric difference: SymDiff = (Official \ OSM) U (OSM \ Official)
    # 2. However, we need to identify groups of adjacent components of the symmetric difference.
    # 3. A component is 'parks_official' if it's only in official.
    # 4. A component is 'parks_osm' if it's only in osm.
    # 5. A component is 'both' if it's a union of adjacent official-only and osm-only parts.
    
    # Wait, "adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature".
    # Example: Official is a square. OSM is the same square shifted 1m right.
    # Official - OSM = a thin strip on the left.
    # OSM - Official = a thin strip on the right.
    # These two strips are adjacent (they share a boundary).
    # If we merge them, we get one feature (the union of the two strips).
    # The prompt says the source for this merged group should be 'both'.

    # Step 1: Get the parts
    # To handle adjacency correctly, we should work with the geometries of the symmetric difference.
    # But we need to know which part came from where.
    
    # Let's use the individual polygons from the input if possible, or just use the difference operations.
    # Actually, the best way to find "adjacent" patches is to take all the pieces of the symmetric difference
    # and cluster them if they touch.
    
    # Let's get the polygons of (Official - OSM) and (OSM - Official).
    # We must ensure we are working with the polygons after the difference.
    
    from shapely.ops import unary_union
    
    # Get individual parts of (Official - OSM)
    # We use the union of all official polygons first to ensure we don't have issues with internal overlaps.
    off_geom = official.geometry.unary_union
    osm_geom = osm.geometry.unary_union
    
    only_off_geom = off_geom.difference(osm_geom)
    only_osm_geom = osm_geom.difference(off_geom)
    
    # Now we need to "merge adjacent disagreement patches".
    # Disagreement patches are the components of only_off_geom and only_osm_geom.
    # Let's extract all individual polygons from both.
    
    def get_polygons(geom):
        if geom.is_empty:
            return []
        if isinstance(geom, Polygon):
            return [geom]
        if isinstance(geom, MultiPolygon):
            return list(geom.geoms)
        return []

    parts_off = get_polygons(only_off_geom)
    parts_osm = get_polygons(only_osm_geom)
    
    # We need to group them. Two parts are in the same group if they touch.
    # This is a graph problem.
    
    all_parts = []
    sources = []
    for p in parts_off:
        all_parts.append(p)
        sources.append('parks_official')
    for p in parts_osm:
        all_parts.append(p)
        sources.append('parks_osm')
        
    if not all_parts:
        # Create empty GeoDataFrames if no disagreement
        gpd.GeoDataFrame(columns=['source'], geometry='geometry', crs=official.crs).to_file("parks_disagreement.geojson", driver='GeoJSON')
        gpd.GeoDataFrame(columns=['source', 'geometry'], geometry='geometry', crs=official.crs).to_file("park_label_anchors.geojson", driver='GeoJSON')
        return

    n = len(all_parts)
    adj = [[] for _ in range(n)]
    for i in range(n):
        for j in range(i + 1, n):
            # "adjacent" usually means sharing a boundary (intersection is not empty and not just a point/line)
            # but for polygons, "touches" in shapely means intersection is non-empty and in the interior of at least one.
            # However, in GIS, "adjacent" often includes sharing an edge. 
            # Shapely's `touches` is exactly "interiors do not intersect, but boundaries do".
            # Since we are dealing with the difference of two sets, the boundaries of these patches
            # MUST be the boundary of the original union or the boundary of the intersection.
            # The patches from `only_off` and `only_osm` will touch along the boundary of the intersection.
            if all_parts[i].touches(all_parts[j]):
                adj[i].append(j)
                adj[j].append(i)

    # Find connected components
    visited = [False] * n
    groups = []
    
    for i in range(n):
        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)

    # For each group, merge geometries and determine source
    final_geoms = []
    final_sources = []
    final_anchors = []

    for group in groups:
        group_geoms = [all_parts[idx] for idx in group]
        merged_geom = unary_union(group_geoms)
        
        # Determine source
        group_sources = set(sources[idx] for idx in group)
        if len(group_sources) > 1:
            source_val = 'both'
        elif len(group_sources) == 1:
            source_val = list(group_sources)[0]
        else:
            source_val = 'unknown' # Should not happen
            
        final_geoms.append(merged_geom)
        final_sources.append(source_val)
        
        # Anchor point: strictly inside
        # representative_point() is guaranteed to be in the interior.
        anchor = merged_geom.representative_point()
        final_anchors.append({'geometry': anchor, 'source': source_val})

    # Create GeoDataFrames
    gdf_disagreement = gpd.GeoDataFrame({'source': final_sources}, geometry=final_geoms, crs=official.crs)
    gdf_disagreement.to_file("parks_disagreement.geojson", driver='GeoJSON')

    # For anchors, we need a GeoDataFrame. 
    # Note: anchors are points, so they should be in a GeoDataFrame.
    gdf_anchors = gpd.GeoDataFrame(final_anchors, crs=official.crs)
    gdf_anchors.to_file("park_label_anchors.geojson", driver='GeoJSON')

solve()