geo-l2-nyc-park-symdiff

Rerun

import geopandas as gpd

def check_layers(path):
    import fiona
    layers = fiona.listlayers(path)
    print(f"Layers in {path}: {layers}")

if __name__ == "__main__":
    check_layers('/work/nyc_parks.gpkg')

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

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

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

    # 2. Get geometry only to simplify processing
    official_geom = official.geometry.unary_union
    osm_geom = osm.geometry.unary_union

    # 3. Find discrepancies
    # Parts in official but not in osm
    diff_official = official_geom.difference(osm_geom)
    # Parts in osm but not in official
    diff_osm = osm_geom.difference(official_geom)

    # A "patch" is a piece of geometry that exists in one but not both.
    # We need to find these and then merge adjacent ones.
    # However, a shifted park's gain and loss might be adjacent.
    # Let's represent the disagreement as a set of polygons.
    
    # To handle "adjacent disagreement patches merged together", we can do the following:
    # We have the set of all areas that are in (Official XOR OSM).
    # Let's find the union of these differences.
    xor_geom = official_geom.symmetric_difference(osm_geom)
    
    # If xor_geom is a single geometry, we might want to break it into its constituent 
    # connected components (polygons/multipolygons).
    # But wait, a "shifted park" might consist of a piece that is 'official' and 
    # a piece that is 'osm', and they might be touching.
    # If they touch, we want them merged into one feature with a specific source?
    # "One MultiPolygon per merged group with a source attribute whose value is exactly one of parks_official, parks_osm, or both"
    # Wait, the instruction says: "exactly one of parks_official, parks_osm, or both".
    # If a group is a merger of an 'official' patch and an 'osm' patch, 
    # what is the source? 
    # "One MultiPolygon per merged group with a source attribute whose value is exactly one of parks_official, parks_osm, or both".
    # If it's a merger of an official part and an osm part, it should be 'both'?
    # No, "exactly one of ...". This implies if it contains both, it's 'both'.
    # Let's re-read: "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 shifted park has a 'gain' (in osm) and a 'loss' (in official), they are adjacent.
    # Merging them results in a single geometry.
    # The source of this merged geometry should be 'both' because it contains parts from both.
    
    # Let's refine the process:
    # 1. Find all 'only official' polygons: official_only = official_geom - osm_geom
    # 2. Find all 'only osm' polygons: osm_only = osm_geom - official_geom
    # 3. Combine them into a single collection of polygons.
    # 4. Dissolve them based on adjacency.
    
    # Actually, the XOR geometry itself is the collection of all "disagreement patches".
    # If we take the connected components of the XOR geometry, we get the "merged groups".
    
    from shapely.geometry import MultiPolygon
    
    def get_connected_components(geometry):
        if geometry.is_empty:
            return []
        if isinstance(geometry, Polygon):
            return [geometry]
        if isinstance(geometry, MultiPolygon):
            return list(geometry.geoms)
        # In case of GeometryCollection
        if hasattr(geometry, 'geoms'):
            components = []
            for geom in geometry.geoms:
                components.extend(get_connected_components(geom))
            return components
        return [geometry]

    components = get_connected_components(xor_geom)
    
    disagreement_features = []
    anchor_features = []

    for comp in components:
        if comp.is_empty:
            continue
            
        # Determine source
        # A component is 'parks_official' if it's entirely within official_geom AND not in osm_geom?
        # No, the component is a subset of the XOR.
        # If the component is entirely contained in official_geom (but not in osm_geom), it's 'parks_official'.
        # If it's entirely contained in osm_geom (but not in official_geom), it's 'parks_osm'.
        # If it contains parts of both, it's 'both'.
        
        # Since the component is part of the symmetric difference, any part of it is either
        # in official OR in osm, but NOT both.
        # So a component can be a mix of 'only official' and 'only osm' parts.
        
        # Let's check parts of the component.
        # We can use intersection to see if the component has parts in official or osm.
        # But we must be careful: the component is already in the XOR, so it cannot be in 'both' (intersection).
        # Wait. If it's in the XOR, it is (Official \ OSM) OR (OSM \ Official).
        # Thus, it cannot have an intersection with (Official INTERSECT OSM).
        # So any part of the component is either in (Official \ OSM) or (OSM \ Official).
        
        # If a component is composed of both 'official only' parts and 'osm only' parts,
        # it means they are touching.
        
        in_official = comp.intersection(official_geom)
        in_osm = comp.intersection(osm_geom)
        
        # Note: Because comp is a subset of (official XOR osm), 
        # in_official and in_osm are disjoint.
        
        if not in_official.is_empty and not in_osm.is_empty:
            source = 'both'
        elif not in_official.is_empty:
            source = 'parks_official'
        elif not in_osm.is_empty:
            source = 'parks_osm'
        else:
            continue # Should not happen

        # Prepare disagreement feature
        # The requirement says "One MultiPolygon per merged group".
        # We ensure comp is treated as MultiPolygon if it's not.
        if isinstance(comp, Polygon):
            poly_to_save = MultiPolygon([comp])
        else:
            poly_to_save = comp
            
        disagreement_features.append({
            'geometry': poly_to_save,
            'source': source
        })
        
        # Prepare anchor point
        # "a label-anchor point that falls strictly inside the group's geometry"
        # representative_point() is a good candidate as it's guaranteed to be inside.
        anchor_point = comp.representative_point()
        anchor_features.append({
            'geometry': anchor_point,
            'source': source
        })

    # 4. Save to GeoJSON
    gdf_disagreement = gpd.GeoDataFrame(disagreement_features, crs=official.crs)
    gdf_disagreement.to_file('/work/parks_disagreement.geojson', driver='GeoJSON')

    gdf_anchors = gpd.GeoDataFrame(anchor_features, crs=official.crs)
    gdf_anchors.to_file('/work/park_label_anchors.geojson', driver='GeoJSON')

    print(f"Created parks_disagreement.geojson with {len(gdf_disagreement)} features.")
    print(f"Created park_label_anchors.geojson with {len(gdf_anchors)} features.")

