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():
    # Load the layers
    official = gpd.read_file("nyc_parks.gpkg", layer="parks_official")
    osm = gpd.read_file("nyc_parks.gpkg", layer="parks_osm")

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

    # We need to find areas covered by ONLY one of the two.
    # A "patch" is a piece of geometry.
    # To do this systematically:
    # 1. Get the union of all boundaries (to find all atomic polygons)
    # 2. Or more simply:
    #    Area in official but not in osm: official.difference(osm)
    #    Area in osm but not in official: osm.difference(official)
    
    # However, 'difference' on a MultiPolygon might return many tiny pieces.
    # The requirement says: "adjacent disagreement patches merged together 
    # so that a shifted park's gain and loss form a single feature"
    # This implies we should look at the symmetric difference: (A - B) U (B - A)
    # BUT we want to keep track of which side it belongs to.
    
    # Let's define the disagreement regions.
    # A disagreement region is a connected component in the set difference.
    # Wait, if a park is shifted, the "gain" (OSM part not in Official) 
    # and the "loss" (Official part not in OSM) are adjacent.
    # They share a boundary.
    
    # Let's construct a geometry that represents all disagreement areas.
    # We want to find connected components of the symmetric difference.
    # But we need to label them with 'source'.
    # If a component contains parts of 'official - osm' AND 'osm - official', 
    # we need to decide if it's 'both'? 
    # The prompt says: "One MultiPolygon per merged group with a `source` 
    # attribute whose value is exactly one of `parks_official`, `parks_osm`, or `both`".
    # Actually, if they are adjacent, they form one feature. 
    # If that feature contains parts from both, the source is 'both'.
    
    # Algorithm:
    # 1. Compute Diff_Official = Official - OSM
    # 2. Compute Diff_OSM = OSM - Official
    # 3. Combine them: All_Diff = Diff_Official.union(Diff_OSM)
    # 4. Find connected components of All_Diff.
    # 5. For each component:
    #    a. Check if it intersects (is contained in) Diff_Official.
    #    b. Check if it intersects (is contained in) Diff_OSM.
    #    c. Assign source:
    #       - If only in Diff_Official -> 'parks_official'
    #       - If only in Diff_OSM -> 'parks_osm'
    #       - If in both -> 'both'
    #    d. Find an anchor point.

    # Step 1 & 2
    diff_official = official.unary_union.difference(osm.unary_union)
    diff_osm = osm.unary_union.difference(official.unary_union)
    
    # Step 3
    all_diff_union = diff_official.union(diff_osm)
    
    # Step 4: Extract individual polygons/components
    # unary_union of a geometry can still be a single MultiPolygon.
    # We want to explode it into individual pieces to find connected components.
    # Actually, unary_union of the differences might already be separated.
    # Let's use a trick: collect all polygons from the difference and group them.
    
    def get_polygons(geom):
        if geom.is_empty:
            return []
        if isinstance(geom, Polygon):
            return [geom]
        if isinstance(geom, MultiPolygon):
            return list(geom.geoms)
        return []

    # Note: official.unary_union might be a MultiPolygon. 
    # We need to handle the case where diff_official or diff_osm are MultiPolygons.
    
    # To properly get connected components, we can use the fact that 
    # shared boundaries mean they are part of the same "disagreement group".
    # The simplest way to find connected components of the symmetric difference 
    # is to dissolve the union of the differences. 
    # But we want to preserve the distinction.
    
    # Let's create a single geometry representing all disagreement.
    # We can use a "buffered" approach or just use the polygons from unary_union.
    # If we take the union of (Diff_Official + Diff_OSM), 
    # any parts that touch will be merged into one polygon.
    
    combined_diff = diff_official.union(diff_osm)
    
    # The 'combined_diff' might be a Polygon or MultiPolygon.
    # If it's a MultiPolygon, each part is a connected component.
    
    components = []
    if combined_diff.is_empty:
        print("No disagreements found.")
        return
    
    if hasattr(combined_diff, 'geoms'):
        components = list(combined_diff.geoms)
    else:
        components = [combined_diff]

    results = []
    for comp in components:
        # Determine source
        in_official = comp.intersects(diff_official)
        in_osm = comp.intersects(diff_osm)
        
        # We use a small buffer or check intersection with a point to be safe,
        # but intersection of geometries is fine. 
        # Since 'comp' is a connected component of the union, 
        # it will either be entirely in one, or it will span the boundary.
        
        # Note: 'comp.intersects(diff_official)' might be true if it just touches.
        # But since 'comp' is a component of the symmetric difference, 
        # if it's 'both', it must have interior points in both.
        # Let's check if the intersection has area.
        
        has_official_area = diff_official.intersection(comp).area > 1e-9
        has_osm_area = diff_osm.intersection(comp).area > 1e-9
        
        if has_official_area and has_osm_area:
            source = 'both'
        elif has_official_area:
            source = 'parks_official'
        elif has_osm_area:
            source = 'parks_osm'
        else:
            continue # Should not happen
            
        # Find anchor point: a point strictly inside
        # representative_point() is a good candidate as it is guaranteed to be inside.
        anchor = comp.representative_point()
        
        results.append({
            'geometry': comp,
            'source': source,
            'anchor': anchor
        })

    # Create GeoDataFrame for disagreements
    disagreement_gdf = gpd.GeoDataFrame(results, crs=official.crs)
    # The requirement asks for 'source' column in the GeoJSON.
    # The 'anchor' is not in the first GeoJSON.
    
    # Output 1: parks_disagreement.geojson
    disagreement_gdf[['geometry', 'source']].to_file("parks_disagreement.geojson", driver="GeoJSON")
    
    # Output 2: park_label_anchors.geojson
    # "exactly one anchor per disagreement group in matching order"
    anchors_gdf = gpd.GeoDataFrame({
        'source': disagreement_gdf['source'],
        'geometry': disagreement_gdf['anchor']
    }, crs=official.crs)
    anchors_gdf.to_file("park_label_anchors.geojson", driver="GeoJSON")
    
    print(f"Processed {len(results)} disagreement groups.")

