geo-l2-nyc-park-symdiff

Rerun

import geopandas as gpd

def main():
    gpkg_path = 'nyc_parks.gpkg'
    
    # Read layers
    official = gpd.read_file(gpkg_path, layer='parks_official')
    osm = gpd.read_file(gpkg_path, layer='parks_osm')
    
    print(f"Official: {len(official)} features")
    print(f"OSM: {len(osm)} features")
    
    # Ensure CRS matches
    if official.crs != osm.crs:
        osm = osm.to_crs(official.crs)

    # We want to find symmetric difference.
    # The problem asks for "every patch claimed by exactly one side".
    # This is the symmetric difference: (Official - OSM) union (OSM - Official)
    # But it's more complex because shapes might be slightly shifted.
    
    # Let's try the standard approach:
    # 1. Union all official parks -> multipolygon_official
    # 2. Union all osm parks -> multipolygon_osm
    # 3. Difference: official_only = multipolygon_official.difference(multipolygon_osm)
    # 4. Difference: osm_only = multipolygon_osm.difference(multipolygon_official)
    
    official_union = official.geometry.unary_union
    osm_union = osm.geometry.unary_union
    
    official_only_geom = official_union.difference(osm_union)
    osm_only_geom = osm_union.difference(official_union)
    
    # Now we have two sets of geometries.
    # However, the prompt says "adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature".
    # This implies that if an official park is slightly larger than an OSM park, the "extra" part of the official park (the gain)
    # and the "missing" part of the OSM park (the loss) are actually the same area? 
    # Wait, no. If a park is shifted, the area where it *was* but isn't anymore (loss) is OSM_only.
    # The area where it *is* but wasn't before (gain) is Official_only.
    # "so that a shifted park's gain and loss form a single feature".
    # This means we need to merge the official_only and osm_only if they are "adjacent" or "connected" in some way.
    # But they are disjoint by definition of symmetric difference.
    # Wait, if they are disjoint, they can't be "merged together" unless they touch.
    # If they touch, it means the boundary between official and osm is a line.
    # If they touch, the union of official_only and osm_only will merge them into one feature if we use buffer(0) or union.
    
    # Let's refine:
    # We want the symmetric difference, but we want to group the "gain" and "loss" parts if they are adjacent.
    # The symmetric difference is (Official \ OSM) U (OSM \ Official).
    # These two sets are disjoint. 
    # If a feature in (Official \ OSM) touches a feature in (OSM \ Official), they should be one feature.
    
    sym_diff_geom = official_only_geom.union(osm_only_geom)
    # Since official_only_geom and osm_only_geom are disjoint, 
    # their union will only merge them if they share a boundary.
    # If they share a boundary, the union of them will result in a single geometry.
    
    # Now we need to identify which part is official and which is osm.
    # 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, "exactly one side" refers to the disagreement patches.
    # If a patch is claimed by "both", it's not a disagreement patch.
    # "I need every patch claimed by exactly one side" -> This is exactly the symmetric difference.
    # "adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature"
    # This implies the "source" for a merged feature (gain + loss) might be... what?
    # If it's a gain and a loss, it's not "exactly one side".
    # Re-reading: "One MultiPolygon per merged group with a `source` attribute whose value is exactly one of `parks_official`, `parks_osm`, or `both`"
    # If the group contains BOTH a gain and a loss, what is its source?
    # "exactly one of `parks_official`, `parks_osm`, or `both`"
    # If it's a shifted park, the area is claimed by Official (gain) and OSM (loss).
    # If the merged feature represents both, then the source is 'both'? 
    # But the prompt says "every patch claimed by exactly one side".
    # If a patch is claimed by both, it's not in the disagreement set.
    # This is slightly contradictory. 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."
    # If a feature is a "shifted park", the "gain" part is claimed by `parks_official` (it's in official but not osm).
    # The "loss" part is claimed by `parks_osm` (it's in osm but not official).
    # If they are merged into one feature, what is its `source`?
    # The prompt says: "a `source` attribute whose value is exactly one of `parks_official`, `parks_osm`, or `both`".
    # If the merged feature contains patches from both, then `source` should be `both`.
    
    # Let's trace:
    # 1. official_only = official_union.difference(osm_union)
    # 2. osm_only = osm_union.difference(official_union)
    # 3. These are the "disagreement patches".
    # 4. We need to merge them if they are adjacent.
    # 5. For each merged group:
    #    - if it only contains parts from official_only -> source = 'parks_official'
    #    - if it only contains parts from osm_only -> source = 'parks_osm'
    #    - if it contains parts from both -> source = 'both'
    
    # Implementation detail:
    # Sym_diff = official_only | osm_only.
    # To merge adjacent ones:
    # We can use a small buffer to merge things that touch.
    # Actually, a simple union of the two sets of polygons (treated as individual features) 
    # where we merge those that touch should work.
    
    # Let's use a more robust way to merge:
    # 1. Get all individual polygons from official_only and osm_only.
    # 2. For each polygon, tag it as 'parks_official' or 'parks_osm'.
    # 3. Merge polygons that touch.
    
    # To get individual polygons from a geometry (which might be a MultiPolygon):
    from shapely.geometry import Polygon, MultiPolygon

    def get_polys(geom):
        if isinstance(geom, Polygon):
            return [geom]
        elif isinstance(geom, MultiPolygon):
            return list(geom.geoms)
        else:
            return []

    official_polys = []
    for p in get_polys(official_only_geom):
        official_polys.append({'geometry': p, 'source': 'parks_official'})
        
    osm_polys = []
    for p in get_polys(osm_only_geom):
        osm_polys.append({'geometry': p, 'source': 'parks_osm'})
        
    all_polys = official_polys + osm_polys
    gdf_all = gpd.GeoDataFrame(all_polys, crs=official.crs)
    
    # To merge adjacent ones, we can use a spatial join or a trick.
    # A common trick is to buffer by a tiny amount and then union.
    # But we want to preserve the source.
    
    # Let's use a clustering approach.
    # 1. Build a graph where nodes are polygons and edges exist if they touch.
    # 2. Find connected components.
    # 3. For each component, union the geometries and determine the source.
    
    import networkx as nx
    
    G = nx.Graph()
    for i in range(len(gdf_all)):
        G.add_node(i)
        
    # Spatial index for efficiency
    sindex = gdf_all.sindex
    for i, row in gdf_all.iterrows():
        possible_matches_index = list(sindex.intersection(row.geometry.bounds))
        for j in possible_matches_index:
            if i < j:
                if row.geometry.touches(gdf_all.iloc[j].geometry) or row.geometry.intersects(gdf_all.iloc[j].geometry):
                    # Using intersects to be safe, but touches is more precise for adjacency.
                    # Given they are from symmetric difference, they shouldn't overlap.
                    # They can only touch at the boundary.
                    G.add_edge(i, j)
                    
    components = list(nx.connected_components(G))
    
    merged_features = []
    for comp in components:
        comp_gdf = gdf_all.iloc[list(comp)]
        merged_geom = comp_gdf.geometry.unary_union
        
        # Determine source
        sources = set(comp_gdf['source'].unique())
        if len(sources) > 1:
            source = 'both'
        else:
            source = list(sources)[0]
            
        # Label anchor: a point strictly inside
        # representative_point() is guaranteed to be inside
        anchor = merged_geom.representative_point()
        
        merged_features.append({
            'geometry': merged_geom,
            'source': source,
            'anchor': anchor
        })

    # Create GeoJSONs
    # parks_disagreement.geojson
    disagreement_gdf = gpd.GeoDataFrame([
        {'geometry': f['geometry'], 'source': f['source']} for f in merged_features
    ], crs=official.crs)
    disagreement_gdf.to_file('parks_disagreement.geojson', driver='GeoJSON')
    
    # park_label_anchors.geojson
    anchors_gdf = gpd.GeoDataFrame([
        {'geometry': f['anchor'], 'source': f['source']} for f in merged_features
    ], crs=official.crs)
    anchors_gdf.to_file('park_label_anchors.geojson', driver='GeoJSON')
    
    print(f"Created parks_disagreement.geojson with {len(disagreement_gdf)} features")
    print(f"Created park_label_anchors.geojson with {len(anchors_gdf)} features")