if __name__ == "__main__":
    solve()

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

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

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

    # 2. Get geometry only to simplify processing
    official_geom = official.geometry.union_all()
    osm_geom = osm.geometry.union_all()

    # 3. Find discrepancies
    # We want the symmetric difference (XOR)
    xor_geom = official_geom.symmetric_difference(osm_geom)
    
    from shapely.geometry import MultiPolygon
    
    def get_connected_components(geometry):
        if geometry.is_empty:
            return []
        if isinstance(geometry, Polygon):
            return [geometry]
        if isinstance(geometry, MultiPolygon):
            return list(geometry.geoms)
        if hasattr(geometry, 'geoms'):
            components = []
            for geom in geometry.geoms:
                components.extend(get_connected_components(geom))
            return components
        return [geometry]

    components = get_connected_components(xor_geom)
    
    disagreement_features = []
    anchor_features = []

    for comp in components:
        if comp.is_empty:
            continue
            
        # Check if this component contains parts of both official and osm
        # We need to check against the original geometries. 
        # However, we must check against the 'original' polygons because 
        # official_geom and osm_geom are already unions.
        # Actually, the symmetric difference handles the logic.
        # A component is either in (Official \ OSM) or (OSM \ Official) OR a combination.
        # To be a combination, it must touch a boundary.
        
        # Let's use a small buffer or checking intersection with the original sets.
        # But since we already have the union, let's intersect the component with official_geom and osm_geom.
        # However, the XOR means component is (Official AND NOT OSM) OR (OSM AND NOT Official).
        # Wait, if a component contains both, it must be that some part is in (Official \ OSM) 
        # and some part is in (OSM \ Official).
        # These two parts are disjoint. For them to be in the same component, 
        # they must touch.
        
        in_official = comp.intersection(official_geom)
        in_osm = comp.intersection(osm_geom)
        
        # We need to check if the intersection is non-trivial.
        # Because they are part of the symmetric difference, in_official and in_osm 
        # cannot overlap.
        
        has_official = not in_official.is_empty and in_official.area > 1e-7
        has_osm = not in_osm.is_empty and in_osm.area > 1e-7
        
        if has_official and has_osm:
            source = 'both'
        elif has_official:
            source = 'parks_official'
        elif has_osm:
            source = 'parks_osm'
        else:
            # This case might happen due to precision issues
            continue

        # Prepare disagreement feature
        if isinstance(comp, Polygon):
            poly_to_save = MultiPolygon([comp])
        else:
            poly_to_save = comp
            
        disagreement_features.append({
            'geometry': poly_to_save,
            'source': source
        })
        
        # Prepare anchor point
        anchor_point = comp.representative_point()
        anchor_features.append({
            'geometry': anchor_point,
            'source': source
        })

    # 4. Save to GeoJSON
    gdf_disagreement = gpd.GeoDataFrame(disagreement_features, crs=official.crs)
    gdf_disagreement.to_file('/work/parks_disagreement.geojson', driver='GeoJSON')

    gdf_anchors = gpd.GeoDataFrame(anchor_features, crs=official.crs)
    gdf_anchors.to_file('/work/park_label_anchors.geojson', driver='GeoJSON')

    print(f"Created parks_disagreement.geojson with {len(gdf_disagreement)} features.")
    print(f"Created park_label_anchors.geojson with {len(gdf_anchors)} features.")

