geo-l2-nyc-park-symdiff

Rerun

Let me start by understanding the problem:

1. We have a GeoPackage `nyc_parks.gpkg` with two layers: `parks_official` and `parks_osm`
2. I need to find areas where they disagree (claimed by exactly one side)
3. Adjacent disagreement patches should be merged together
4. Output a MultiPolygon per merged group with a `source` attribute (either `parks_official`, `parks_osm`, or `both`)
5. Also output label anchor points (one per group, strictly inside the geometry)

Let me first explore the data to understand what we're working with.

Let me understand the task:

1. We have two layers in `nyc_parks.gpkg`: `parks_official` and `parks_osm`
2. We need to find areas claimed by exactly one side (not both) - i.e., the symmetric difference
3. Adjacent disagreement patches should be merged together
4. Each merged group should have:
   - A `source` attribute: `parks_official` (only in official), `parks_osm` (only in OSM), or `both` (wait, the instruction says "every patch claimed by exactly one side" but then says source can be `both` too... Let me re-read)

"every patch claimed by exactly one side, with adjacent disagreement patches merged together" - So we're looking at the symmetric difference, where a patch is in one layer but not the other. But adjacent disagreement patches (some in official, some in OSM) that are adjacent should be merged... Wait, that doesn't make sense. Let me re-read.

"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."

So if a park in official is at location A, and in OSM it's at location B (shifted), then:
- Area in official but not OSM = loss (official source)
- Area in OSM but not official = gain (OSM source)
- If these are adjacent (because of the shift), they should be merged into one feature with source... hmm, but what source?

Wait, "so that a shifted park's gain and loss form a single feature rather than two" - this means adjacent patches that disagree (one in official only, one in OSM only) that touch should be merged into a single feature. But what `source` would it have?

Let me re-read: "One MultiPolygon per merged group with a source attribute whose value is exactly one of parks_official, parks_osm, or both"

Ah, so the source can be `both` too. That makes sense - if the gain and loss from a shifted park are adjacent and we merge them, the merged feature contains areas from both sources, so source = `both`.

Let me think about the algorithm:

1. Compute the symmetric difference (areas in official only, areas in OSM only)
2. For each polygon in official-only and each polygon in OSM-only, check adjacency (touching boundaries)
3. Merge adjacent disagreement patches into groups
4. For each merged group:
   - If all patches are official-only → source = "parks_official"
   - If all patches are OSM-only → source = "parks_osm"
   - If mixed → source = "both"
5. Create a label anchor point strictly inside each group's geometry

Let me write a Python script to do this.

Actually, let me think more carefully. The symmetric difference of two polygon layers gives:
- official \ osm (areas in official not in OSM)
- osm \ official (areas in OSM not in official)