if __name__ == "__main__":
    main()

import geopandas as gpd
import networkx as nx
from shapely.geometry import Polygon, MultiPolygon

def main():
    gpkg_path = 'nyc_parks.gpkg'
    
    # Read layers
    official = gpd.read_file(gpkg_path, layer='parks_official')
    osm = gpd.read_file(gpkg_path, layer='parks_osm')
    
    print(f"Official: {len(official)} features")
    print(f"OSM: {len(osm)} features")
    
    # Ensure CRS matches
    if official.crs != osm.crs:
        osm = osm.to_crs(official.crs)

    # Use unary_union to get a single geometry for each set
    official_union = official.geometry.unary_union
    osm_union = osm.geometry.unary_union
    
    # The disagreement patches are the symmetric difference
    # official_only = official_union - osm_union
    # osm_only = osm_union - official_union
    official_only_geom = official_union.difference(osm_union)
    osm_only_geom = osm_union.difference(official_union)
    
    def get_polys(geom):
        if isinstance(geom, Polygon):
            return [geom]
        elif isinstance(geom, MultiPolygon):
            return list(geom.geoms)
        else:
            return []

    all_polys_data = []
    
    # Get official patches
    for p in get_polys(official_only_geom):
        all_polys_data.append({'geometry': p, 'source': 'parks_official'})
        
    # Get osm patches
    for p in get_polys(osm_only_geom):
        all_polys_data.append({'geometry': p, 'source': 'parks_osm'})
        
    if not all_polys_data:
        print("No disagreement found.")
        return

    gdf_all = gpd.GeoDataFrame(all_polys_data, crs=official.crs)
    
    # Merge adjacent ones
    # Build a graph of adjacency
    G = nx.Graph()
    n = len(gdf_all)
    for i in range(n):
        G.add_node(i)
        
    sindex = gdf_all.sindex
    for i in range(n):
        geom_i = gdf_all.iloc[i].geometry
        # Get potential neighbors using spatial index
        possible_matches_index = list(sindex.intersection(geom_i.bounds))
        for j in possible_matches_index:
            if i < j:
                geom_j = gdf_all.iloc[j].geometry
                # They are disagreement patches, so they shouldn't overlap, only touch.
                # Using intersects/touches/overlaps to catch adjacency.
                # A small buffer might be safer for floating point precision if they should touch.
                if geom_i.intersects(geom_j):
                    G.add_edge(i, j)
                    
    components = list(nx.connected_components(G))
    
    merged_features = []
    for comp in components:
        comp_indices = list(comp)
        comp_gdf = gdf_all.iloc[comp_indices]
        merged_geom = comp_gdf.geometry.unary_union
        
        # Determine source
        sources = set(comp_gdf['source'].unique())
        if len(sources) > 1:
            source = 'both'
        else:
            source = list(sources)[0]
            
        # Label anchor: a point strictly inside
        anchor = merged_geom.representative_point()
        
        merged_features.append({
            'geometry': merged_geom,
            'source': source,
            'anchor': anchor
        })

    # Create GeoJSONs
    # 1. parks_disagreement.geojson
    disagreement_data = []
    for f in merged_features:
        disagreement_data.append({'geometry': f['geometry'], 'source': f['source']})
    disagreement_gdf = gpd.GeoDataFrame(disagreement_data, crs=official.crs)
    disagreement_gdf.to_file('parks_disagreement.geojson', driver='GeoJSON')
    
    # 2. park_label_anchors.geojson
    anchor_data = []
    for f in merged_features:
        anchor_data.append({'geometry': f['anchor'], 'source': f['source']})
    anchors_gdf = gpd.GeoDataFrame(anchor_data, crs=official.crs)
    anchors_gdf.to_file('park_label_anchors.geojson', driver='GeoJSON')
    
    print(f"Created parks_disagreement.geojson with {len(disagreement_gdf)} features")
    print(f"Created park_label_anchors.geojson with {len(anchors_gdf)} features")

