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 that are claimed by exactly one side (disagreement areas)
3. Adjacent disagreement patches should be merged together
4. Output two GeoJSON FeatureCollections:
   - `parks_disagreement.geojson` with a `source` attribute (`parks_official`, `parks_osm`, or `both`)
   - `park_label_anchors.geojson` with same `source` attribute, one anchor per disagreement group

Let me first look at what's in the GeoPackage.

Let me understand the data:

- `parks_official`: 1380 features, CRS=EPSG:6539 (NY State Plane Long Island, feet)
- `parks_osm`: 1372 features, CRS=EPSG:6539

Both have a `park_id` column. Let me check if the park_ids align.

The task:
1. Find areas claimed by exactly one side (symmetric difference).
2. Merge adjacent disagreement patches together so that a shifted park's gain and loss form a single feature.
3. One MultiPolygon per merged group with a `source` attribute: `parks_official`, `parks_osm`, or `both`.
4. For each group, a label-anchor point strictly inside the group's geometry.
5. Output: `parks_disagreement.geojson` and `park_label_anchors.geojson`.

Let me think about the approach:

1. First, dissolve each layer into a single multipolygon (all parks merged).
2. Compute the symmetric difference = (official U osm) - (official ∩ osm).
   - Actually, we want areas claimed by exactly one side. So:
   - `only_official = official - osm` (areas in official but not OSM)
   - `only_osm = osm - official` (areas in OSM but not official)
3. These are the "disagreement" patches. But we need to merge adjacent disagreement patches.
   - Wait, what does "adjacent disagreement patches" mean? If a park is shifted, the gain and loss are separate polygons. They may be adjacent (touching) and should be merged.
   - So we need to take all the symmetric difference polygons, group them by connectivity (touching), and merge each group into a single MultiPolygon.
   - Each group gets a `source` attribute based on what kind of polygons are in it. But actually, each polygon in the symmetric difference is either `only_official` or `only_osm`. If we merge adjacent ones, a merged group could contain both types. But wait - can an `only_official` polygon touch an `only_osm` polygon? In theory yes if the park shifted.
   
   Actually, let me re-think. "Every patch claimed by exactly one side" - this is the symmetric difference. Then "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 was shifted, there's an area in official that's not in OSM (loss) and an area in OSM that's not in official (gain). These two patches might be adjacent (touching at a vertex or edge). If they're adjacent, merge them into one.

   For grouping: use unary_union on the symmetric difference, then... actually we need to group by touching/adjacent relationship.

   Approach:
   - Compute `only_official = official.difference(osm)`
   - Compute `only_osm = osm.difference(official)`
   - Combine them into a single GeoSeries
   - Dissolve/merge touching polygons into groups
   - For each group, create a MultiPolygon
   - Assign `source` based on... hmm. If a group contains only `only_official` parts, source is `parks_official`. If only `only_osm`, source is `parks_osm`. If both, what? 

   Actually, re-reading: "a source attribute whose value is exactly one of parks_official, parks_osm, or both". So if a merged group contains only polygons from official, it's `parks_official`. If only from OSM, it's `parks_osm`. If mixed (an adjacent gain and loss merged together), it would be... well, that doesn't directly make sense. Let me think more carefully.

   Actually maybe "both" means the area is in both sources (intersection). But the task says "every patch claimed by exactly one side" - that's the symmetric difference. And then "adjacent disagreement patches merged together" - so we merge the symmetric difference polygons. The source would then be... 

   Hmm, 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 (those exact strings, matching the input GPKG layer names)."

   So `source` can be `parks_official`, `parks_osm`, or `both`. If a group came entirely from the official side (i.e., area in official but not OSM), source = `parks_official`. If entirely from OSM side, source = `parks_osm`. If the group has a mix (gain + loss merged together), source = `both`.

   Wait but the prompt says "every patch claimed by exactly one side" - that's symmetric difference, not intersection. So `both` wouldn't apply to symmetric difference... unless merged groups can contain both types.

   I think the intent is:
   1. Start with symmetric difference polygons (each is either `only_official` or `only_osm`).
   2. Merge adjacent/touching polygons into groups.
   3. For each group:
       - If all polygons are `only_official` → source = `parks_official`
       - If all polygons are `only_osm` → source = `parks_osm`
       - If mixed (gain+loss adjacent) → source = `both`

   This makes sense! A shifted park creates both a gain (OSM side) and a loss (official side) that are adjacent, and when merged they represent the full "disagreement zone" for that shifted park. For that merged group, `both` is appropriate.

   Let me code this up.

   For the label anchor: a point strictly inside the group's geometry. I can use `representative_point()` which is guaranteed to be inside the polygon.

   Let me code the solution.