if __name__ == "__main__":
    solve()

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

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

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

    # 1. Compute the actual difference areas
    # We want to find where they don't overlap.
    # official_only = Official - OSM
    # osm_only = OSM - Official
    
    # Using union_all() if available (newer geopandas) or unary_union
    try:
        official_geom = official.union_all()
    except AttributeError:
        official_geom = official.unary_union
        
    try:
        osm_geom = osm.union_all()
    except AttributeError:
        osm_geom = osm.unary_union

    diff_official = official_geom.difference(osm_geom)
    diff_osm = osm_geom.difference(official_geom)
    
    # 2. The disagreement area is the union of these differences
    # We want to group adjacent parts. 
    # If a piece of 'official_only' is adjacent to a piece of 'osm_only', 
    # they belong to the same disagreement group.
    
    # To find these groups, we take the union of the two differences.
    # This will merge any parts that share a boundary.
    combined_diff = diff_official.union(diff_osm)
    
    if combined_diff.is_empty:
        print("No disagreements found.")
        return

    # 3. Extract individual components (connected components)
    if hasattr(combined_diff, 'geoms'):
        components = list(combined_diff.geoms)
    else:
        components = [combined_diff]

    results = []
    for comp in components:
        # Determine source:
        # Use intersection with a small area to verify if it contains parts of both
        # or just one.
        
        # We check if the component has area in diff_official and diff_osm.
        # We use a tiny buffer or intersection to ensure we are looking at the interior.
        has_official_area = diff_official.intersection(comp).area > 1e-7
        has_osm_area = diff_osm.intersection(comp).area > 1e-7
        
        if has_official_area and has_osm_area:
            source = 'both'
        elif has_official_area:
            source = 'parks_official'
        elif has_osm_area:
            source = 'parks_osm'
        else:
            # This might happen due to floating point precision near boundaries
            # but in a symmetric difference union, it shouldn't.
            continue
            
        # Find anchor point: a point strictly inside
        anchor = comp.representative_point()
        
        results.append({
            'geometry': comp,
            'source': source,
            'anchor': anchor
        })

    if not results:
        print("No disagreement components extracted.")
        return

    # Create GeoDataFrame for disagreements
    disagreement_gdf = gpd.GeoDataFrame(results, crs=official.crs)
    
    # Output 1: parks_disagreement.geojson
    # Requirement: "One MultiPolygon per merged group with a `source` attribute"
    # We need to ensure geometry is MultiPolygon if it contains multiple parts,
    # but 'comp' is already a single connected component. 
    # However, to be safe and follow "One MultiPolygon per merged group", 
    # if comp is a Polygon, we can keep it or cast it. 
    # Standard practice is to allow both, but if the user insists on MultiPolygon:
    def ensure_multipolygon(geom):
        if isinstance(geom, Polygon):
            return MultiPolygon([geom])
        return geom

    disagreement_gdf['geometry'] = disagreement_gdf['geometry'].apply(ensure_multipolygon)
    disagreement_gdf[['geometry', 'source']].to_file("parks_disagreement.geojson", driver="GeoJSON")
    
    # Output 2: park_label_anchors.geojson
    # "exactly one anchor per disagreement group in matching order"
    anchors_gdf = gpd.GeoDataFrame({
        'source': disagreement_gdf['source'],
        'geometry': disagreement_gdf['anchor']
    }, crs=official.crs)
    anchors_gdf.to_file("park_label_anchors.geojson", driver="GeoJSON")
    
    print(f"Processed {len(results)} disagreement groups.")
    print(f"Sources found: {disagreement_gdf['source'].unique()}")

