geo-l2-nyc-park-symdiff

Rerun

Let me start by understanding the problem:

1. We have a GeoPackage `nyc_parks.gpkg` with layers: `parks_official`, `parks_osm`, and `both`
2. Need to find areas where only one source claims a park (disagreement)
3. Merge adjacent disagreement patches so shifted parks form single features
4. Output two GeoJSON files:
   - `parks_disagreement.geojson`: MultiPolygons with `source` attribute
   - `park_label_anchors.geojson`: Points with `source` attribute, one per group

Let me first inspect the data to understand its structure.

Let me understand the problem:

1. We have two layers in a GeoPackage: `parks_official` and `parks_osm`
2. We need to find areas where they disagree - i.e., areas claimed by exactly one side
3. Adjacent disagreement patches should be merged together
4. For each merged group, we need:
   - A `source` attribute: exactly one of `parks_official`, `parks_osm`, or `both`
   - A label-anchor point strictly inside the group's geometry

Let me think about the approach:

- "Areas claimed by exactly one side" = symmetric difference between the two layers
- But also, areas claimed by both = intersection
- The user wants: every patch claimed by exactly one side, with adjacent disagreement patches merged together. So a shifted park creates a "gain" (new area OSM has that official doesn't) and a "loss" (area official has that OSM doesn't). These are adjacent and should be merged.

Wait, let me 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."

So the approach:
1. Compute the symmetric difference of the two layers (areas in one but not the other)
2. But we need to categorize each piece: is it from official-only or osm-only?
3. "Adjacent disagreement patches" - if an area from official-only touches an area from osm-only, they should be merged into a single feature with source "both"? Wait no...

Actually, let me think about this differently. The user says:
- "every patch claimed by exactly one side" - so we're looking at the symmetric difference
- "adjacent disagreement patches merged together" - so if an official-only patch touches an osm-only patch, merge them
- source values: parks_official, parks_osm, or both

So the source attribute depends on who claims it. If it's from official-only, source = "parks_official". If OSM-only, source = "parks_osm". If it's a merged group containing both official-only and osm-only parts that are adjacent, source = "both".

Wait, but they also say "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 an official park and an OSM park cover similar but offset areas, the official-only part (loss) and OSM-only part (gain) should be merged together. Since that merged group contains both official-only and osm-only areas, the source would be "both".

So the final three sources are:
- `parks_official` - areas only in official, not touching any OSM-only area
- `parks_osm` - areas only in OSM, not touching any official-only area
- `both` - merged groups that contain both official-only and OSM-only components that are adjacent

Also areas that are both (intersection) - actually no, the user said "every patch claimed by exactly one side". So the intersection (areas in both) is NOT a disagreement patch. The disagreement patches are only the symmetric difference areas.

But then, if I have two polygonal differences that are adjacent (one official-only, one osm-only), I merge them and source is "both".

Let me think about the algorithm:

1. Compute: `official_only = official.difference(osm)` → multipolygons of areas only in official
2. Compute: `osm_only = osm.difference(official)` → multipolygons of areas only in OSM
3. Dissolve/union all these single-part polygons together
4. But we need to track which source(s) each part comes from
5. Merge adjacent parts regardless of source
6. For the final merged multipolygon, determine if it contains parts from official_only, osm_only, or both → set source accordingly

Actually, let me think more carefully. The user wants the "symmetric difference" areas, but then merge adjacent disagreement patches.

Step 1: Compute the official-only areas (dissolved/union of official minus osm)
Step 2: Compute the osm-only areas (dissolved/union of osm minus official)
Step 3: Explode all these into individual polygons
Step 4: For each polygon, tag it as `parks_official` or `parks_osm`
Step 5: Merge adjacent polygons (regardless of their tag)
Step 6: For each merged group, if it contains only parks_official polygons, source = "parks_official"; if only parks_osm, source = "parks_osm"; if both, source = "both"
Step 7: Create a label anchor point inside each merged multipolygon

For the label anchor, I need a point strictly inside the geometry. I can use `representative_point()` from shapely which guarantees it's within the geometry (though it might be near the boundary). Or I could use `point_on_surface()`. Let me use `representative_point()` which is guaranteed to be within.

Let me write the Python script.