```python
import geopandas as gpd
from shapely.ops import unary_union, polygonize
from shapely.geometry import MultiPolygon, Polygon
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')

# Reproject to a suitable projected CRS for area operations (already in EPSG:6539)

# Dissolve each into a single geometry
official_union = unary_union(official.geometry)
osm_union = unary_union(osm.geometry)

# Symmetric difference components
official_only = official_union.difference(osm_union)  # in official, not in OSM
osm_only = osm_union.difference(official_union)  # in OSM, not in official

# We need to handle these as collections of polygons
# Use polygonize or just iterate over geoms

def get_polygons(geom):
    """Extract individual polygons from a geometry."""
    if geom is None or geom.is_empty:
        return []
    if geom.geom_type == 'Polygon':
        return [geom]
    elif geom.geom_type == 'MultiPolygon':
        return list(geom.geoms)
    else:
        return []

official_only_polys = get_polygons(official_only)
osm_only_polys = get_polygons(osm_only)

# Tag them
tags = []
for p in official_only_polys:
    tags.append({'geometry': p, 'source_side': 'parks_official'})
for p in osm_only_polys:
    tags.append({'geometry': p, 'source_side': 'parks_osm'})

# Merge adjacent/touching polygons
# Use unary_union which merges touching polygons
# But we need to track source_side across merged groups

# Build a GeoDataFrame
gdf = gpd.GeoDataFrame(tags, crs=official.crs)

# Use spatial join to find touching relationships
# Perform a dissolve by connectivity using unary_union
# The trick: buffer(0) to snap, then use polygonize + spatial join

# Alternative: use a graph approach - for each polygon, find all polygons it touches
# Then find connected components

from shapely.ops import unary_union

# We can merge all polygons (across both sides) with a small negative buffer to snap 
# Then dissolve them by connectivity

# Actually simpler: use unary_union on all polygons, and for each resulting polygon
# figure out which source_side polygons contributed to it

# Step 1: Merge all polygons
all_polys = gdf.geometry.tolist()
merged = unary_union(all_polys)

# If merged is a Polygon, make it a list
merged_polys = get_polygons(merged)

# Now for each merged polygon, determine which original tagged polygons intersect/touch it
groups = []
for mp in merged_polys:
    # Find all original polygons that intersect this merged polygon
    intersecting = gdf[gdf.intersects(mp)]
    sources = intersecting['source_side'].unique()
    
    if len(sources) == 1:
        source = sources[0]
    else:
        source = 'both'
    
    groups.append({'geometry': mp, 'source': source})

# Create disagreement GeoDataFrame
disagreement = gpd.GeoDataFrame(groups, crs=official.crs)
disagreement = disagreement[['source', 'geometry']]

# Create label anchors
anchors = disagreement.copy()
anchors['geometry'] = anchors['geometry'].representative_point()

# Save
disagreement.to_file('/work/parks_disagreement.geojson', driver='GeoJSON')
anchors.to_file('/work/park_label_anchors.geojson', driver='GeoJSON')
```

Wait, there's a subtle issue with using `unary_union` on all polygons - it would merge touching polygons regardless of source_side. That's exactly what we want! Adjacent polygons (whether same side or different side) get merged. Then for each merged polygon, we check which source_side polygons contribute to it.

But actually, there could be a case where two `parks_official` polygons are adjacent but not touching an `osm_only` polygon - they should be merged together into one group. That's fine with unary_union.

Also, I need to think about whether `unary_union` correctly merges polygons that only touch at a point... In shapely, unary_union by default only merges polygons that overlap or share edges, not those that only touch at a point. But the task says "adjacent" and "touching" - let me check if touching at a vertex counts.

"adjacent disagreement patches merged together" - I think edge-adjacent is what's meant. unary_union handles that fine.

But wait, there's another subtlety. Two polygons from the same source could be in different "groups" if they don't touch each other. That's fine - they'd be separate features in the output.

Let me also handle the case where `unary_union` returns a GeometryCollection or other types.

Let me also consider: what if a polygon from official_only and a polygon from osm_only are adjacent but at a single point? In standard topology, they wouldn't merge. But if the task intends "touching at a point" to also merge... let me check if shapely's unary_union merges polygons that touch at a point. I believe it does NOT by default. The polygons would remain separate but in the same group. That seems fine.

Actually, hmm, let me reconsider. If two polygons only touch at a point, unary_union won't merge them into a single polygon. They'd stay as separate polygons in a MultiPolygon. But I said `merged = unary_union(all_polys)` - if they only touch at a point, they'd be a MultiPolygon with two polygons in it. And then I'd iterate over the individual polygons. So they'd end up as separate groups.

But I think "adjacent" means sharing an edge. Let me go with unary_union as is.

Wait, actually I realize my approach of iterating over `unary_union` result's polygons is wrong. If unary_union merges touching polygons into one polygon, great. But if it doesn't merge them (e.g., point-touching), the resulting collection would have them as separate polygons. I'd lose the connection.

Let me think of a better approach.

Better approach using graph connected components:
1. Build a graph where each polygon is a node.
2. Add an edge between nodes if their polygons touch/intersect.
3. Find connected components.
4. For each component, merge all polygons into a single MultiPolygon.
5. Determine source based on the unique sources in the component.

This handles both edge-touching and point-touching adjacency.

```python
import networkx as nx
from shapely import strtree

# Build STRtree for spatial indexing
tree = strtree.STRtree(gdf.geometry.tolist())

G = nx.Graph()
G.add_nodes_from(range(len(gdf)))

for i, geom in enumerate(gdf.geometry):
    # Find potential neighbors
    candidates = tree.query(geom)
    for j in candidates:
        if i != j and geom.touches(gdf.geometry.iloc[j]):
            G.add_edge(i, j)

# Find connected components
components = list(nx.connected_components(G))

groups = []
for comp in components:
    indices = list(comp)
    polys = [gdf.geometry.iloc[i] for i in indices]
    sources = [gdf['source_side'].iloc[i] for i in indices]
    
    merged = unary_union(polys)
    if merged.geom_type == 'Polygon':
        merged = MultiPolygon([merged])
    elif merged.geom_type == 'MultiPolygon':
        pass
    else:
        # GeometryCollection etc - extract polygons
        merged = MultiPolygon([g for g in merged.geoms if g.geom_type == 'Polygon'])
    
    unique_sources = set(sources)
    if len(unique_sources) == 1:
        source = list(unique_sources)[0]
    else:
        source = 'both'
    
    groups.append({'geometry': merged, 'source': source})
```

This is cleaner. Let me use this approach.

Actually, I don't need networkx - I can implement a simple BFS/DFS on the adjacency graph using the STRtree. But networkx is cleaner and probably installed. Let me check.

Actually, let me keep it simple and not rely on networkx. I can use a BFS on the graph.