A "shifted park" means the park exists in both datasets but at slightly different locations. The area where they differ forms a ring-like shape. The gain (OSM-only) and loss (official-only) patches might touch each other (they'd be adjacent along the boundary).

So the algorithm:
1. Compute `official_only = official.difference(osm_union)` (or use overlay)
2. Compute `osm_only = osm.difference(official_union)`
3. Label each polygon with its source
4. Merge adjacent polygons that touch each other, combining their sources
5. For each merged group, determine the source:
   - If all constituent polygons have the same source, use that
   - If mixed, use "both"

Let me implement this.

Actually, I need to be careful with the union/difference operations. Let me use geopandas overlay operations.

Let me write the script:

```python
import geopandas as gpd
from shapely.ops import unary_union, polygonize
from shapely.geometry import MultiPolygon, Polygon, Point
import numpy as np

# Read layers
official = gpd.read_file('nyc_parks.gpkg', layer='parks_official')
osm = gpd.read_file('nyc_parks.gpkg', layer='parks_osm')

# Ensure same CRS - they should be, both EPSG:6539

# Compute the symmetric difference
# official_only = areas in official not in OSM
# osm_only = areas in OSM not in official

# Union all OSM geometries for difference
osm_union = unary_union(osm.geometry)
official_union = unary_union(official.geometry)

# Dissolve each layer first to avoid fragmenting
official_dissolved = official.dissolve()
osm_dissolved = osm.dissolve()

# Now compute
official_only = official_dissolved.geometry.iloc[0].difference(osm_union)
osm_only = osm_dissolved.geometry.iloc[0].difference(official_union)
```

Wait, but dissolving might lose the granularity. Let me think again.

Actually, the simpler approach: compute the overlay and get the symmetric difference.

```python
# Use overlay to find where they differ
from geopandas import overlay

# Areas in official not in osm
official_only = overlay(official, osm, how='difference')
official_only['source'] = 'parks_official'

# Areas in osm not in official
osm_only = overlay(osm, official, how='difference')
osm_only['source'] = 'parks_osm'
```

Then merge adjacent polygons...

Actually, the overlay difference approach should work well. Let me think about the merging step.

For merging adjacent patches:
- We have a GeoDataFrame with polygons from both official_only and osm_only
- We want to merge any polygons that touch (share a boundary)
- For each merged group, determine the source

I can do this with unary union + polygonization, but that's complex. Another approach:

1. Create a GeoSeries of all disagreement polygons with their source
2. Use spatial join/intersection to find which polygons touch each other
3. Build a graph of touching polygons
4. Find connected components
5. Merge each component

Let me implement this step by step.

Actually, a simpler approach: 
1. For each polygon, buffer by 0 to clean geometry, then check adjacency
2. Use the unary_union approach: dissolve all disagreement polygons, then the resulting MultiPolygon's parts naturally merge adjacent ones. But I need to track source per part.

Hmm, let me think differently. 

The approach I'll use:

1. Compute `official_only` using overlay difference
2. Compute `osm_only` using overlay difference
3. Combine them into one GeoDataFrame with a `source` column
4. Use spatial adjacency (touches) to find connected components
5. Merge each component, determine the combined source
6. Compute label anchors

Let me code this up.

For the label anchors: I need a point that falls strictly inside the group's geometry. I can use `representative_point()` from shapely, which guarantees the point is within the geometry (unlike centroid which might be outside).

Let me write the script.

Actually, let me think about edge cases. The `difference` operation between two polygon layers can produce many small sliver polygons due to slight boundary misalignments. But the data is in the same CRS and is park data, so hopefully the boundaries are relatively clean.

Let me code this properly.

```python
import geopandas as gpd
from shapely.ops import unary_union
from shapely.geometry import MultiPolygon, Polygon, Point
import numpy as np
import json

# Read layers
official = gpd.read_file('/work/nyc_parks.gpkg', layer='parks_official')
osm = gpd.read_file('/work/nyc_parks.gpkg', layer='parks_osm')

# Ensure same CRS
crs = official.crs

# Compute symmetric difference
official_only = gpd.overlay(official, osm, how='difference')
official_only['source'] = 'parks_official'

osm_only = gpd.overlay(osm, official, how='difference')
osm_only['source'] = 'parks_osm'

# Combine
disagreement = gpd.GeoDataFrame(pd.concat([official_only[['source', 'geometry']], 
                                            osm_only[['source', 'geometry']]], ignore_index=True),
                                 crs=crs)

# Now merge adjacent polygons
# Build adjacency graph
# For each pair of polygons, check if they touch
# Use R-tree spatial index for efficiency

from shapely import STRtree

# Explode MultiPolygons to individual polygons
# ... or we can work with them as is

# Let's use a simpler approach: buffer by epsilon and union
# Actually, let me use the graph approach
```

Let me think about a cleaner approach. I'll work polygon by polygon.

Actually, the simplest robust approach might be:

1. Create two dissolved layers: official_dissolved and osm_dissolved
2. Compute the symmetric difference as dissolved_official.symmetric_difference(dissolved_osm)
3. This gives a MultiPolygon. Each polygon in it is either:
   - In official only (inside official, outside osm)
   - In OSM only (inside osm, outside official)
4. But this loses the distinction of source... The symmetric difference just gives geometry, not attribution.

OK so let me go back to the overlay approach but with a better merge strategy.

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

# Read
official = gpd.read_file('/work/nyc_parks.gpkg', layer='parks_official')
osm = gpd.read_file('/work/nyc_parks.gpkg', layer='parks_osm')
crs = official.crs

# Dissolve to get the full coverage
official_dissolved = gpd.GeoDataFrame({'geometry': [unary_union(official.geometry)]}, crs=crs)
osm_dissolved = gpd.GeoDataFrame({'geometry': [unary_union(osm.geometry)]}, crs=crs)

# Compute differences
official_only = gpd.overlay(official_dissolved, osm_dissolved, how='difference')
official_only['source'] = 'parks_official'

osm_only = gpd.overlay(osm_dissolved, official_dissolved, how='difference')
osm_only['source'] = 'parks_osm'

# Combine
all_disagreement = pd.concat([official_only[['source', 'geometry']], 
                               osm_only[['source', 'geometry']]], ignore_index=True)
disagreement_gdf = gpd.GeoDataFrame(all_disagreement, crs=crs)

# Now merge adjacent/touching polygons
# Strategy: use unary_union to merge touching polygons
# But we need to track source...

# Alternative: dissolve by source first, then look at adjacency between sources
```

Hmm, this is getting complex. Let me think of a different approach.

The key insight: "so that a shifted park's gain and loss form a single feature rather than two" - this means that if an official-only polygon touches an OSM-only polygon, they should be merged into one feature with source = "both".

So the algorithm is:
1. Find all disagreement polygons (official-only and OSM-only)
2. Build adjacency graph: edges between polygons that touch
3. Find connected components
4. For each component:
   - Merge all geometries
   - Determine source: "parks_official" if all are official-only, "parks_osm" if all are OSM-only, "both" if mixed
   - Compute label anchor

Let me code this properly.

```python
import geopandas as gpd
import pandas as pd
from shapely.ops import unary_union
from shapely.geometry import MultiPolygon, Polygon, Point, shape, mapping
from shapely import STRtree
import json

# Read
official = gpd.read_file('/work/nyc_parks.gpkg', layer='parks_official')
osm = gpd.read_file('/work/nyc_parks.gpkg', layer='parks_osm')
crs = official.crs

# Compute differences using dissolved layers for efficiency
official_union = unary_union(official.geometry)
osm_union = unary_union(osm.geometry)

# We can use shapely directly
official_only_geom = official_union.difference(osm_union)
osm_only_geom = osm_union.difference(official_union)

# Convert to GeoDataFrame with individual polygons
def explode_multi_to_rows(geom, source):
    rows = []
    if geom is None or geom.is_empty:
        return rows
    if geom.geom_type == 'MultiPolygon':
        for poly in geom.geoms:
            rows.append({'source': source, 'geometry': poly})
    elif geom.geom_type == 'Polygon':
        rows.append({'source': source, 'geometry': geom})
    return rows

rows = []
rows.extend(explode_multi_to_rows(official_only_geom, 'parks_official'))
rows.extend(explode_multi_to_rows(osm_only_geom, 'parks_osm'))

disagreement_gdf = gpd.GeoDataFrame(rows, crs=crs)
print(f"Number of disagreement polygons: {len(disagreement_gdf)}")

# Now build adjacency graph
# Use spatial index
sindex = disagreement_gdf.sindex
n = len(disagreement_gdf)

# adjacency as list of sets
adj = [set() for _ in range(n)]

for i, row in disagreement_gdf.iterrows():
    # Find candidates that might touch this polygon
    candidates = list(sindex.intersection(row.geometry.bounds))
    for j in candidates:
        if j <= i:
            continue
        # Check if they touch or intersect
        if row.geometry.touches(disagreement_gdf.iloc[j].geometry) or row.geometry.intersects(disagreement_gdf.iloc[j].geometry):
            # Actually, for "adjacent disagreement patches", they should touch or be very close
            # Let's use touches for strict adjacency
            if row.geometry.touches(disagreement_gdf.iloc[j].geometry):
                adj[i].add(j)
                adj[j].add(i)

# Find connected components
visited = [False] * n
components = []
for i in range(n):
    if not visited[i]:
        component = []
        stack = [i]
        while stack:
            node = stack.pop()
            if not visited[node]:
                visited[node] = True
                component.append(node)
                for neighbor in adj[node]:
                    if not visited[neighbor]:
                        stack.append(neighbor)
        components.append(component)

print(f"Number of connected components: {len(components)}")

# Merge each component
results = []
for comp in components:
    geoms = [disagreement_gdf.iloc[i].geometry for i in comp]
    sources = [disagreement_gdf.iloc[i]['source'] for i in comp]
    
    merged_geom = unary_union(geoms)
    
    # Determine source
    unique_sources = set(sources)
    if len(unique_sources) == 1:
        source = list(unique_sources)[0]
    else:
        source = 'both'
    
    # Label anchor: representative point inside the geometry
    anchor = merged_geom.representative_point()
    
    results.append({'source': source, 'geometry': merged_geom, 'anchor': anchor})

# Create output GeoDataFrames
...
```

Wait, I need to be careful about what "touches" means. Two polygons "touch" if they share a boundary but their interiors don't intersect. For a shifted park, the gain (OSM-only) and loss (official-only) patches would typically touch along a boundary line.

But there could be cases where they're so shifted that they don't touch - in that case they'd remain separate features.

Also, I should consider using a small buffer to handle floating-point precision issues with `touches`. Let me use `touches` with a small tolerance.

Actually, let me also think about whether I need to use `buffer(0)` to clean geometries first.

Let me write the full script and test it.

Actually, wait - I also need to consider that there might be multiple disconnected disagreement patches that all belong to the same "group" (like one park that shifted). But those wouldn't be adjacent, so they'd be separate groups. That's fine as per the instructions.

Let me also think: when doing the difference between dissolved layers, we get the symmetric difference. But `official_union.difference(osm_union)` gives official \ osm, which could be multiple disconnected polygons. These could be:
- Parts of parks that exist in official but not at all in OSM
- Parts of parks that exist in both but the boundaries are slightly different (shifted)

The adjacency merging would catch the shifted case where the official-only and OSM-only patches touch.

Let me also consider: what if there's a park that's completely missing from one dataset? Then the entire park polygon appears in the difference, and since there's no corresponding OSM-only patch touching it, it stays as a single-source feature. Good.

Let me write the complete script now.

One more thought: I should use `buffer(0)` on the geometries to clean any topological issues, and I should also consider using a small epsilon for the touches check.

Let me also consider using `intersects` instead of `touches` - because if there's a slight numerical issue, `touches` might fail. Let me use a small buffer: `geom_a.buffer(1e-8).intersects(geom_b.buffer(1e-8))` or something similar.

Actually, let me use a robust approach: consider two polygons adjacent if their intersection area is negligible (below a threshold like 1 square meter) but they are within a small distance of each other. Or simply use `touches` after `buffer(0)` on both.

Let me code it up and test.

Let me think about this more carefully. The task says:

1. "every patch claimed by exactly one side" - this is the symmetric difference
2. "with adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two" - this means if there's a park that shifted, there's an area lost by one side and gained by the other. They should be merged if they touch.

But wait - I interpreted "adjacent disagreement patches" as any patches that touch (geometrically adjacent), regardless of source. But the user says "a shifted park's gain and loss form a single feature" - meaning an area lost by official (claimed only by OSM) and an area gained by official (claimed only by official) that are adjacent should be merged. That's what I did - I merged touching polygons regardless of source.

But there's a subtlety: when I merge polygons from different sources (e.g., a `parks_official` polygon touching a `parks_osm` polygon), the resulting source is `'both'`. That seems correct - if a park shifted, one area was lost and one was gained, and they're adjacent, so the merged group represents "this whole area was in dispute" - i.e., both sides claimed some part but not all.

Let me verify this is correct with the logic. Let's say park A in official shifted to location B in OSM. The official-only area is where official has park but OSM doesn't (the old location). The OSM-only area is where OSM has park but official doesn't (the new location). If these are adjacent (touching), they'd be merged into one "both" group. That makes sense - it represents the area of the shift.

Wait, but the user says "every patch claimed by exactly one side" - that implies we first find areas claimed by exactly one side (official - OSM) and (OSM - official). Then merge adjacent ones. The `both` case emerges from adjacent patches from different sources being merged. That sounds right.

Let me also double-check that I'm not missing the "claimed by both" areas. The user says "every patch claimed by exactly one side" - so I intentionally exclude areas claimed by both. The intersection (both) is implicit in the surrounding area not being disagreement. Good.

However, I realize there might be a subtle issue. If an area is claimed by both (intersection), then within the symmetric difference, we have only areas claimed by one side. When two adjacent patches from different sides merge, the merged feature has no area that's claimed by both - it's just two adjacent areas each claimed by one side. So labeling it `both` is a bit misleading... but it reflects that the merged group contains contributions from both sources.

Actually, let me think about this differently. The user says:

- "every patch claimed by exactly one side" - so the basic units are official-only and osm-only
- "adjacent disagreement patches merged together" - patches that are adjacent are merged
- "a shifted park's gain and loss form a single feature" - so if official has a park here, OSM has a park there that's shifted, the area where official has it but OSM doesn't (loss from official perspective) and the area where OSM has it but official doesn't (gain from OSM perspective) should be one feature

After merging, if a group contains only official-only patches, source = 'parks_official'. If only osm-only patches, source = 'parks_osm'. If both types, source = 'both'.

Hmm, but "both" could be misinterpreted as meaning the area is claimed by both sides, when actually the merged feature contains patches claimed by one side or the other, not both. But the user said "a shifted park's gain and loss form a single feature" - so this merged feature represents a shifted situation. And its source is "both" - which means both sources have something to say about this area. I think that's correct.

Let me also verify the geometry types more thoroughly and make sure the output is clean.

Let me check if there are any edge cases I'm missing.