Actually wait - I need to be careful with the geometry operations. Let me plan:

1. Read the two layers from the GPKG
2. Dissolve both layers (union all geometries within each layer)
3. Compute symmetric difference parts
4. Actually, I should work with individual polygons from both sides

Let me think about this differently. The goal is:
- Find all areas that are in one layer but not the other
- Merge adjacent such areas

A simpler approach:
1. Explode both layers into individual polygons
2. Compute `official_only_polys = official_polygons - osm_union` 
3. Compute `osm_only_polys = osm_polygons - official_union`
4. Tag each resulting polygon with its source
5. Merge adjacent/touching polygons together (regardless of source)
6. Determine source of each merged group

But actually, since we're dealing with multipolygons, let me use a different approach:

1. Create a dissolved/unioned geometry for official and osm
2. Compute:
   - `official_only = official_union.difference(osm_union)`
   - `osm_only = osm_union.difference(official_union)`
3. Explode these into individual polygons
4. Tag each polygon
5. Merge touching polygons from different sources
6. Set source based on composition

For merging adjacent polygons, I can use `unary_union` with different source tags... but that merges everything. I need a more careful approach.

One way: 
1. Collect all polygons from official_only and osm_only
2. Use a spatial join / adjacency check to find which polygons touch
3. Group connected components
4. Merge within each group

For the connected components approach, I can:
- Build a graph where nodes are polygons and edges indicate touching
- Find connected components
- Union the polygons in each component

Let me write the script:

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

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

# Dissolve each layer
official_union = unary_union(official.geometry)
osm_union = unary_union(osm.geometry)

# Compute symmetric difference parts
official_only = official_union.difference(osm_union)  # areas only in official
osm_only = osm_union.difference(official_union)       # areas only in OSM

# Function to extract individual polygons from a geometry
def extract_polygons(geom, source):
    polys = []
    if geom is None or geom.is_empty:
        return polys
    if geom.geom_type == 'Polygon':
        polys.append((geom, source))
    elif geom.geom_type == 'MultiPolygon':
        for g in geom.geoms:
            polys.append((g, source))
    return polys

# Extract and tag polygons
polygons_with_source = []
polygons_with_source.extend(extract_polygons(official_only, 'parks_official'))
polygons_with_source.extend(extract_polygons(osm_only, 'parks_osm'))

print(f"Total polygons: {len(polygons_with_source)}")

# Build adjacency graph
n = len(polygons_with_source)
adj = {i: [] for i in range(n)}

for i in range(n):
    for j in range(i+1, n):
        if polygons_with_source[i][0].touches(polygons_with_source[j][0]):
            adj[i].append(j)
            adj[j].append(i)

# Find connected components
visited = set()
components = []
for i in range(n):
    if i not in visited:
        stack = [i]
        comp = set()
        while stack:
            node = stack.pop()
            if node not in visited:
                visited.add(node)
                comp.add(node)
                for nb in adj[node]:
                    if nb not in visited:
                        stack.append(nb)
        components.append(comp)

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

# Merge polygons in each component
results = []
for comp in components:
    comp_polys = [polygons_with_source[i][0] for i in comp]
    comp_sources = set(polygons_with_source[i][1] for i in comp)
    
    # Merge
    merged = unary_union(comp_polys)
    
    # Determine source
    if len(comp_sources) == 1:
        source = comp_sources.pop()
    else:
        source = 'both'
    
    # Get label anchor
    if merged.geom_type == 'MultiPolygon':
        # Use representative_point of the largest polygon
        largest = max(merged.geoms, key=lambda g: g.area)
        anchor = largest.representative_point()
    else:
        anchor = merged.representative_point()
    
    results.append({
        'geometry': merged,
        'source': source,
        'anchor': anchor
    })

# Write to GeoJSON - disagreement polygons
features_disagreement = []
features_anchors = []
for i, r in enumerate(results):
    feat_disagreement = {
        'type': 'Feature',
        'properties': {'source': r['source']},
        'geometry': r['geometry'].__geo_interface__
    }
    features_disagreement.append(feat_disagreement)
    
    feat_anchor = {
        'type': 'Feature',
        'properties': {'source': r['source']},
        'geometry': r['anchor'].__geo_interface__
    }
    features_anchors.append(feat_anchor)