Let me also be careful about `touches` vs `intersects`. The task says "adjacent" - I think `touches` (share boundary but don't overlap) is the right relationship, because symmetric difference polygons can't overlap (they're from the difference operation).

Wait, actually, can two symmetric difference polygons overlap? If official_only and osm_only are both computed as differences against the other, then they shouldn't overlap by definition. And within each source side, the polygons come from a unary_union difference, so they won't overlap either. So `touches` is the right check.

But what if a park shifted such that the gain polygon touches the loss polygon at multiple points along an edge? That's exactly `touches`.

Let me also think about using `buffer(0)` to clean up any geometry issues before processing.

Let me write the full script now.

Oh wait, another consideration: `official_only` and `osm_only` from `difference()` might have slivers or artifacts. I should clean them up, perhaps by removing very small polygons.

Let me also check if `networkx` is available.

Actually, I'll just implement a simple connected components without networkx.

Let me write the full script.

Let me think about the problem:

1. We have two layers: `parks_official` and `parks_osm` in EPSG:6539 (NY state plane, feet).
2. We need to find areas where only one side has parks (disagreement), merge adjacent disagreement patches, and assign a `source` attribute.
3. Also need label-anchor points inside each merged group.

The approach:
- Dissolve all parks in each layer into a single multi-polygon (or union them).
- Compute the symmetric difference (areas in official but not in OSM, and areas in OSM but not in official).
- Wait, but there might be cases where both layers have geometries but they overlap or don't overlap. Let me think more carefully.

Actually, looking at the data, they both have the same `park_id` values. So each park is identified by `park_id`. The question is about which geometries are claimed by one dataset vs the other. So for a given park_id, they might have different geometries.

Let me reconsider. The instruction says "every patch claimed by exactly one side". So:

1. Union/dissolve all official parks into one geometry.
2. Union/dissolve all OSM parks into one geometry.
3. The symmetric difference between these two unioned geometries gives us the areas claimed by exactly one side.
4. But we need to split the symmetric difference into:
   - Areas that are in official but not in OSM → source = 'parks_official'
   - Areas that are in OSM but not in official → source = 'parks_osm'

Then "merge adjacent disagreement patches" - if a park's boundary shifted, the gain on one side and loss on the other might be adjacent. So we need to take all the disagreement polygons (both official-only and OSM-only) and merge any that are adjacent (touch) into a single MultiPolygon per group.

For each merged group, we need a label anchor point that falls strictly inside.

Let me code this up.

The steps:
1. Load both layers, dissolve each into a single geometry (union).
2. Compute `official_only = official.difference(osm)` - areas only in official
3. Compute `osm_only = osm.difference(official)` - areas only in OSM
4. Combine these into a single GeoDataFrame with a `source` column.
5. For "adjacent disagreement patches merged together": use spatial graph clustering - any polygons that touch (including along edges) should be merged into one group.
6. For each group, create a MultiPolygon of all constituent polygons, and compute a label anchor point that is strictly inside.

For the label anchor, I can use `representative_point()` from shapely which guarantees a point within the geometry.

Let me write the script.

Actually, let me think about this more carefully. When we dissolve the official parks, we might create a massive multi-polygon. Then difference with OSM. This should work.

Wait - but what about areas where both have parks? The symmetric difference won't include those. And the areas where both agree don't need to be in the output.

Let me also think about performance - 1380 and 1372 features, so the dissolve should be fine.

Let me code this up:

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

# Load the layers
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 each into a single geometry (union)
official_union = unary_union(official.geometry)
osm_union = unary_union(osm.geometry)

# Compute symmetric difference parts
official_only_geom = official_union.difference(osm_union)
osm_only_geom = osm_union.difference(official_union)

# Also compute the areas both agree
both_geom = official_union.intersection(osm_union)
```

Hmm, but the `source = 'both'` was mentioned in the prompt. 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. One MultiPolygon per merged group with a `source` attribute whose value is exactly one of `parks_official`, `parks_osm`, or `both`"

Wait, so there's a `source` of `both` too? That seems odd given "every patch claimed by exactly one side". Let me re-read.

Actually reading again: "I need every patch claimed by exactly one side" - so it's the disagreement areas only. But then the source attribute can be one of three values including 'both'. 

Hmm, maybe the 'both' is for when after merging adjacent patches, some patches that were on the boundary between official-only and osm-only get merged together? No, that doesn't make sense.

Wait, let me read the exact wording more carefully:

"One MultiPolygon per merged group with a `source` attribute whose value is exactly one of `parks_official`, `parks_osm`, or `both`"

So the source value is one of those three strings. For disagreement patches, it would be either `parks_official` (only in official) or `parks_osm` (only in OSM). But what about `both`?

Maybe the `both` source is for areas where both agree? But the prompt says "every patch claimed by exactly one side". 

Or maybe after merging adjacent disagreement patches, a merged group could contain both official-only and osm-only parts that are adjacent? In that case, the merged group would have source 'both'? That seems a bit odd but let me re-read:

"adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two"

Yes! This confirms it. If a park shifted, the area gained by the official side and the area lost (now only in OSM) might be adjacent. So these get merged into one group with source 'both'.

So the logic is:
1. Compute official_only and osm_only polygons
2. Build a graph of adjacency - connect any polygons that touch
3. Each connected component becomes a group
4. For each group, if it contains official_only polygons, it's marked 'parks_official'. If only osm_only, it's 'parks_osm'. If both, it's 'both'.

Actually wait, that doesn't make sense either. Let me re-read one more time.

"I need every patch claimed by exactly one side" - so it's the disagreement areas.

"adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two" - so if an official-only patch is adjacent to an OSM-only patch (because a park shifted), they get merged.

"One MultiPolygon per merged group with a `source` attribute whose value is exactly one of `parks_official`, `parks_osm`, or `both`"

So for each merged group, we assign one of three source values. If the group contains only official-only polygons → 'parks_official'. If only OSM-only → 'parks_osm'. If it contains both types (because they're adjacent, e.g., a shifted park) → 'both'.

That makes perfect sense! Let me implement this.

Steps:
1. Union-dissolve each layer.
2. Compute official_only = official_union - osm_union (difference)
3. Compute osm_only = osm_union - official_union (difference)
4. Extract individual polygons from each result (handling MultiPolygon)
5. Assign source type to each polygon
6. Build adjacency graph where polygons touch (using DE-9IM: touches or intersects? For merging "adjacent", touching should be enough, but a park shift could also overlap slightly. Let me use touches + overlaps/within for robustness... Actually let me just use intersects since we're dealing with disjoint patches anyway after the difference operation. Wait, difference guarantees no overlap. So touching is sufficient for adjacency.)

Hmm actually, when we do `official_union.difference(osm_union)`, the result won't overlap with `osm_only`. But two patches within `official_only` could be adjacent. And an `official_only` patch could be adjacent to an `osm_only` patch. So we need to group all touching polygons regardless of type.

7. For each connected component (graph cluster), collect all polygons, merge into a MultiPolygon, and determine the source.
8. For each group, compute a representative point inside.

Let me write the script.

Actually, I need to think about what "strictly inside" means for the anchor point. Shapely's `representative_point()` guarantees a point within the geometry (not on boundary). That should work.

Let me code this up:

```python
import geopandas as gpd
import networkx as nx
from shapely.geometry import MultiPolygon, Polygon, shape
from shapely.ops import unary_union
from shapely import wkt

# Load
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

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

# Disagreement areas
official_only = official_union.difference(osm_union)
osm_only = osm_union.difference(official_union)

# Helper 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 p in geom.geoms:
            polys.append((p, source))
    return polys

polygons = []
polygons.extend(extract_polygons(official_only, 'parks_official'))
polygons.extend(extract_polygons(osm_only, 'parks_osm'))

# Build adjacency graph
G = nx.Graph()
for i, (poly, src) in enumerate(polygons):
    G.add_node(i, geometry=poly, source=src)

# Add edges for touching polygons
for i in range(len(polygons)):
    for j in range(i+1, len(polygons)):
        if polygons[i][0].touches(polygons[j][0]):
            G.add_edge(i, j)

# Find connected components
components = list(nx.connected_components(G))

# Build output features
groups = []
for comp in components:
    group_polys = []
    sources = set()
    for idx in comp:
        group_polys.append(polygons[idx][0])
        sources.add(polygons[idx][1])
    
    merged = MultiPolygon(group_polys)
    
    if len(sources) == 2:
        source = 'both'
    else:
        source = list(sources)[0]
    
    groups.append({'geometry': merged, 'source': source})

# Create GeoDataFrames
gdf_disagreement = gpd.GeoDataFrame(groups, crs=crs, geometry='geometry')

# Label anchors - use representative_point
gdf_anchors = gdf_disagreement.copy()
gdf_anchors['geometry'] = gdf_anchors['geometry'].representative_point()

# Write GeoJSON
gdf_disagreement.to_file('/work/parks_disagreement.geojson', driver='GeoJSON')
gdf_anchors.to_file('/work/park_label_anchors.geojson', driver='GeoJSON')
```

This looks good. Let me write it and run it.

But wait, I should think about performance. `official_union.difference(osm_union)` could be slow on complex geometries but with ~1300 features it should be fine.

Also, I should be careful about the graph - if there are 1000+ individual polygons, the nested loop for adjacency testing is O(n^2). We can use spatial indexing to speed it up. Let me use geopandas' sjoin for this.

Let me use a more efficient approach with an R-tree.

Actually, let me first check how many individual polygons we get from the difference operation.

Let me just write and run the script. If it's too slow I'll optimize.

Let me understand the problem:

1. We have two layers: `parks_official` and `parks_osm`
2. We need to find areas of disagreement - areas claimed by exactly one side
3. Adjacent disagreement patches should be merged together (so if a park shifted, its gain and loss merge into one feature)
4. Each merged group gets a `source` attribute: `parks_official`, `parks_osm`, or `both`
5. Need label-anchor points strictly inside each group

Let me think about the approach:

The user says "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 we need Symmetric Difference of the two layers - areas that are in official OR OSM but not both. And then merge adjacent/touching parts.

Wait, but also areas claimed by both should be included if I read it right. Let me re-read:

"I need every patch claimed by exactly one side"

So areas where both layers agree (intersection) are not needed. We want the symmetric difference.

"with adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two."

So if official has A and OSM has B, where A and B are overlapping but shifted, then:
- A - B is a gain/loss area (official only)
- B - A is a gain/loss area (osm only)

But if these two patches are adjacent (share a boundary), they should be merged into one MultiPolygon.

Wait, but they have different sources! How can a single merged feature have two different sources?

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

Hmm, so the merged groups can have source `both`? That doesn't make sense if we're only looking at disagreement patches claimed by exactly one side.

Wait, let me re-read again: "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"

I think the key insight is that when a park boundary shifts:
- The area where the old boundary was but not the new one is claimed by official only
- The area where the new boundary is but not the old one is claimed by OSM only

These two areas share a boundary (they're adjacent). So when we merge adjacent disagreement patches, a shifted park's gain (OSM-only) and loss (Official-only) become one merged feature.

But what source would such a merged feature have? It would be a mix of both... So the source would be... hmm.

Actually, I think the interpretation is different. Let me think about this as a GIS operation:

1. Compute the symmetric difference between the two layers (union minus intersection).
2. The symmetric difference gives us polygons that are in exactly one of the layers.
3. But some of these polygons might be adjacent (share a boundary). Merge those.
4. For each merged polygon, determine whether it came from official, osm, or both.

Wait, but if we merge an official-only polygon with an adjacent OSM-only polygon, the merged result contains areas from both sources. So the source would be... both?

Hmm, but "every patch claimed by exactly one side" - a merged feature wouldn't be claimed by exactly one side.

Let me reconsider. Perhaps the answer is simpler:

We take the symmetric difference (areas claimed by exactly one of the two datasets). These are polygons that are in official XOR osm. Then we merge any adjacent/touching polygons within this set. Each resulting merged polygon could contain pieces that were original from official (official-only portions) or from osm (osm-only portions). For a merged polygon that contains both types of pieces, the source would be... both?

Actually, looking at "parks_official, parks_osm, or both" - these are the layer names. So the source indicates which original layer(s) contributed to this merged group.

Let me think about this differently. Let's think about it as an overlay operation:

1. Union of both layers
2. Assign each union polygon a "source" based on which layers contain it
3. Filter to polygons where source is only one layer (not both)
4. Dissolve adjacent polygons with the same source
5. But also... "adjacent disagreement patches merged" - does this mean adjacent patches with DIFFERENT sources get merged too?

"a shifted park's gain and loss form a single feature rather than two" - yes, gain and loss have different sources (official only vs osm only), but they're adjacent along the shifted boundary, so they should merge.

So: 
- Take symmetric difference (official - osm) and (osm - official)
- Union these together
- Dissolve (merge) all touching/adjacent polygons regardless of source
- For each merged group, determine its source: if it contains areas from both original layers, source = "both"; if only from official, source = "parks_official"; if only from osm, source = "parks_osm"

Wait, actually let me re-read: "the `source` attribute whose value is exactly one of `parks_official`, `parks_osm`, or `both`"

And "every patch claimed by exactly one side"

I think the source "both" would mean the merged group contains areas from both official and osm (i.e., it's the merged gain+loss of a shifted park). The source "parks_official" means the area is only from official (no corresponding OSM park at all). Similarly for "parks_osm".

Actually, I think there's another way to read this:

Maybe the user wants the full picture - not just the XOR. They want:
- Areas claimed by official only → source = parks_official
- Areas claimed by OSM only → source = parks_osm
- Areas claimed by both → source = both

But then merging adjacent disagreement patches... 

Hmm, let me think again. "every patch claimed by exactly one side" - this sounds like XOR. But "adjacent disagreement patches merged" and "a shifted park's gain and loss form a single feature rather than two" - this echoes the XOR interpretation where the gain (OSM-only) and loss (Official-only) of a shifted park should merge into one.

But then the source of that merged feature would be ambiguous... unless we label it "both".

Actually, I think the most useful interpretation is:

1. Compute the symmetric difference (XOR) = (Official ∪ OSM) - (Official ∩ OSM)
2. Merge all adjacent/touching polygons in the symmetric difference
3. For each merged group, determine its source based on what original areas it contains:
   - If it contains only Official-only area → "parks_official"
   - If it contains only OSM-only area → "parks_osm"
   - If it contains both → "both"

This makes sense because when a park shifts, the old location (Official-only) and new location (OSM-only) would be adjacent and merge, and the source would be "both" indicating the merged group contains disagreement from both sides.

Let me code this up.

Actually, another thought: maybe I should just compute the full overlay:
1. Compute Official - OSM → official_only
2. Compute OSM - Official → osm_only
3. Compute Official ∩ OSM → both
4. Merge all adjacent features from all three groups together
5. For each merged group, determine which of the three categories contributed

But "every patch claimed by exactly one side" - this excludes "both" areas from the disagreement set.

Hmm, let me just go with my initial interpretation since it seems the most coherent:

1. Compute symmetric difference (Official XOR OSM)
2. Dissolve adjacent polygons
3. Tag each merged polygon with source
4. Generate label anchors

Let me code this.

Actually wait - the user says "with adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two."

A shifted park: Official has park A at location X, OSM has park A at location Y (shifted).
- Official-only area: X - Y (where the park was officially but OSM doesn't have it)
- OSM-only area: Y - X (where OSM has it but official doesn't)

These areas are adjacent (they share a boundary along the overlap). Merging them yields a single feature covering both X-Y and Y-X, which encompasses the whole shifted area.

The source of this merged feature would be "both" because it contains areas from both official and osm.

For a park that exists ONLY in official (no OSM counterpart), the official-only area wouldn't be adjacent to any OSM-only area, so it stays as its own feature with source "parks_official".

OK let me implement.

Also, I need to check the CRS - it's EPSG:6539 (NY state plane, feet). For a "label anchor point strictly inside the group's geometry" I'll use `representative_point()` which guarantees being inside.

Let me write the script.

Wait, the dissolve (merge adjacent polygons) is key. In shapely/geopandas, I can do:

1. Compute official_only = official.difference(osm.unary_union) ... no, this is per park_id. 

Actually, I need to think about this differently. The two layers don't necessarily have matching park_ids. Let me look at the data more carefully.

Actually, from the output above, most park_ids are shared (1360 out of 1380/1372). So there's some alignment.

Let me think about the right approach:

Option A: Park-by-park approach
- For each park_id that exists in both, compute symmetric difference
- For parks only in official, those are official-only
- For parks only in OSM, those are OSM-only
- Then merge adjacent/all features

But this doesn't handle cases where parks don't match by ID.

Option B: Full overlay approach
- Union both layers
- For each piece of the union, determine if it's in official, osm, or both
- Filter to pieces that are in exactly one
- Merge adjacent pieces

Let me go with Option B as it's cleaner.

Actually, let me think more carefully. The full overlay (union of both layers) produces many small slivers. Let me think about an alternative.

**Symmetric difference approach:**
1. Compute the union of all official geometries → `official_union`
2. Compute the union of all OSM geometries → `osm_union`
3. `official_only = official_union - osm_union` (differenced)
4. `osm_only = osm_union - official_union`
5. Combine official_only and osm_only
6. Merge adjacent polygons
7. Label with source

This is simpler and cleaner. Let me use this approach.

But wait, step 6: "merge adjacent polygons" - I need to dissolve/union all the polygons and then explode multi-part to single part, then merge adjacent ones.

Actually, in geopandas, if I just `dissolve()` all features and then `explode()`, I'll get individual polygons. But that doesn't merge adjacent ones unless they're touching.

The proper way: 
- After computing official_only and osm_only, combine them into one GeoDataFrame
- Use `dissolve()` with no group to get a single MultiPolygon
- Then `explode()` to get individual polygons (which are now non-overlapping but may be adjacent)
- To merge adjacent polygons, use `unary_union()` which merges touching polygons
- Then `explode()` again

Wait, unary_union already merges touching polygons. So:

1. official_only = official_union.difference(osm_union)
2. osm_only = osm_union.difference(official_union)
3. Combine into one geometry collection
4. unary_union → this merges touching polygons into MultiPolygons
5. For each resulting polygon/MultiPolygon, determine source

But if I unary_union official_only and osm_only together, everything gets merged into one big thing... 

No wait, `unary_union` takes a list of geometries and merges them topologically. If two polygons are adjacent (share a boundary), they get merged/joined. If they're disjoint, they stay as separate parts of a MultiPolygon.

So `unary_union([official_only, osm_only])` would merge adjacent official-only and osm-only polygons together! That's exactly what we want.

Then:
- If a resulting MultiPolygon contains area from both official and osm → source = "both"
- If only from official → "parks_official"
- If only from osm → "parks_osm"

But how do I know which parts came from which? I need to keep track.

Alternative approach: 
1. Compute official_only
2. Compute osm_only
3. For official_only, set source = "parks_official"
4. For osm_only, set source = "parks_osm"
5. Combine into one GeoDataFrame
6. Use `dissolve()` with some grouping... 

Hmm, this is tricky. Let me think differently.

Actually, I think the best approach is:

1. Compute the overlay union of the two layers to get all the pieces
2. For each piece, determine source
3. Dissolve by source... no, because we want adjacent pieces with DIFFERENT sources to merge too.

OK, let me try a different approach:

1. Compute official_only (official - osm)
2. Compute osm_only (osm - official)
3. Combine them
4. Create a spatial overlay (intersection) of each piece with the original official and osm polygons to track source
5. Merge adjacent polygons

Actually, maybe simpler:

Step 1: official_only = official_union - osm_union
Step 2: osm_only = osm_union - official_union 
Step 3: Tag official_only with source='parks_official', osm_only with source='parks_osm'
Step 4: Combine into one GeoSeries
Step 5: unary_union → this merges touching polygons regardless of their original source
Step 6: Explode into individual geometries (each is a Polygon or MultiPolygon)
Step 7: For each geometry, determine source by checking intersection with official_only and osm_only

For step 7, for each merged geometry G:
- If G.intersects(official_only_area) and G.intersects(osm_only_area) → "both"
- If G.intersects(official_only_area) only → "parks_official"
- If G.intersects(osm_only_area) only → "parks_osm"

But this has issues with precision/rounding. Better to check representative points.

Actually, I think the cleanest approach is:

1. Compute the symmetric difference between the two layers at the collection level
2. Merge adjacent pieces
3. Determine source for each merged piece

Let me just code it up and handle it step by step.

Here's my plan:

```python
import geopandas as gpd
from shapely.ops import unary_union, polygonize
from shapely.geometry import MultiPolygon

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

# Get the union of all geometries in each layer
official_union = unary_union(official.geometry)
osm_union = unary_union(osm.geometry)

# Symmetric difference components
official_only = official_union.difference(osm_union)
osm_only = osm_union.difference(official_union)

# Also compute the intersection (both)
both_areas = official_union.intersection(osm_union)

# Now, we want "every patch claimed by exactly one side" = official_only + osm_only
# Merge adjacent patches

# Create a combined GeoSeries with source tracking
# We need to handle the case where official_only and osm_only are MultiPolygons
```

Hmm, but tracking source through unary_union is tricky since unary_union merges everything without preserving attributes.

Let me try a different approach using polygonization:

1. Take official_only as one set, osm_only as another
2. Label each with source
3. Use unary_union on both to merge adjacent ones
4. Explode result
5. For each exploded feature, intersect it with the original labeled pieces to determine source

Actually, let me think about this more concretely.

After step 3, `unary_union([official_only, osm_only])` gives a merged collection where:
- Disjoint polygons stay separate (as MultiPolygon parts)
- Touching polygons (official-only adjacent to osm-only) get merged into MultiPolygon with merged boundaries

Wait, `unary_union` on touching polygons actually merges them into a single Polygon (dissolving the shared boundary). It doesn't keep them as separate parts of a MultiPolygon.

So after unary_union, if an official-only polygon and an osm-only polygon are adjacent, they'd be merged into a single Polygon. Then I can't tell which parts were official vs osm just by looking at the geometry.

But I CAN compare: for a merged polygon P, does P.intersection(official_only) have area > 0? Does P.intersection(osm_only) have area > 0?

Since official_only and osm_only are disjoint (by construction), the intersection of P with each tells us the source.

But floating point issues... Let me use a tolerance.

OK let me just write the code and test it.

Actually, I realize there might be issues with unary_union merging touching polygons from different sources into a single Polygon. Let me check if that's what we want.

Yes, "a shifted park's gain and loss form a single feature" - so if the gain (OSM-only) and loss (Official-only) are adjacent, they should merge into one Polygon (not MultiPolygon). And the source would be "both" since it contains areas from both.

OK, let me code it.

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

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

official_union = unary_union(official.geometry)
osm_union = unary_union(osm.geometry)

official_only = official_union.difference(osm_union)
osm_only = osm_union.difference(official_union)

# Merge adjacent disagreement patches
merged = unary_union([official_only, osm_only])

# Handle empty geometries
if merged.is_empty:
    merged_geoms = []
elif merged.geom_type == 'Polygon':
    merged_geoms = [merged]
elif merged.geom_type == 'MultiPolygon':
    merged_geoms = list(merged.geoms)

# For each merged geometry, determine source
results = []
for geom in merged_geoms:
    # Check overlap with official_only and osm_only
    has_official = geom.intersection(official_only).area > 1e-6
    has_osm = geom.intersection(osm_only).area > 1e-6
    
    if has_official and has_osm:
        source = 'both'
    elif has_official:
        source = 'parks_official'
    else:
        source = 'parks_osm'
    
    results.append({'geometry': geom, 'source': source})

# Create GeoDataFrame
result_gdf = gpd.GeoDataFrame(results, crs=official.crs)
```

Wait, but `official_only` is the result of `difference()`. If official has many non-overlapping parks, `official_only` is a MultiPolygon with many parts. Similarly for `osm_only`.

When I do `unary_union([official_only, osm_only])`, this takes the two (Multi)Polygons and merges them. Since official_only and osm_only are disjoint (they don't overlap, they can only touch at boundaries), the result is essentially these geometries merged together with shared boundaries dissolved.

Let me make sure this works as I expect. If official_only has parts A, B, C and osm_only has parts D, E, F where A touches D but all others are disjoint, then:
- A and D merge into a single Polygon P1
- B remains as Polygon P2
- C remains as Polygon P3
- E remains as Polygon P4
- F remains as Polygon P5

Then P1 would have source "both", P2 would be "parks_official", etc.

This is exactly what we want!

Let me code this properly.

Actually, I also need to handle the case where the merged result is a MultiPolygon but some of its sub-geometries should have different sources. In a MultiPolygon from unary_union, each sub-geometry is a polygon, and they could have different sources. But since we're using "one MultiPolygon per merged group", I think each merged group is one Polygon (or MultiPolygon if there are disjoint parts that somehow ended up together).

Wait, unary_union returns either a Polygon or MultiPolygon. If the result is a MultiPolygon with parts that are all disjoint, each part is a separate "merged group." So I should handle each part individually.

Let me refine:

```python
merged = unary_union([official_only, osm_only])
if merged.geom_type == 'MultiPolygon':
    for part in merged.geoms:
        # ...
elif merged.geom_type == 'Polygon':
    # handle single polygon
```

OK, I think my approach is correct. Let me write the full script.

For label anchors, I'll use `geom.representative_point()` which returns a point guaranteed to be inside the geometry.

Let me write the full solution.

Let me understand the problem:

1. We have two layers in a GeoPackage: `parks_official` and `parks_osm`.
2. Both have park polygons with a `park_id` field.
3. We need to find areas where the parks disagree - i.e., areas that are in one layer but not the other.
4. "Adjacent disagreement patches merged together" - so if there are multiple disconnected patches of disagreement that are adjacent (touching), they should be merged into one group.
5. Each group gets a `source` attribute: `parks_official` (area only in official), `parks_osm` (area only in OSM), or `both` (overlap area... wait).

Actually, re-reading: "I need every patch claimed by exactly one side" - so we want areas that are in official only (official - OSM) and areas that are in OSM only (OSM - official). These are the "disagreement" patches.

Then "with adjacent disagreement patches merged together" - if a park shifted, the area lost (only in official) and the area gained (only in OSM) might be adjacent, so they should be merged.

Wait, let me re-read: "so that a shifted park's gain and loss form a single feature rather than two" - so if a park shifted, the old area (only in official) and new area (only in OSM) are adjacent, and we want to merge them into one feature with source... hmm, but they have different sources.

Actually, I think they mean: find all areas that are in exactly one of the two sources. Take both the "official-only" patches and "OSM-only" patches. Then, any adjacent patches (regardless of source) get merged into a single MultiPolygon. So a shifted park might produce two adjacent disagreement patches (one official-only where the park was, one OSM-only where the park moved to), and they get merged into one feature.

But what would the `source` attribute be for a merged group containing both official-only and OSM-only patches? Let me re-read...

"For each group I also need a label-anchor point that falls strictly inside the group's geometry."

I think the approach is:
1. Compute `official_only = official - osm` (differenced)
2. Compute `osm_only = osm - official` (differenced)
3. Combine both into one set of disagreement patches
4. Merge/disjoint any patches that are adjacent (touching)
5. Assign source based on which layer(s) contributed

Actually wait, if we merge an official-only patch with an OSM-only patch, the source would be both? No... the prompt says "source attribute whose value is exactly one of parks_official, parks_osm, or both". Hmm.

Let me re-read one more time: "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 the "disagreement patches" are areas claimed by exactly one side. But then adjacent ones (which could be official-only adjacent to OSM-only) get merged. The merged group would have source = both if it contains both types.

Actually, I think a cleaner interpretation: the "disagreement" is the symmetric difference (official XOR osm). Then we group these into contiguous clusters. If a cluster was formed from both official-only and OSM-only areas, its source is `both`; if purely official-only, source is `parks_official`; if purely OSM-only, source is `parks_osm`.

Wait, but that doesn't quite make sense either. Let me think again about what "source" means in context.

Actually, I think the simplest correct interpretation:

1. Compute the symmetric difference between official and OSM parks.
2. But before computing, we need to figure out which official parks correspond to which OSM parks (by park_id, since they share the same UUIDs).
3. For each pair of matching parks (same park_id), compute their symmetric difference to find disagreement areas.
4. For parks only in one source, the entire park is a disagreement area.
5. Then, merge any adjacent disagreement patches regardless of source.

Hmm, but that's complex. Let me think about it differently.

Actually, I think the approach is:

1. Union all official park polygons
2. Union all OSM park polygons
3. Compute:
   - `only_official = union_official - union_osm` (areas only in official)
   - `only_osm = union_osm - union_official` (areas only in OSM)
4. These are the "disagreement patches"
5. Now merge adjacent patches: take all disagreement patches and do a dissolve/unary_union on adjacency
6. For each merged group, determine if it contains only official-only parts, only OSM-only parts, or both

Actually, let me think more carefully. The problem says "a shifted park's gain and loss form a single feature." So if park A shifted north, the southern part that was in official but not OSM, and the northern part that was in OSM but not official, are adjacent and should be one feature.

So:
- Compute `official_only = official_union - osm_union`
- Compute `osm_only = osm_union - official_union`
- Label official_only polygons with source='parks_official'
- Label osm_only polygons with source='parks_osm'
- Now, merge adjacent polygons regardless of source
- If a merged group contains both types, label as 'both'

For "merge adjacent": we need to find which polygons touch/overlap and group them. Use `unary_union` (which merges all adjacent things) won't work because we'd lose the source information. We need to union within groups that are connected.

One approach: take all disagreement polygons, assign each a unique ID, then use `dissolve` with `by='connected_component'` where we compute the connected components of the adjacency graph.

Actually, `unary_union` merges everything that intersects. But that changes the geometry. Let me think of a better approach.

We can use `shapely.ops.unary_union` to union the geometries, but we need to first assign source info. 

Better approach:
1. Create GeoDataFrame with all disagreement areas, each with a source attribute.
2. Use spatial join / intersection to find which polygons are adjacent (touching or overlapping).
3. Build a graph of adjacency.
4. Find connected components.
5. For each component, union all geometries into a single MultiPolygon, and set source based on the sources of the components.

Let me implement this.

Actually, I think the simplest approach:

1. `official_union = official.unary_union`
2. `osm_union = osm.unary_union`
3. `only_official = official_union.difference(osm_union)` -> this gives a MultiPolygon/GeometryCollection
4. `only_osm = osm_union.difference(official_union)` -> same
5. Convert each to individual Polygons
6. Combine into one GeoDataFrame with source labels
7. Compute adjacency and merge

For step 7, I can:
- Compute a buffer of 0 on each polygon (to handle exact adjacency)
- Then use `unary_union` which merges intersecting/touching polygons
- Then explode back to individual polygons/MultiPolygons
- For each resulting merged shape, determine its source by checking which original source polygons contributed

Actually wait, the buffer(0) + unary_union approach loses source info. Let me use a different approach.

```python
import networkx as nx

# All disagreement polygons as individual rows
gdf = GeoDataFrame with columns: id, source, geometry

# Build adjacency
adjacency = []
for i in range(len(gdf)):
    for j in range(i+1, len(gdf)):
        if gdf.iloc[i].geometry.touches(gdf.iloc[j].geometry) or gdf.iloc[i].geometry.intersects(gdf.iloc[j].geometry):
            adjacency.append((i, j))

G = nx.Graph()
G.add_edges_from(adjacency)
components = list(nx.connected_components(G))

# For each component, merge geometries and determine source
for comp in components:
    merged = unary_union([gdf.iloc[i].geometry for i in comp])
    sources = set(gdf.iloc[i].source for i in comp)
    if len(sources) == 1:
        source = sources.pop()
    else:
        source = 'both'
```

But this could be slow for many polygons. Let me check how many disagreement patches there might be.

Actually, with 1380 and 1372 features, the symmetric difference probably isn't huge. Let me just write the script.

Also, for the label anchor point - "a label-anchor point that falls strictly inside the group's geometry". I can use `representative_point()` from shapely, which guarantees a point within the geometry.

Let me write the full solution.

One more thing to consider: when computing `official_union - osm_union`, the result might be a GeometryCollection with mixed types. I need to extract only polygons/multipolygons.

Let me also consider: maybe I should do this per-park-id rather than globally. If two parks are adjacent in both official and OSM, their difference might accidentally create a single patch spanning both parks. But the problem says "a shifted park's gain and loss" - so it's per-park context.

Hmm, but if I do it per park_id:
- For matching park_ids, compute the difference between official and OSM polygons for that specific park
- For park_ids only in official, the whole polygon is official-only
- For park_ids only in OSM, the whole polygon is OSM-only

This seems cleaner and more correct. Let me do this.

Actually, there are 1380 official and 1372 OSM features with 1360 matching park_ids. So 20 official-only and 12 OSM-only.

For matching park_ids (1360 pairs), I need to compute the symmetric difference area for each pair.

For the 20 official-only: the entire polygon is a disagreement patch (source: parks_official)
For the 12 OSM-only: the entire polygon is a disagreement patch (source: parks_osm)
For each matching pair: 
  - official_diff = official_geom - osm_geom (source: parks_official)
  - osm_diff = osm_geom - official_geom (source: parks_osm)
  - These could be empty if the geometries match exactly

Then, merge adjacent disagreement patches.

This approach seems right. Let me code it up.

Actually, I just realized - for the "both" source: when would something have source "both"? Let me re-read...

"with a source attribute whose value is exactly one of parks_official, parks_osm, or both"

So source can be:
- `parks_official` - the merged group contains only official-only parts
- `parks_osm` - the merged group contains only OSM-only parts
- `both` - the merged group contains both official-only AND OSM-only parts (i.e., a shifted park where we merged the loss and gain)

This makes sense. If a park shifted, the old location (official-only) and new location (OSM-only) are adjacent, they get merged, and the merged feature has source 'both'.

Let me code this up now.

Let me think carefully about the requirements again.

"We 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 the task says "patches claimed by exactly one side" - but then says to merge adjacent disagreement patches, with a source attribute of `parks_official`, `parks_osm`, or `both`.

Wait, 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`"

So if we merge an official-only patch adjacent to an OSM-only patch, the merged group gets 'both'. That makes sense as a combined disagreement zone.

But wait - actually re-reading more carefully: "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"

When a park is shifted, the area that's only in official (loss) and the area that's only in OSM (gain) would be adjacent. They merge into one feature. That feature should get 'both' as its source because it represents both a gain and a loss.

But what about patches that are ONLY official-only or ONLY OSM-only and not adjacent to the other? Those would have source 'parks_official' or 'parks_osm'.

My current implementation does this correctly. Let me check if there could be any edge cases.

Wait, but there's another important consideration. If a polygon in official and OSM overlap exactly, the difference would be empty, so no patch. That's fine.

But what about a park that exists ONLY in official and is adjacent to a park that exists ONLY in OSM? They'd get merged and become 'both'. Hmm, is that desired? The task says "a shifted park's gain and loss" - suggesting it's the SAME park that shifted. So if park X exists in both but its geometry differs, the difference areas might be adjacent. And for parks only in one source, they'd be standalone patches.

My approach handles this fine - the merging is purely geometric adjacency regardless of whether they came from the same park_id or not.

But wait - let me think about the "both" cases more carefully. If an official-only area is adjacent to an OSM-only area (from the same or different parks), they'd merge into a 'both' group. That seems correct per the spec.

Actually, hmm, I need to be careful about what "every patch claimed by exactly one side" means. This says every PATCH (individual disagreement piece) is claimed by exactly one side. Then "adjacent disagreement patches merged together" - so we merge adjacent patches (each claimed by one side) into groups. The group's source is 'both' if it contains patches from both sides.

Let me also think about the label-anchor requirement: "a label-anchor point that falls strictly inside the group's geometry". I used `representative_point()` which is guaranteed to be within the geometry, not just the centroid. Good.

Let me run it and check the output.

Let me verify the strict containment condition more carefully - `representative_point()` should return a point strictly inside, but let me double check.

Actually, `representative_point()` in Shapely returns a point guaranteed to be within the object, not on its boundary. So it should be strictly inside. But the check I wrote might give false negatives for MultiPolygons where the point is in one polygon but `contains` checks the whole MultiPolygon. Let me also double check that the `both` source is correct - it means patches that come from both sides touching each other.

The solution looks correct. Let me verify the files exist.