if __name__ == "__main__":
    solve()

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

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

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

    # 1. Compute the actual difference areas
    try:
        official_geom = official.union_all()
    except AttributeError:
        official_geom = official.unary_union
        
    try:
        osm_geom = osm.union_all()
    except AttributeError:
        osm_geom = osm.unary_union

    # Diff_Official: in official but not in osm
    # Diff_OSM: in osm but not in official
    diff_official = official_geom.difference(osm_geom)
    diff_osm = osm_geom.difference(official_geom)
    
    # 2. Grouping adjacent disagreement patches.
    # We want to find connected components of (Diff_Official UNION Diff_OSM).
    # If a piece of Diff_Official is adjacent to Diff_OSM, they form one component.
    # To find these, we can use the union of the two.
    # However, 'difference' results can sometimes have shared boundaries that are 
    # technically not "touching" due to precision. 
    # But since these are derived from the same input geometries, 
    # the boundaries should align.
    
    combined_diff_geom = diff_official.union(diff_osm)
    
    if combined_diff_geom.is_empty:
        print("No disagreements found.")
        return

    # 3. Extract individual components
    if hasattr(combined_diff_geom, 'geoms'):
        components = list(combined_diff_geom.geoms)
    else:
        components = [combined_diff_geom]

    results = []
    for comp in components:
        # Use a small buffer or intersection to determine if it's 'both'
        # The component is a connected component of the symmetric difference.
        # If it contains parts from BOTH diff_official and diff_osm, it's 'both'.
        
        # Important: check intersection with the actual difference geometries.
        # Intersection with the original diff_official might pick up boundary lines.
        # We check if the intersection has area.
        inter_official = diff_official.intersection(comp)
        inter_osm = diff_osm.intersection(comp)
        
        has_official_area = inter_official.area > 1e-7
        has_osm_area = inter_osm.area > 1e-7
        
        if has_official_area and has_osm_area:
            source = 'both'
        elif has_official_area:
            source = 'parks_official'
        elif has_osm_area:
            source = 'parks_osm'
        else:
            continue
            
        # Find anchor point: a point strictly inside the component.
        anchor = comp.representative_point()
        
        results.append({
            'geometry': comp,
            'source': source,
            'anchor': anchor
        })

    if not results:
        print("No disagreement components extracted.")
        return

    # Create GeoDataFrame for disagreements
    disagreement_gdf = gpd.GeoDataFrame(results, crs=official.crs)
    
    # Ensure geometry is MultiPolygon as requested (or just Polygon/MultiPolygon)
    # The requirement: "One MultiPolygon per merged group"
    # Note: if a component is a Polygon, we can wrap it in a MultiPolygon.
    def to_multipolygon(geom):
        if isinstance(geom, Polygon):
            return MultiPolygon([geom])
        return geom

    disagreement_gdf['geometry'] = disagreement_gdf['geometry'].apply(to_multipolygon)
    
    # Output 1: parks_disagreement.geojson
    # Need to make sure 'source' column is exactly as required.
    disagreement_gdf[['geometry', 'source']].to_file("parks_disagreement.geojson", driver="GeoJSON")
    
    # Output 2: park_label_anchors.geojson
    # "exactly one anchor per disagreement group in matching order"
    anchors_gdf = gpd.GeoDataFrame({
        'source': disagreement_gdf['source'],
        'geometry': disagreement_gdf['anchor']
    }, crs=official.crs)
    anchors_gdf.to_file("park_label_anchors.geojson", driver="GeoJSON")
    
    print(f"Processed {len(results)} disagreement groups.")
    print(f"Sources found: {disagreement_gdf['source'].unique()}")