if __name__ == "__main__":
    main()

import geopandas as gpd
import networkx as nx
from shapely.geometry import Polygon, MultiPolygon

def main():
    gpkg_path = '/work/nyc_parks.gpkg'
    
    # Read layers
    official = gpd.read_file(gpkg_path, layer='parks_official')
    osm = gpd.read_file(gpkg_path, layer='parks_osm')
    
    print(f"Official: {len(official)} features")
    print(f"OSM: {len(osm)} features")
    
    # Ensure CRS matches
    if official.crs != osm.crs:
        osm = osm.to_crs(official.crs)

    # Use unary_union to get a single geometry for each set
    official_union = official.geometry.unary_union
    osm_union = osm.geometry.unary_union
    
    # The disagreement patches are the symmetric difference
    official_only_geom = official_union.difference(osm_union)
    osm_only_geom = osm_union.difference(official_union)
    
    def get_polys(geom):
        if isinstance(geom, Polygon):
            return [geom]
        elif isinstance(geom, MultiPolygon):
            return list(geom.geoms)
        else:
            return []

    all_polys_data = []
    
    # Get official patches
    for p in get_polys(official_only_geom):
        all_polys_data.append({'geometry': p, 'source': 'parks_official'})
        
    # Get osm patches
    for p in get_polys(osm_only_geom):
        all_polys_data.append({'geometry': p, 'source': 'parks_osm'})
        
    if not all_polys_data:
        print("No disagreement found.")
        return

    gdf_all = gpd.GeoDataFrame(all_polys_data, crs=official.crs)
    
    # Merge adjacent ones
    G = nx.Graph()
    n = len(gdf_all)
    for i in range(n):
        G.add_node(i)
        
    sindex = gdf_all.sindex
    for i in range(n):
        geom_i = gdf_all.iloc[i].geometry
        possible_matches_index = list(sindex.intersection(geom_i.bounds))
        for j in possible_matches_index:
            if i < j:
                geom_j = gdf_all.iloc[j].geometry
                if geom_i.intersects(geom_j):
                    G.add_edge(i, j)
                    
    components = list(nx.connected_components(G))
    
    merged_features = []
    for comp in components:
        comp_indices = list(comp)
        comp_gdf = gdf_all.iloc[comp_indices]
        merged_geom = comp_gdf.geometry.unary_union
        
        sources = set(comp_gdf['source'].unique())
        if len(sources) > 1:
            source = 'both'
        else:
            source = list(sources)[0]
            
        anchor = merged_geom.representative_point()
        
        merged_features.append({
            'geometry': merged_geom,
            'source': source,
            'anchor': anchor
        })

    # 1. parks_disagreement.geojson
    disagreement_data = []
    for f in merged_features:
        disagreement_data.append({'geometry': f['geometry'], 'source': f['source']})
    disagreement_gdf = gpd.GeoDataFrame(disagreement_data, crs=official.crs)
    disagreement_gdf.to_file('/work/parks_disagreement.geojson', driver='GeoJSON')
    
    # 2. park_label_anchors.geojson
    anchor_data = []
    for f in merged_features:
        anchor_data.append({'geometry': f['anchor'], 'source': f['source']})
    anchors_gdf = gpd.GeoDataFrame(anchor_data, crs=official.crs)
    anchors_gdf.to_file('/work/park_label_anchors.geojson', driver='GeoJSON')
    
    print(f"Created /work/parks_disagreement.geojson with {len(disagreement_gdf)} features")
    print(f"Created /work/park_label_anchors.geojson with {len(anchors_gdf)} features")

