geo-l2-nyc-park-symdiff

Rerun

# /// script
# requires-python = ">=3.12"
# dependencies = ["geopandas", "shapely", "pyproj"]
# ///

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

def get_label_anchor(geom):
    """Returns a point strictly inside the geometry."""
    if geom.is_empty:
        return None
    
    # Try representative_point first as it's guaranteed to be inside
    point = geom.representative_point()
    
    # If for some reason it's still not strictly inside (e.g. extremely thin/complex)
    # we'll just return it. For most polygons, representative_point is perfect.
    return point

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

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

    # 2. Identify discrepancies
    # We need to find patches that are in official but NOT in osm,
    # and patches that are in osm but NOT in official.
    
    # Union of all official geometries
    official_union = official.geometry.unary_union
    # Union of all osm geometries
    osm_union = osm.geometry.unary_union

    # Patches only in official: official - osm
    only_official_geom = official_union.difference(osm_union)
    
    # Patches only in osm: osm - official
    only_osm_geom = osm_union.difference(official_union)

    # 3. Combine and dissolve adjacent discrepancies
    # The prompt says: "adjacent disagreement patches merged together so that a shifted 
    # park's gain and loss form a single feature rather than two."
    # Wait, if a park is shifted, the gain and loss are adjacent. 
    # If we merge them, we get the "disagreement area".
    
    # Actually, "adjacent disagreement patches merged together" implies we should look 
    # at the set of all disagreement areas (the symmetric difference) 
    # and merge any that touch.
    
    disagreement_geom = official_union.symmetric_difference(osm_union)
    
    # Since symmetric_difference returns a geometry that might be a MultiPolygon,
    # we can use unary_union on the individual parts to ensure adjacency is respected.
    # However, symmetric_difference already returns the union of the disjoint parts.
    # To "merge adjacent disagreement patches", we need to ensure that if a patch 
    # of 'only_official' touches a patch of 'only_osm', they are merged.
    
    # Let's refine:
    # A 'disagreement' is any area where they don't overlap perfectly.
    # The prompt wants "every patch claimed by exactly one side, with adjacent 
    # disagreement patches merged together".
    
    # This means:
    # 1. Calculate the symmetric difference (the area of disagreement).
    # 2. The problem is that the symmetric difference doesn't tell us which side "claimed" it.
    #    But the prompt says: "One MultiPolygon per merged group with a source attribute 
    #    whose value is exactly one of parks_official, parks_osm, or both".
    #    WAIT. If they are merged, how can the source be "exactly one of"?
    #    Ah: "One MultiPolygon per merged group with a source attribute whose value is exactly 
    #    one of parks_official, parks_osm, or both".
    #    If a merged group contains both official-only and osm-only parts, the source must be 'both'? 
    #    No, the prompt says: "every patch claimed by exactly one side".
    #    This suggests the patches themselves are the objects of interest.
    #    "with adjacent disagreement patches merged together so that a shifted park's gain and 
    #    loss form a single feature rather than two."
    
    # Let's re-read carefully: "I need every patch claimed by exactly one side, with adjacent 
    # disagreement patches merged together so that a shifted park's gain and loss form a 
    # single feature rather than two."
    # This implies if an official park is shifted, the 'gap' (official but not osm) and 
    # the 'overlap' (osm but not official) are adjacent and should be one feature.
    # BUT then the 'source' cannot be "exactly one of parks_official, parks_osm, or both".
    # Let's look at the source values: `parks_official`, `parks_osm`, or `both`.
    
    # If I have a group of polygons that touch, and some are official-only and some are osm-only,
    # and I merge them, the resulting MultiPolygon is "both"? 
    # Or does "both" mean something else?
    
    # Let's re-read: "One MultiPolygon per merged group with a source attribute whose value 
    # is exactly one of parks_official, parks_osm, or both (those exact strings, 
    # matching the input GPKG layer names)."
    
    # If a merged group contains parts from both, the source is 'both'. 
    # If it contains parts only from official, it's 'parks_official'.
    # If it contains parts only from osm, it's 'parks_osm'.

    # Wait, if they are merged, they must be adjacent. 
    # If they are adjacent, and they are "disagreement patches", 
    # then a shifted park's gain (official only) and loss (osm only) 
    # will definitely touch. 
    # Thus they will be merged into one feature.
    # And that feature's source will be 'both'.

    # 1. Get the patches.
    # patches_official = official_union - osm_union
    # patches_osm = osm_union - official_union
    
    # 2. We need to identify which patches are adjacent.
    # The easiest way is to take the union of (patches_official) and (patches_osm)
    # but we need to keep track of the source.
    
    # Actually, if we merge all adjacent disagreement patches, 
    # we can just take the `symmetric_difference` and then 
    # dissolve them into connected components? 
    # No, `symmetric_difference` already gives us the disjoint parts.
    # We need to "merge adjacent disagreement patches".
    # This is basically finding the connected components of the disagreement area.
    
    # Let's try this:
    # 1. Compute symmetric difference: `sym_diff = official_union ^ osm_union`
    # 2. This `sym_diff` is a collection of polygons. 
    # 3. Some are official-only, some are osm-only.
    # 4. We need to group these polygons into connected components.
    #    Two polygons are in the same component if they touch.
    # 5. For each component, determine if it's 'parks_official', 'parks_osm', or 'both'.

    # To find connected components of polygons:
    # We can use a graph where each polygon is a node and edges exist if they touch.
    # Or, more simply, we can use `unary_union` on all polygons, but that just gives 
    # the same union. 
    # Wait, if they touch, `unary_union` will merge them. 
    # If I have a list of polygons that touch, `unary_union` will return a single 
    # MultiPolygon (or Polygon) representing their union.
    # That's exactly what "merged together" means.

    # So the algorithm:
    # 1. Find all "official-only" polygons.
    # 2. Find all "osm-only" polygons.
    # 3. We need to group them if they touch.
    # 4. How to group them? 
    #    We can use `shapely.ops.unary_union` on the union of ALL disagreement polygons?
    #    No, that just gives the whole symmetric difference.
    #    We want to group them into *connected components*.

    # To find connected components:
    # Use a library or a trick. A trick: 
    # Create a collection of all the individual polygons (from both official-only and osm-only).
    # Then use a graph approach.

    from shapely.ops import unary_union
    from collections import deque

    # Step 1: Get the individual polygons from the symmetric difference
    # (using the fact that symmetric difference might return a MultiPolygon)
    
    # Let's get the constituent polygons of the symmetric difference.
    # First, we need to know which part is which.
    
    # official_only = official_union.difference(osm_union)
    # osm_only = osm_union.difference(official_union)
    
    # A better way:
    # 1. Get all the "atoms" of the disagreement.
    # We can use the union of (official_union) and (osm_union) boundaries? No.
    
    # Let's use the difference approach.
    # official_only_geom = official_union.difference(osm_union)
    # osm_only_geom = osm_union.difference(official_union)
    
    # To get individual polygons from a MultiPolygon:
    def get_polygons(geom):
        if geom.is_empty:
            return []
        if geom.geom_type == 'Polygon':
            return [geom]
        if geom.geom_type == 'MultiPolygon':
            return list(geom.geoms)
        return []

    # Note: symmetric difference might contain holes, but unary_union handles that.
    # However, official_only_geom might be a MultiPolygon.
    
    # Let's build a list of (geometry, source)
    parts = []
    for p in get_polygons(official_union.difference(osm_union)):
        parts.append({'geometry': p, 'source': 'parks_official'})
    for p in get_polygons(osm_union.difference(official_union)):
        parts.append({'geometry': p, 'source': 'parks_osm'})

    # Now we need to group parts that touch.
    # We'll use a Disjoint Set Union (DSU) or a simple adjacency graph.
    
    n = len(parts)
    adj = [[] for _ in range(n)]
    for i in range(n):
        for j in range(i + 1, n):
            if parts[i]['geometry'].touches(parts[j]['geometry']) or \
               parts[i]['geometry'].intersects(parts[j]['geometry']):
                # Actually, 'intersects' is safer, but 'touches' is for boundary.
                # Since they are from a symmetric difference, they can only touch at boundaries.
                # But if they share an edge, 'intersects' is true.
                adj[i].append(j)
                adj[j].append(i)

    # Find connected components
    visited = [False] * n
    components = []
    for i in range(n):
        if not visited[i]:
            comp = []
            q = deque([i])
            visited[i] = True
            while q:
                curr = q.popleft()
                comp.append(curr)
                for neighbor in adj[curr]:
                    if not visited[neighbor]:
                        visited[neighbor] = True
                        q.append(neighbor)
            components.append(comp)

    # For each component, merge geometries and determine source
    final_features = []
    for comp_indices in components:
        comp_geoms = [parts[idx]['geometry'] for idx in comp_indices]
        merged_geom = unary_union(comp_geoms)
        
        # Determine source
        sources = set(parts[idx]['source'] for idx in comp_indices)
        if len(sources) > 1:
            source_val = 'both'
        else:
            # If it's 'parks_official' or 'parks_osm'
            # Note: the prompt says "exactly one of parks_official, parks_osm, or both"
            # and "matching the input GPKG layer names".
            # The layer names were 'parks_official' and 'parks_osm'.
            source_val = list(sources)[0]
            
        final_features.append({
            'geometry': merged_geom,
            'source': source_val
        })

    # 4. Create GeoDataFrames
    disagreement_gdf = gpd.GeoDataFrame(final_features, crs=official.crs)
    
    # Sort to ensure deterministic order for anchors (optional but good)
    disagreement_gdf = disagreement_gdf.sort_values('source').reset_index(drop=True)

    # 5. Create anchors
    anchors = []
    for idx, row in disagreement_gdf.iterrows():
        anchor_point = get_label_anchor(row['geometry'])
        anchors.append({
            'geometry': anchor_point,
            'source': row['source']
        })
    
    anchors_gdf = gpd.GeoDataFrame(anchors, crs=disagreement_gdf.crs)

    # 6. Write outputs
    disagreement_gdf.to_file("parks_disagreement.geojson", driver="GeoJSON")
    anchors_gdf.to_file("park_label_anchors.geojson", driver="GeoJSON")
    
    print(f"Created {len(disagreement_gdf)} disagreement features.")
    print(f"Created {len(anchors_gdf)} anchors.")

if __name__ == "__main__":
    solve()