if __name__ == "__main__":
    solve()

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

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

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

    # 1. Compute the actual difference areas
    try:
        official_geom = official.union_all()
    except AttributeError:
        official_geom = official.unary_union
        
    try:
        osm_geom = osm.union_all()
    except AttributeError:
        osm_geom = osm.unary_union

    # Diff_Official: in official but not in osm
    # Diff_OSM: in osm but not in official
    diff_official = official_geom.difference(osm_geom)
    diff_osm = osm_geom.difference(official_geom)
    
    # 2. Grouping adjacent disagreement patches.
    # We want to find connected components of (Diff_Official UNION Diff_OSM).
    # To find these, we can use the union of the two.
    combined_diff_geom = diff_official.union(diff_osm)
    
    if combined_diff_geom.is_empty:
        print("No disagreements found.")
        return

    # 3. Extract individual components
    if hasattr(combined_diff_geom, 'geoms'):
        components = list(combined_diff_geom.geoms)
    else:
        components = [combined_diff_geom]

    results = []
    for comp in components:
        # We need to check if this component 'comp' belongs to 'official', 'osm', or 'both'.
        # A component is 'both' if it contains area from both diff_official and diff_osm.
        
        # The intersection might be a bit finicky with precision.
        # Let's use a point-in-polygon approach on the representative point 
        # or just check the area of intersection.
        
        # To avoid precision issues, let's use a small buffer when checking intersection
        # but we actually want to know if the *interior* of 'comp' overlaps with the 
        # interior of the diffs.
        
        # Intersection area check:
        inter_official_area = diff_official.intersection(comp).area
        inter_osm_area = diff_osm.intersection(comp).area
        
        # If the component is 'both', it must have non-zero area in both.
        # Using a very small epsilon.
        eps = 1e-6
        if inter_official_area > eps and inter_osm_area > eps:
            source = 'both'
        elif inter_official_area > eps:
            source = 'parks_official'
        elif inter_osm_area > eps:
            source = 'parks_osm'
        else:
            # Try a representative point check if area check fails due to precision
            rp = comp.representative_point()
            # We need to be careful, representative_point is inside the geometry.
            # If it's on the boundary of the diffs, it might not work.
            # But since 'comp' is a component of the union, if it's a 'both' component,
            # its interior must contain points from both.
            # We'll check if the rp is in diff_official or diff_osm.
            # We use a tiny buffer to include points on the boundary.
            if diff_official.intersects(rp):
                source = 'parks_official'
            elif diff_osm.intersects(rp):
                source = 'parks_osm'
            else:
                continue # Should not happen
            
        # Find anchor point: a point strictly inside the component.
        anchor = comp.representative_point()
        
        results.append({
            'geometry': comp,
            'source': source,
            'anchor': anchor
        })

    if not results:
        print("No disagreement components extracted.")
        return

    # Create GeoDataFrame for disagreements
    disagreement_gdf = gpd.GeoDataFrame(results, crs=official.crs)
    
    # Ensure geometry is MultiPolygon as requested
    def to_multipolygon(geom):
        if isinstance(geom, Polygon):
            return MultiPolygon([geom])
        return geom

    disagreement_gdf['geometry'] = disagreement_gdf['geometry'].apply(to_multipolygon)
    
    # Output 1: parks_disagreement.geojson
    disagreement_gdf[['geometry', 'source']].to_file("parks_disagreement.geojson", driver="GeoJSON")
    
    # Output 2: park_label_anchors.geojson
    anchors_gdf = gpd.GeoDataFrame({
        'source': disagreement_gdf['source'],
        'geometry': disagreement_gdf['anchor']
    }, crs=official.crs)
    anchors_gdf.to_file("park_label_anchors.geojson", driver="GeoJSON")
    
    print(f"Processed {len(results)} disagreement groups.")
    print(f"Sources found: {disagreement_gdf['source'].unique()}")