if __name__ == "__main__":
    solve()

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

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

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

    # 2. Get geometry only to simplify processing
    official_geom = official.geometry.union_all()
    osm_geom = osm.geometry.union_all()

    # 3. Find discrepancies
    # XOR is the set of all areas that are in exactly one of the two.
    xor_geom = official_geom.symmetric_difference(osm_geom)
    
    from shapely.geometry import MultiPolygon
    
    def get_connected_components(geometry):
        if geometry.is_empty:
            return []
        if isinstance(geometry, Polygon):
            return [geometry]
        if isinstance(geometry, MultiPolygon):
            return list(geometry.geoms)
        if hasattr(geometry, 'geoms'):
            components = []
            for geom in geometry.geoms:
                components.extend(get_connected_components(geom))
            return components
        return [geometry]

    # Get components of the XOR.
    # However, if an 'official' patch and an 'osm' patch are adjacent,
    # they might be part of the same connected component in the XOR 
    # IF they touch. But in XOR, they only exist where one is present and the other isn't.
    # If an official polygon and an osm polygon share a boundary, 
    # the XOR geometry will have those two polygons touching at that boundary.
    # Thus they will be part of the same connected component.
    
    components = get_connected_components(xor_geom)
    
    disagreement_features = []
    anchor_features = []

    for comp in components:
        if comp.is_empty:
            continue
            
        # Determine source
        # A component is 'both' if it contains parts of both official and osm.
        # Because it's from the XOR, it cannot contain parts that are in BOTH official and osm.
        # It contains parts that are in (Official \ OSM) and parts that are in (OSM \ Official).
        
        in_official = comp.intersection(official_geom)
        in_osm = comp.intersection(osm_geom)
        
        # Check if both intersections are non-empty (using area to avoid precision issues)
        has_official = in_official.area > 1e-6
        has_osm = in_osm.area > 1e-6
        
        if has_official and has_osm:
            source = 'both'
        elif has_official:
            source = 'parks_official'
        elif has_osm:
            source = 'parks_osm'
        else:
            continue

        # Prepare disagreement feature
        # Ensure it's a MultiPolygon as requested
        if isinstance(comp, Polygon):
            poly_to_save = MultiPolygon([comp])
        else:
            poly_to_save = comp
            
        disagreement_features.append({
            'geometry': poly_to_save,
            'source': source
        })
        
        # Prepare anchor point
        anchor_point = comp.representative_point()
        anchor_features.append({
            'geometry': anchor_point,
            'source': source
        })

    # 4. Save to GeoJSON
    gdf_disagreement = gpd.GeoDataFrame(disagreement_features, crs=official.crs)
    # Sort to ensure matching order for anchors (though representative_point is computed per feature)
    gdf_disagreement = gdf_disagreement.sort_index()
    gdf_disagreement.to_file('/work/parks_disagreement.geojson', driver='GeoJSON')

    gdf_anchors = gpd.GeoDataFrame(anchor_features, crs=official.crs)
    gdf_anchors = gdf_anchors.sort_index()
    gdf_anchors.to_file('/work/park_label_anchors.geojson', driver='GeoJSON')

    print(f"Created parks_disagreement.geojson with {len(gdf_disagreement)} features.")
    print(f"Created park_label_anchors.geojson with {len(gdf_anchors)} features.")

if __name__ == "__main__":
    solve()