if __name__ == "__main__":
    main()

import geopandas as gpd
import networkx as nx
from shapely.geometry import Polygon, MultiPolygon
import os

def main():
    gpkg_path = '/work/nyc_parks.gpkg'
    
    if not os.path.exists(gpkg_path):
        print(f"Error: {gpkg_path} not found")
        return

    # Read layers
    official = gpd.read_file(gpkg_path, layer='parks_official')
    osm = gpd.read_file(gpkg_path, layer='parks_osm')
    
    print(f"Official: {len(official)} features")
    print(f"OSM: {len(osm)} features")
    
    # Ensure CRS matches
    if official.crs != osm.crs:
        osm = osm.to_crs(official.crs)

    # Use unary_union to get a single geometry for each set
    official_union = official.geometry.unary_union
    osm_union = osm.geometry.unary_union
    
    # The disagreement patches are the symmetric difference
    official_only_geom = official_union.difference(osm_union)
    osm_only_geom = osm_union.difference(official_union)
    
    def get_polys(geom):
        if isinstance(geom, Polygon):
            return [geom]
        elif isinstance(geom, MultiPolygon):
            return list(geom.geoms)
        else:
            return []

    all_polys_data = []
    
    # Get official patches
    for p in get_polys(official_only_geom):
        if not p.is_empty:
            all_polys_data.append({'geometry': p, 'source': 'parks_official'})
        
    # Get osm patches
    for p in get_polys(osm_only_geom):
        if not p.is_empty:
            all_polys_data.append({'geometry': p, 'source': 'parks_osm'})
        
    if not all_polys_data:
        print("No disagreement found.")
        return

    gdf_all = gpd.GeoDataFrame(all_polys_data, crs=official.crs)
    
    # Merge adjacent ones
    G = nx.Graph()
    n = len(gdf_all)
    for i in range(n):
        G.add_node(i)
        
    sindex = gdf_all.sindex
    for i in range(n):
        geom_i = gdf_all.iloc[i].geometry
        possible_matches_index = list(sindex.intersection(geom_i.bounds))
        for j in possible_matches_index:
            if i < j:
                geom_j = gdf_all.iloc[j].geometry
                if geom_i.intersects(geom_j):
                    G.add_edge(i, j)
                    
    components = list(nx.connected_components(G))
    
    merged_features = []
    for comp in components:
        comp_indices = list(comp)
        comp_gdf = gdf_all.iloc[comp_indices]
        merged_geom = comp_gdf.geometry.unary_union
        
        sources = set(comp_gdf['source'].unique())
        if len(sources) > 1:
            source = 'both'
        else:
            source = list(sources)[0]
            
        anchor = merged_geom.representative_point()
        
        merged_features.append({
            'geometry': merged_geom,
            'source': source,
            'anchor': anchor
        })

    # 1. parks_disagreement.geojson
    disagreement_data = []
    for f in merged_features:
        disagreement_data.append({'geometry': f['geometry'], 'source': f['source']})
    disagreement_gdf = gpd.GeoDataFrame(disagreement_data, crs=official.crs)
    disagreement_gdf.to_file('/work/parks_disagreement.geojson', driver='GeoJSON')
    
    # 2. park_label_anchors.geojson
    anchor_data = []
    for f in merged_features:
        anchor_data.append({'geometry': f['anchor'], 'source': f['source']})
    anchors_gdf = gpd.GeoDataFrame(anchor_data, crs=official.crs)
    anchors_gdf.to_file('/work/park_label_anchors.geojson', driver='GeoJSON')
    
    print(f"Created /work/parks_disagreement.geojson with {len(disagreement_gdf)} features")
    print(f"Created /work/park_label_anchors.geojson with {len(anchors_gdf)} features")