if __name__ == "__main__":
    solve()

from shapely.geometry import Polygon; p1 = Polygon([(0,0), (1,0), (1,1), (0,1)]); p2 = Polygon([(1,0), (2,0), (2,1), (1,1)]); union = p1.union(p2); print('Union intersects p1:', union.intersection(p1).area > 1e-7); print('Union intersects p2:', union.intersection(p2).area > 1e-7)"

from shapely.geometry import Polygon
p1 = Polygon([(0,0), (1,0), (1,1), (0,1)])
p2 = Polygon([(1,0), (2,0), (2,1), (1,1)])
union = p1.union(p2)
print('Union area:', union.area)
print('Union intersects p1:', union.intersection(p1).area > 1e-7)
print('Union intersects p2:', union.intersection(p2).area > 1e-7)

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

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

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

    # 1. Compute the actual difference areas
    try:
        official_geom = official.union_all()
    except AttributeError:
        official_geom = official.unary_union
        
    try:
        osm_geom = osm.union_all()
    except AttributeError:
        osm_geom = osm.unary_union

    # Diff_Official: in official but not in osm
    # Diff_OSM: in osm but not in official
    diff_official = official_geom.difference(osm_geom)
    diff_osm = osm_geom.difference(official_geom)
    
    # 2. Grouping adjacent disagreement patches.
    # We want to find connected components of (Diff_Official UNION Diff_OSM).
    combined_diff_geom = diff_official.union(diff_osm)
    
    if combined_diff_geom.is_empty:
        print("No disagreements found.")
        return

    # 3. Extract individual components
    if hasattr(combined_diff_geom, 'geoms'):
        components = list(combined_diff_geom.geoms)
    else:
        components = [combined_diff_geom]

    results = []
    for comp in components:
        # A component is 'both' if it contains area from both diff_official and diff_osm.
        inter_official_area = diff_official.intersection(comp).area
        inter_osm_area = diff_osm.intersection(comp).area
        
        # Check for 'both' first
        if inter_official_area > 1e-6 and inter_osm_area > 1e-6:
            source = 'both'
        elif inter_official_area > 1e-6:
            source = 'parks_official'
        elif inter_osm_area > 1e-6:
            source = 'parks_osm'
        else:
            # Fallback: check representative point
            rp = comp.representative_point()
            if diff_official.contains(rp):
                source = 'parks_official'
            elif diff_osm.contains(rp):
                source = 'parks_osm'
            else:
                # If we still can't tell, maybe it's a line/point but we are dealing with polygons
                continue
            
        anchor = comp.representative_point()
        results.append({
            'geometry': comp,
            'source': source,
            'anchor': anchor
        })

    if not results:
        print("No disagreement components extracted.")
        return

    disagreement_gdf = gpd.GeoDataFrame(results, crs=official.crs)
    
    def to_multipolygon(geom):
        if isinstance(geom, Polygon):
            return MultiPolygon([geom])
        return geom

    disagreement_gdf['geometry'] = disagreement_gdf['geometry'].apply(to_multipolygon)
    disagreement_gdf[['geometry', 'source']].to_file("parks_disagreement.geojson", driver="GeoJSON")
    
    anchors_gdf = gpd.GeoDataFrame({
        'source': disagreement_gdf['source'],
        'geometry': disagreement_gdf['anchor']
    }, crs=official.crs)
    anchors_gdf.to_file("park_label_anchors.geojson", driver="GeoJSON")
    
    print(f"Processed {len(results)} disagreement groups.")
    print(f"Sources found: {disagreement_gdf['source'].unique()}")