fc_disagreement = {
    'type': 'FeatureCollection',
    'features': features_disagreement
}

fc_anchors = {
    'type': 'FeatureCollection',
    'features': features_anchors
}

# Save to files
with open('/work/parks_disagreement.geojson', 'w') as f:
    json.dump(fc_disagreement, f)

with open('/work/park_label_anchors.geojson', 'w') as f:
    json.dump(fc_anchors, f)

print("Done!")
```

Wait, actually I need to be more careful. The `representative_point()` method returns a point that is guaranteed to be within the geometry. Let me use that.

Also, I need to be careful about the coordinate system - it's in feet (NY Long Island ftUS), but GeoJSON technically should be in WGS84 (EPSG:4326). However, the user didn't specify a CRS transformation, and they're asking for GeoJSON. Let me check - GeoJSON spec says coordinates should be in WGS84. But some applications accept projected CRS too. Let me keep them in the original CRS since the user didn't ask for reprojection... Actually, GeoJSON RFC 7946 requires WGS84. But many tools accept any. I'll keep them as-is since the user didn't mention reprojection, and all the analysis is in the projected CRS.

Hmm, actually let me re-read the prompt. No mention of CRS transformation. I'll leave as-is in the projected CRS coordinates.

Let me also consider: what if there are no disagreement areas (all parks match perfectly)? That's unlikely but I should handle it gracefully.

Let me write and run the script.

Let me verify the requirements more carefully:

1. "every patch claimed by exactly one side" - I used symmetric difference (official_only + osm_only). That's correct for areas claimed by exactly one side.

2. "with adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two" - I built an adjacency graph using `.touches()` and merged connected components. But wait - I'm merging polygons that are from *different* source sides (official_only and osm_only) into "both" groups. Let me re-read the requirement.

The user wants: patches claimed by exactly one side (official_only or osm_only), then merge adjacent disagreement patches together. So if a park shifted, the area where official claims it but OSM doesn't (official_only) touches the area where OSM claims it but official doesn't (osm_only) - these two patches should be merged into one feature with source "both".

My current approach does exactly this. For groups containing patches from both sources, source = "both". For groups from only one source, source is that source. Good.

3. "One MultiPolygon per merged group" - I used unary_union which can produce a Polygon or MultiPolygon. Actually the requirement says "One MultiPolygon per merged group" - so even if the result is a single Polygon, I need to wrap it in a MultiPolygon to be safe. Let me check.

Actually, looking more carefully at the requirement: "One MultiPolygon per merged group" - I think they want the output geometry to always be a MultiPolygon. Let me fix that to ensure consistency.

4. "label-anchor point that falls strictly inside the group's geometry" - I used `representative_point()` which is guaranteed to be inside the geometry. Good.

5. Both files need `source` column (as a property in GeoJSON).

6. Anchors in matching order - since I build both lists in the same component loop, they're in matching order.

Let me also fix the MultiPolygon requirement and double-check that `representative_point()` is strictly inside. Per Shapely docs, `representative_point()` returns a point guaranteed to be within the object (i.e., on the boundary or interior). Actually, let me check: shapely's `representative_point()` says "Returns a point that is guaranteed to be within the object, cheaply." It should be strictly inside (not on boundary). But to be safe, let's use `.point_on_surface()` instead which is documented to be "a point guaranteed to be within the object" - actually those are the same thing. Let me use `.representative_point()`.

Actually, the user says "strictly inside the group's geometry" - meaning interior, not boundary. `representative_point()` does guarantee it's within the object (including boundary). Hmm, actually in shapely, `representative_point()` uses the point_on_surface algorithm which does return a point in the interior, not on the boundary. Let me verify with a quick test.

Actually, I'll just ensure the output is correct and also wrap polygons in multipolygons as requested.

Let me update the script to ensure MultiPolygon output.

Everything checks out:

1. Disagreement features are all MultiPolygon ✓
2. Anchor features are all Point ✓
3. Sources match between the two files in order ✓
4. All anchor points are strictly inside their group geometry ✓
5. The source values are exactly `parks_official`, `parks_osm`, or `both` ✓

Let me do a final check that the files exist and look reasonable.