if __name__ == "__main__":
    main()

import geopandas as gpd
import networkx as nx
from shapely.geometry import Polygon, MultiPolygon
import os

def main():
    gpkg_path = '/work/nyc_parks.gpkg'
    
    if not os.path.exists(gpkg_path):
        print(f"Error: {gpkg_path} not found")
        return

    # Read layers
    official = gpd.read_file(gpkg_path, layer='parks_official')
    osm = gpd.read_file(gpkg_path, layer='parks_osm')
    
    print(f"Official: {len(official)} features")
    print(f"OSM: {len(osm)} features")
    
    # Ensure CRS matches
    if official.crs != osm.crs:
        osm = osm.to_crs(official.crs)

    # Use unary_union to get a single geometry for each set
    official_union = official.geometry.unary_union
    osm_union = osm.geometry.unary_union
    
    # The disagreement patches are the symmetric difference
    official_only_geom = official_union.difference(osm_union)
    osm_only_geom = osm_union.difference(official_union)
    
    def get_polys(geom):
        if isinstance(geom, Polygon):
            return [geom]
        elif isinstance(geom, MultiPolygon):
            return list(geom.geoms)
        else:
            return []

    all_polys_data = []
    
    # Get official patches
    for p in get_polys(official_only_geom):
        if not p.is_empty:
            all_polys_data.append({'geometry': p, 'source': 'parks_official'})
        
    # Get osm patches
    for p in get_polys(osm_only_geom):
        if not p.is_empty:
            all_polys_data.append({'geometry': p, 'source': 'parks_osm'})
        
    if not all_polys_data:
        print("No disagreement found.")
        return

    gdf_all = gpd.GeoDataFrame(all_polys_data, crs=official.crs)
    
    # Merge adjacent ones
    G = nx.Graph()
    n = len(gdf_all)
    for i in range(n):
        G.add_node(i)
        
    sindex = gdf_all.sindex
    for i in range(n):
        geom_i = gdf_all.iloc[i].geometry
        possible_matches_index = list(sindex.intersection(geom_i.bounds))
        for j in possible_matches_index:
            if i < j:
                geom_j = gdf_all.iloc[j].geometry
                if geom_i.intersects(geom_j):
                    G.add_edge(i, j)
                    
    components = list(nx.connected_components(G))
    
    merged_features = []
    for comp in components:
        comp_indices = list(comp)
        comp_gdf = gdf_all.iloc[comp_indices]
        merged_geom = comp_gdf.geometry.unary_union
        
        sources = set(comp_gdf['source'].unique())
        if len(sources) > 1:
            source = 'both'
        else:
            source = list(sources)[0]
            
        anchor = merged_geom.representative_point()
        
        merged_features.append({
            'geometry': merged_geom,
            'source': source,
            'anchor': anchor
        })

    # 1. parks_disagreement.geojson
    disagreement_data = []
    for f in merged_features:
        disagreement_data.append({'geometry': f['geometry'], 'source': f['source']})
    disagreement_gdf = gpd.GeoDataFrame(disagreement_data, crs=official.crs)
    disagreement_gdf.to_file('/work/parks_disagreement.geojson', driver='GeoJSON')
    
    # 2. park_label_anchors.geojson
    anchor_data = []
    for f in merged_features:
        anchor_data.append({'geometry': f['anchor'], 'source': f['source']})
    anchors_gdf = gpd.GeoDataFrame(anchor_data, crs=official.crs)
    anchors_gdf.to_file('/work/park_label_anchors.geojson', driver='GeoJSON')
    
    print(f"Created /work/parks_disagreement.geojson with {len(disagreement_gdf)} features")
    print(f"Created /work/park_label_anchors.geojson with {len(anchors_gdf)} features")