if __name__ == "__main__":
    solve()

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

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

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

    # 1. Compute the actual difference areas
    try:
        official_geom = official.union_all()
    except AttributeError:
        official_geom = official.unary_union
        
    try:
        osm_geom = osm.union_all()
    except AttributeError:
        osm_geom = osm.unary_union

    diff_official = official_geom.difference(osm_geom)
    diff_osm = osm_geom.difference(official_geom)
    
    # 2. Grouping adjacent disagreement patches.
    # The union of differences might not merge polygons that only share a boundary 
    # due to floating point issues. 
    # A robust way: buffer the differences slightly, union them, then union that back 
    # with the original differences. 
    # Wait, a simpler way: use a small buffer to join them.
    
    # We use a tiny buffer to ensure they merge.
    # Then we use that merged geometry to find the components.
    # Finally, we use the original diff_official and diff_osm to find the source.
    
    # Buffer by a very small amount to merge adjacent things
    # 1e-4 is usually safe for metric coords (NYC uses UTM or similar)
    buffer_val = 0.01 
    combined_diff_geom = diff_official.buffer(buffer_val).union(diff_osm.buffer(buffer_val))
    
    # Now, the components of combined_diff_geom are our "merged groups".
    # But we must ensure we don't include the buffer in our final geometry.
    # So for each component in combined_diff_geom, we take its intersection 
    # with the actual (unbuffered) union of differences.
    
    actual_disagreement_union = diff_official.union(diff_osm)

    if combined_diff_geom.is_empty:
        print("No disagreements found.")
        return

    if hasattr(combined_diff_geom, 'geoms'):
        buffered_components = list(combined_diff_geom.geoms)
    else:
        buffered_components = [combined_diff_geom]

    results = []
    for b_comp in buffered_components:
        # The actual geometry for this group is the intersection 
        # of the buffered component and the actual differences.
        comp = b_comp.intersection(actual_disagreement_union)
        
        if comp.is_empty or comp.area < 1e-7:
            continue

        # Determine source
        inter_official_area = diff_official.intersection(comp).area
        inter_osm_area = diff_osm.intersection(comp).area
        
        if inter_official_area > 1e-6 and inter_osm_area > 1e-6:
            source = 'both'
        elif inter_official_area > 1e-6:
            source = 'parks_official'
        elif inter_osm_area > 1e-6:
            source = 'parks_osm'
        else:
            # Fallback
            rp = comp.representative_point()
            if diff_official.contains(rp):
                source = 'parks_official'
            elif diff_osm.contains(rp):
                source = 'parks_osm'
            else:
                continue
            
        anchor = comp.representative_point()
        results.append({
            'geometry': comp,
            'source': source,
            'anchor': anchor
        })

    if not results:
        print("No disagreement components extracted.")
        return

    disagreement_gdf = gpd.GeoDataFrame(results, crs=official.crs)
    
    def to_multipolygon(geom):
        if isinstance(geom, Polygon):
            return MultiPolygon([geom])
        return geom

    disagreement_gdf['geometry'] = disagreement_gdf['geometry'].apply(to_multipolygon)
    disagreement_gdf[['geometry', 'source']].to_file("parks_disagreement.geojson", driver="GeoJSON")
    
    anchors_gdf = gpd.GeoDataFrame({
        'source': disagreement_gdf['source'],
        'geometry': disagreement_gdf['anchor']
    }, crs=official.crs)
    anchors_gdf.to_file("park_label_anchors.geojson", driver="GeoJSON")
    
    print(f"Processed {len(results)} disagreement groups.")
    print(f"Sources found: {disagreement_gdf['source'].unique()}")

if __name__ == "__main__":
    solve()