if __name__ == "__main__":
    main()

import geopandas as gpd
import networkx as nx
from shapely.geometry import Polygon, MultiPolygon
import os

def main():
    gpkg_path = '/work/nyc_parks.gpkg'
    
    if not os.path.exists(gpkg_path):
        print(f"Error: {gpkg_path} not found")
        return

    # Read layers
    official = gpd.read_file(gpkg_path, layer='parks_official')
    osm = gpd.read_file(gpkg_path, layer='parks_osm')
    
    print(f"Official: {len(official)} features")
    print(f"OSM: {len(osm)} features")
    
    # Ensure CRS matches
    if official.crs != osm.crs:
        osm = osm.to_crs(official.crs)

    # Use unary_union to get a single geometry for each set
    official_union = official.geometry.unary_union
    osm_union = osm.geometry.unary_union
    
    # The disagreement patches are the symmetric difference
    official_only_geom = official_union.difference(osm_union)
    osm_only_geom = osm_union.difference(official_union)
    
    def get_polys(geom):
        if isinstance(geom, Polygon):
            return [geom]
        elif isinstance(geom, MultiPolygon):
            return list(geom.geoms)
        else:
            return []

    all_polys_data = []
    
    # Get official patches
    for p in get_polys(official_only_geom):
        if not p.is_empty:
            all_polys_data.append({'geometry': p, 'source': 'parks_official'})
        
    # Get osm patches
    for p in get_polys(osm_only_geom):
        if not p.is_empty:
            all_polys_data.append({'geometry': p, 'source': 'parks_osm'})
        
    if not all_polys_data:
        print("No disagreement found.")
        return

    gdf_all = gpd.GeoDataFrame(all_polys_data, crs=official.crs)
    
    # Merge adjacent ones
    G = nx.Graph()
    n = len(gdf_all)
    for i in range(n):
        G.add_node(i)
        
    sindex = gdf_all.sindex
    for i in range(n):
        geom_i = gdf_all.iloc[i].geometry
        possible_matches_index = list(sindex.intersection(geom_i.bounds))
        for j in possible_matches_index:
            if i < j:
                geom_j = gdf_all.iloc[j].geometry
                if geom_i.intersects(geom_j):
                    G.add_edge(i, j)
                    
    components = list(nx.connected_components(G))
    
    merged_features = []
    for comp in components:
        comp_indices = list(comp)
        comp_gdf = gdf_all.iloc[comp_indices]
        merged_geom = comp_gdf.geometry.unary_union
        
        sources = set(comp_gdf['source'].unique())
        if len(sources) > 1:
            source = 'both'
        else:
            source = list(sources)[0]
            
        anchor = merged_geom.representative_point()
        
        merged_features.append({
            'geometry': merged_geom,
            'source': source,
            'anchor': anchor
        })

    # 1. parks_disagreement.geojson
    disagreement_data = []
    for f in merged_features:
        disagreement_data.append({'geometry': f['geometry'], 'source': f['source']})
    disagreement_gdf = gpd.GeoDataFrame(disagreement_data, crs=official.crs)
    disagreement_gdf.to_file('/work/parks_disagreement.geojson', driver='GeoJSON')
    
    # 2. park_label_anchors.geojson
    anchor_data = []
    for f in merged_features:
        anchor_data.append({'geometry': f['anchor'], 'source': f['source']})
    anchors_gdf = gpd.GeoDataFrame(anchor_data, crs=official.crs)
    anchors_gdf.to_file('/work/park_label_anchors.geojson', driver='GeoJSON')
    
    print(f"Created /work/parks_disagreement.geojson with {len(disagreement_gdf)} features")
    print(f"Created /work/park_label_anchors.geojson with {len(anchors_gdf)} features")

if __name__ == "__main__":
    main()