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 claimed by exactly one side (disagreement areas)
3. Adjacent disagreement patches should be merged together
4. Each merged group gets a `source` attribute: `parks_official`, `parks_osm`, or `both`
5. For each group, I need a label-anchor point strictly inside the geometry
6. Output: `parks_disagreement.geojson` and `park_label_anchors.geojson`

Let me first explore 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. I need to find areas that are in one but not the other (symmetric difference), but merge adjacent disagreement patches
3. Each merged group should be a MultiPolygon with a `source` attribute: `parks_official` (in official only), `parks_osm` (in OSM only), or `both` (in both, so not a disagreement - wait, re-reading...)

Wait, re-reading the task: "I need every patch claimed by exactly one side" - so only areas that are in one but not the other.

"with adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two" - A park that shifts will have some area gained (new area not in official) and some area lost (old area not in OSM). If these are adjacent (touching), they should be merged into one feature.

Let me think about this more carefully.

For each park area, we have:
- Areas in `parks_official` only (source = `parks_official`)
- Areas in `parks_osm` only (source = `parks_osm`)

But the task says "adjacent disagreement patches merged together" - so if a `parks_official`-only patch is adjacent to a `parks_osm`-only patch, they get merged into a single feature with what source?

Wait, I think the instruction is saying: "every patch claimed by exactly one side" - so only the symmetric difference. And then if two such patches are adjacent (one being official-only and the other being OSM-only), they get merged together.

But what source would the merged feature have? Let me re-read...

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

Hmm, if we merge an official-only area and an OSM-only area together, the source... could be something other than those three? No, it says exactly one of those three strings.

Wait, maybe the "both" case is actually part of the disagreement. Let me think about this differently.

Actually, I think the "disagreement" is NOT the symmetric difference. Let me re-read:

"I need every patch claimed by exactly one side" - this is the symmetric difference (official-only ∪ osm-only)

"with adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two" - So an official-only patch that touches an osm-only patch should be merged into one feature.

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

So the merged group could have source = `both` if it contains both official-only and osm-only areas? That would make sense for the "shifted park" scenario.

Let me think about this algorithm:

1. Compute official-only areas: `parks_official` - `parks_osm` (source = `parks_official`)
2. Compute OSM-only areas: `parks_osm` - `parks_official` (source = `parks_osm`)
3. Now, merge adjacent patches that have different sources.

For step 3, we could:
- Take all the patches (from both sides), dissolve them with adjacency.
- For each connected group, check what sources are present.
- If only one source present: that's the source.
- If both sources present: source = `both`.

This way, a shifted park where the gain (OSM-only) touches the loss (official-only) becomes a single feature with source = `both`.

Let me implement this.

Actually, wait. Let me reconsider. "every patch claimed by exactly one side" - this is the symmetric difference. Then "adjacent disagreement patches merged together" - so if an official-only polygon touches an OSM-only polygon, they become one group. 

For each group, if it contains only official-only areas → source = `parks_official`
If it contains only OSM-only areas → source = `parks_osm`
If it contains both → source = `both`

Then for each group, I need a label anchor point that falls strictly inside the group's geometry.

Let me write the script.

Algorithm:
1. Compute the symmetric difference
2. Dissolve/merge adjacent patches
3. Assign source based on which original sources are present
4. Compute label anchor points (representative points that fall strictly inside)

Let me think about how to do step 2. I'll use `unary_union` with `dissolve` on a combined GeoDataFrame.

Actually, let me think about this more carefully with shapely/geopandas.

Step 1: Compute official-only and OSM-only areas.
```python
off_only = gpd.overlay(off, osm, how='difference')
osm_only = gpd.overlay(osm, off, how='difference')
```

But overlay can be slow for 1380 polygons. Let me think of a better approach.

Actually, I think using `unary_union` to merge all official polygons, merge all OSM polygons, then take symmetric difference might be cleaner at the layer level.

```python
off_union = off.geometry.unary_union
osm_union = osm.geometry.unary_union
off_only = off_union.difference(osm_union)
osm_only = osm_union.difference(off_union)
```

Then I need to explode these into individual polygons, assign source tags, merge adjacent ones, etc.

Let me think about the merging step more carefully.

After getting `off_only` and `osm_only` as MultiPolygon/GeometryCollection, I can explode them into individual polygons, assign source tags to each polygon based on which set they came from.

Then dissolve polygons that touch each other regardless of source using`dissolveby` with adjacency — actually geopandas dissolve aggregates by attribute and dissolves touching polygons within each group. But since touching polygons with different sources should ALSO be merged into one group with source 'both', we need to handle adjacency across source labels.

One approach:
1. Gather all polygons with their source tags
2. Use `unary_union` on ALL polygons withST_ClusterDBSCAN or similar to find groups of adjacent polygons regardless of source tag
Actually simpler approach using GeoPandas spatial operations:

Let me use `shapely.ops.unary_union` but that would merge everything including across sources. Instead:

1. Create a GeoDataFrame with all polygons, each with a source tag
2. Use DBSCAN clustering with eps=0 (touching) to cluster adjacent polygons
3. Dissolve by cluster, then assign source based on what's in the cluster

Let me code this up.

For DBSCAN clustering of touching polygons:
```python
from sklearn.cluster import DBSCAN
import numpy as np

But I need to compute a distance matrix where touching = distance 0.

Actually, a simpler approach: use `dissolve` with a custom grouping based on spatial adjacency.

Let me think of another approach using `shapely` ops:

```python
from shapely.ops import unary_union, polygonize
```

Hmm, actually the simplest approach might be:

1. Take all polygons (both off_only and osm_only), each with source tag.
2. Use `dissolve` with a connectivity graph.

Let me try a different approach. Since the number of patches might not be that large:

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

# Get all patches with source
patches = []
for geom in off_only.geoms:  # if MultiPolygon
    patches.append((geom, 'parks_official'))
for geom in osm_only.geoms:
    patches.append((geom, 'parks_osm'))

# For each patch, find which group it belongs to by adjacency
# Use a graph approach
```

Let me think about this differently. Actually, let me use a union-dissolve approach:

All patches together: `all_union = unary_union(all_patches)`. This dissolves everything into multipolygons where each multipolygon is a connected group. Then I can iterate over each connected component in the result, and for each, check which source(s) the original patches had.

Wait, but `unary_union` of all patches will just merge them all including across sources into disconnected multipolygons. That's exactly what I want! Each sub-polygon in the result represents a connected component of adjacent patches. But I need to know which sources contributed to each component.

Hmm, but when I unary_union them, I lose track of which source each part came from.

Alternative approach: use spatial indexing to find which original patches overlap with which component of the unary_union.

```python
all_merged = unary_union([geom for geom, src in patches])
# all_merged is a MultiPolygon of connected groups
for component in all_merged.geoms:
    # Find which patches intersect this component
    sources = set()
    for geom, src in patches:
        if component.intersects(geom):
            sources.add(src)
```

But this might be a bit brute-force. Let me think if there's a cleaner way.

Actually, I think using `geopandas` and `dissolve` with DBSCAN clustering is the cleanest.

```python
from sklearn.cluster import DBSCAN
import numpy as np

# Build GeoDataFrame
gdf = gpd.GeoDataFrame({'source': sources, 'geometry': geometries})
gdf.crs = off.crs

# Get centroids or representative points for clustering
# Use spatial index and the geometry's boundary distance
coords = np.array([(geom.centroid.x, geom.centroid.y) for geom in gdf.geometry])

# But centroid-based clustering won't correctly identify touching polygons...
```

DBSCAN with distance = 0 doesn't work for touching polygons with the centroid approach. I really need a proper adjacency matrix.

Let me use a different approach. I'll build an adjacency graph manually using spatial index:

```python
from shapely import STRtree

tree = STRtree(geometries)
adjacency = {i: set() for i in range(len(geometries))}
for i, g in enumerate(geometries):
    candidates = tree.query(g)
    for j in candidates:
        if i != j and g.touches(geometries[j]):
            adjacency[i].add(j)
```

Then find connected components in this graph, dissolve by component.

This is doable. Let me write the full script.

Actually, an even simpler approach: use GeoPandas' `dissolve` with `by` parameter on a spatial attribute. No, that doesn't work directly.

Or use `geopandas.sjoin` or overlay to find intersections...

OK let me just go with the graph approach. Let me write the full script.

For the label anchors: I need a point "strictly inside" the group's geometry. I can use `representative_point()` which guarantees it's inside, or `point_on_surface()`.

Let me write the script now.

Actually wait - I want to use `STRtree` from shapely. Let me check which version is available.

Also, thinking about the symmetrical difference computation more carefully.

Instead of overlay which can be slow and complex, let me use `unary_union` and `difference`:

```python
off_union = off.geometry.unary_union
osm_union = osm.geometry.unary_union

off_only_geom = off_union.difference(osm_union)
osm_only_geom = osm_union.difference(off_union)
```

Then explode these into individual polygon patches.

Let me write the code:

```python
import geopandas as gpd
from shapely.geometry import shape, Polygon, MultiPolygon
from shapely import STRtree
import numpy as np

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

# Get union of each
off_union = off.geometry.unary_union
osm_union = osm.geometry.unary_union

# Symmetric difference components
off_only = off_union.difference(osm_union)
osm_only = osm_union.difference(off_union)

# Collect individual polygons with source
patches = []  # list of (geometry, source)
if off_only.geom_type == 'MultiPolygon':
    for geom in off_only.geoms:
        patches.append((geom, 'parks_official'))
elif off_only.geom_type == 'Polygon':
    patches.append((off_only, 'parks_official'))

if osm_only.geom_type == 'MultiPolygon':
    for geom in osm_only.geoms:
        patches.append((geom, 'parks_osm'))
elif osm_only.geom_type == 'Polygon':
    patches.append((osm_only, 'parks_osm'))

print(f"Total patches: {len(patches)}")

# Build adjacency graph
geometries = [p[0] for p in patches]
sources = [p[1] for p in patches]
n = len(geometries)

tree = STRtree(geometries)
adj = {i: set() for i in range(n)}

for i, g in enumerate(geometries):
    # query returns indices of geometries that intersect
    candidates = tree.query(g)
    for j in candidates:
        if i != j and g.touches(geometries[j]):
            adj[i].add(j)
            adj[j].add(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 merged groups: {len(components)}")

# Create result features
result_geoms = []
result_sources = []

for comp in components:
    comp_geoms = [geometries[i] for i in comp]
    comp_sources = set(sources[i] for i in comp)
    
    # Dissolve geometries
    dissolved = unary_union(comp_geoms)
    
    # Determine source
    if len(comp_sources) == 1:
        src = list(comp_sources)[0]
    else:
        src = 'both'
    
    result_geoms.append(dissolved)
    result_sources.append(src)

# Create GeoDataFrame for disagreement polygons
gdf_disagreement = gpd.GeoDataFrame({
    'source': result_sources,
    'geometry': result_geoms
}, crs=off.crs)

# Compute label anchors (representative point inside geometry)
anchors = []
for geom in gdf_disagreement.geometry:
    anchors.append(geom.representative_point())

gdf_anchors = gpd.GeoDataFrame({
    'source': gdf_disagreement['source'].values,
    'geometry': anchors
}, crs=off.crs)

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

print("Done!")
```

Wait, I need to also handle the case where the symmetric difference produces a GeometryCollection. Let me also handle empty results.

Also, I need to make sure `representative_point()` returns a point strictly inside. Let me check: `representative_point()` from shapely returns a point that is guaranteed to be within the geometry. Yes, it's documented as "Returns a point that is guaranteed to be within the geometry."

Let me also think about edge cases:
- If there are 0 patches, we should handle that gracefully.
- The STRtree query returns indices... let me check the shapely version.

Let me just write the script and run it.

One concern: `STRtree` in newer shapely (2.x) returns indices, in older shapely (1.x) returns objects. Let me check what's installed.

Let me understand the problem:

1. We have two layers: `parks_official` and `parks_osm` in a GPKG
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 - so if a park shifted, the gain and loss areas merge into one feature
4. The result should have a `source` attribute: `parks_official` (only official claims it), `parks_osm` (only OSM claims it), or `both` (both claim it)
5. Need label-anchor points strictly inside each group's geometry
6. Output two GeoJSON FeatureCollections

Let me think about the approach:

1. Compute symmetric difference between the two layers (areas where they disagree)
2. Also compute intersection (areas where they agree - "both")
3. But the instruction says "every patch claimed by exactly one side, with adjacent disagreement patches merged together"

Wait, 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 shifted - the official has area A, OSM has area B overlapping partially. The symmetric difference has:
- A \ B (official only)
- B \ A (osm only)

But if these are adjacent (touching), they should be merged into a single feature with source? Wait, that doesn't make sense - they have different sources.

Actually, I think the goal is: for a shifted park, there's area A\B (loss) and area B\A (gain). Both are disagreement patches. If they touch/are adjacent, they should be merged into one feature. But what would the source be?

Hmm, let me re-read: "every patch claimed by exactly one side... with adjacent disagreement patches merged together... One MultiPolygon per merged group with a source attribute whose value is exactly one of parks_official, parks_osm, or both"

So the groups have source = parks_official OR parks_osm OR both.

Wait, I think maybe the approach is:
- Compute the union of official and osm
- For each piece of the union, determine which source(s) claim it
- If both claim it → source = both
- If only official claims it → source = parks_official
- If only osm claims it → source = parks_osm
- Then merge adjacent pieces that have the same source

Actually wait - "adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature." This means if a park in the official is shifted relative to OSM, the gain (only OSM) and loss (only official) patches are adjacent and should be merged.

But they have different sources! Unless the merged group has some other source...

Hmm, or maybe it's: merge all adjacent disagreement patches regardless of source? Then what's the source attribute?

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

I think maybe the approach is:
1. Compute the full geometry union of both layers
2. For each polygon piece in the union, tag it with which source(s) contain it
3. Dissolve/merge adjacent pieces that share the same source tag
4. Output the result

Let me think again about what "adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two" means.

If park P in official is in location L1, and the same park in OSM is in location L2 (slightly shifted), then:
- L1 \ L2 = areas only in official (loss) 
- L2 \ L1 = areas only in OSM (gain)

These two areas are adjacent (they touch along the boundary). The user wants them merged into one feature. But what source would that merged feature have? It can't be "both" because it's only claimed by one side each.

Hmm, unless the interpretation is: we're looking at the two layers as a whole, doing an overlay (union), and for each piece of the pie, it's claimed by some subset of layers. Then we merge adjacent pieces with the same source. 

So the approach:
1. Overlay = union of official and OSM, producing fragments
2. For each fragment, check if it's contained by official only, OSM only, or both
3. Dissolve adjacent fragments with same source classification
4. Output

Let's think about what "adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature" means in this framework.

A shifted park: Official has polygon A, OSM has polygon B (overlapping). The union gives us:
- A\B → only official
- A∩B → both
- B\A → only OSM

If we merge adjacent same-source fragments, A\B and B\A are NOT same source, so they wouldn't merge. But the user wants them to merge...

Wait, maybe I'm misreading. 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."

"disagreement patches" = patches claimed by exactly one side. So A\B and B\A are both disagreement patches. If they're adjacent (touching), they should be merged. The merged feature would contain areas claimed by exactly one side, but parts have different sources. So what source attribute?

Hmm, maybe the merged feature would have source = some other value? No, the instruction says exactly one of the three strings.

Or maybe I need to think of it differently. Perhaps "adjacent disagreement patches" means patches of disagreement that are adjacent to each OTHER, regardless of which side they belong to. And when merged, the source attribute just describes what the merged area is... 

Actually, let me reconsider. Maybe the interpretation is:
- For each park (by park_id), compute its official footprint and OSM footprint
- Take the symmetric difference (disagreement)
- Adjacent pieces of the symmetric difference (even from different sources) merge into one feature
- But... the source would still be ambiguous

Actually, I think I might be overcomplicating it. Let me re-read:

"I need every patch claimed by exactly one side"

OK so both disagreement types.

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

So the two types of disagreement (official-only and osm-only) are adjacent and should become ONE merged feature.

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

So each merged group gets ONE of these three source values.

For a merged group containing both official-only AND osm-only areas, the source can't be just one of the layers. Maybe it's "both"? No, "both" means both claim it, but for disagreement areas, only one side claims any given point.

Hmm, I think I might be misunderstanding the merging. Let me try a different interpretation:

Maybe we should:
1. Compute three sets: official-only, osm-only, both (intersection)
2. Merge adjacent features within each set
3. Output with source = parks_official, parks_osm, both respectively

But this doesn't match "adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature."

Let me try yet another interpretation. Maybe the user wants:

For each park:
- If official and OSM geometries differ:
  - Take the symmetric difference (disagreement area)
  - This may have multiple parts. Merge them all into one MultiPolygon
  - Source = ??? 

Hmm, the source for a symmetric difference that contains both official-only and OSM-only areas...

Wait, maybe the merging is across BOTH layers, not per park. And the "gain and loss" refers to the fact that when you look at the symmetric difference of the full union of both layers, adjacent disagreement patches (official-only touching osm-only) should be merged.

Let me look at it from a purely geometric perspective:
- Compute the full union of both layers (dissolve each layer, then union)
- For each fragment, label it with its source
- Dissolve adjacent fragments with... hmm

Actually, I think I need to approach this differently. Let me try another interpretation:

1. Take the union overlay of official and OSM (without dissolving within each layer first)
2. Each output piece gets attribute: which layer(s) contain it
3. Merge adjacent pieces with the SAME source classification
4. Output the three classes: official-only, osm-only, both

But this means a shifted park's A\B and B\A would NOT be merged (different sources). So this can't be right.

Let me try yet another interpretation:

For the symmetric difference (disagreement):
- All official-only patches from all parks
- All OSM-only patches from all parks
- If any official-only patch touches any OSM-only patch, merge them
- The merged group's source... hmm

OK, I think I'm going to try a different approach that might work. Let me interpret the task as:

1. Compute the union overlay (all-pairs intersection) between the two layers
2. Tag each fragment: which layer(s) it belongs to
3. Dissolve/merge adjacent fragments that have the same tag
4. This gives us three sets: official-only, osm-only, both
5. For the disagreement (official-only + osm-only), if an official-only piece and osm-only piece are adjacent, merge them
6. The source for such a merged piece... let me check the instruction again

"I need every patch claimed by exactly one side, with adjacent disagreement patches merged together"

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

So each merged group gets ONE of these three values. If official-only and osm-only areas are merged together, what source do we assign?

Actually, I think the answer might be simpler. Let me re-read:

"I need every patch claimed by exactly one side" - this means we're only interested in disagreement, not the "both" areas.

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

So we take all disagreement patches (A\B and B\A across all features), merge adjacent ones, and the result is one MultiPolygon per merged group. But what source?

Hmm wait, maybe the merged group has both pieces from official-only and osm-only, and the source attribute simply doesn't apply cleanly. But the instruction says source must be one of those three.

Let me go with a practical approach:

1. Compute the union overlay between official and OSM layers (using overlay with union mode)
2. For each fragment, determine which layers contain it
3. Classify each fragment: "official_only", "osm_only", "both"
4. Dissolve/merge adjacent fragments with the same classification
5. Filter to only disagreement (official_only + osm_only)
6. Within the disagreement set, dissolve/merge ALL adjacent features regardless of classification (this is the "adjacent disagreement patches merged together" part)
7. For these merged groups, I need to determine the source attribute

For step 7, maybe the source attribute for a merged group is just a label of convenience. Let me think...

Actually maybe I'm way overthinking this. Let me try the simplest plausible interpretation:

1. Compute the overlay union of the two layers
2. Each fragment has one of three types: official-only, osm-only, both
3. Merge adjacent fragments regardless of type
4. For each group, determine if it has only official-only → source=parks_official; only osm-only → source=parks_osm; contains both types or the "both" type → source=both

Wait, that could work! If a shifted park's A\B (official-only) and B\A (osm-only) are adjacent, they merge into one group. That group contains both official-only and osm-only areas, so its source would be "both". That makes sense!

Let me test this logic:
- Official-only area: only claimed by official, not by OSM
- OSM-only area: only claimed by OSM, not by official
- When merged: contains areas claimed by official (the official-only part) AND areas claimed by OSM (the osm-only part), so the group as a whole is claimed by both sources → "both"

Yes, I think this is the interpretation! The merged group is "claimed by both" in the sense that it contains territory that each side claims, even though no single point is claimed by both.

So the algorithm:
1. Compute union overlay of official and OSM
2. Classify each fragment
3. Dissolve/merge adjacent fragments with the SAME classification into groups
4. NOW, dissolve/merge adjacent groups regardless of classification (only for official-only and osm-only groups? or all three?)
   - Actually, "adjacent disagreement patches merged together" - so official-only and osm-only groups merge together. What about "both" groups? I think "both" groups are not "disagreement patches", so they stay separate.
5. After merging official-only and osm-only groups that touch, these merged groups get source = "both"
6. Remaining official-only groups → source = "parks_official"
7. Remaining osm-only groups → source = "parks_osm"
8. "both" groups that were already "both" → source = "both"

Actually, this is getting complex. Let me simplify:

1. Compute union overlay → fragments with classifications
2. First pass: dissolve adjacent fragments with the same classification → groups
3. Second pass: merge adjacent groups where one is official-only and the other is osm-only (disagreement patches)
   - These merged groups get source = "both"
4. Output all groups

Let me implement this.

Actually, let me think about the "both" areas separately. The "both" areas (intersection) are areas of agreement - both layers claim them. These are NOT disagreement patches. The user asks for "every patch claimed by exactly one side" with disagreement patches merged. So "both" areas are already separate and don't get merged with anything.

Wait, the user says "every patch claimed by exactly one side" - so disagreement only. But then also says source can be "both"... Let me re-read more carefully.

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

I think "every patch claimed by exactly one side" refers to the disagreement patches, and the "adjacent disagreement patches merged together" refers to merging official-only and osm-only patches that touch. The resulting merged groups could have source = "both" (meaning the group contains claims from both sides, even though no single point is claimed by both).

But what about the actual "both" areas (where both layers agree)? I think those should also be in the output with source = "both".

So the full set is:
1. Areas where both agree → source = "both" (not merged with anything - they're not disagreement)
2. Areas where only official claims it → source = "parks_official" (but merged with adjacent osm-only areas)
3. Areas where only OSM claims it → source = "parks_osm" (but merged with adjacent official-only areas)

And after merging official-only and osm-only groups that touch, they become source = "both".

Let me implement this step by step in Python.

Steps:
1. Read both layers
2. Compute the union overlay
3. For each fragment, determine source
4. Dissolve same-source fragments
5. Merge adjacent official-only and osm-only groups

Let me write the script.

Actually, I realize the union overlay approach might be complex with geopandas. Let me think about a simpler approach using shapely operations.

Approach with shapely/geopandas:
1. Dissolve each layer into a single multipolygon
2. Compute:
   - both = intersection(dissolved_official, dissolved_osm)
   - official_only = difference(dissolved_official, dissolved_osm) 
   - osm_only = difference(dissolved_osm, dissolved_official)
3. Now, each of these is a single MultiPolygon
4. But we need to explode and merge adjacent disagreement patches

Actually, dissolving each layer first loses per-park detail. Let me not dissolve individual layers first.

Alternative approach using overlay:

```python
# Union overlay
union = gpd.overlay(off, osm, how='union')
# Now we have fragments with overlapping attributes from both layers
# Classify each fragment
```

But `gpd.overlay` with `how='union'` requires both layers to have the same columns. Let me just use the geometric operations.

Let me try a different approach using shapely directly:

1. For each park in official, check if it exists in OSM by park_id
2. For matching parks, compute differences and intersection
3. For parks only in one layer, they're entirely disagreement
4. Merge adjacent disagreement patches

Actually, the simplest approach:

1. Compute the full overlay (union) using geopandas
2. For each fragment, determine which layer(s) contain it
3. Dissolve by classification
4. Merge adjacent official-only and osm-only fragments

Let me write this step by step.

Actually, let me use `gpd.overlay(off, osm, how='union')`. Both layers have the same columns (park_id, park_class), so the overlay will create fragments with attributes from both sides. But the column names will conflict - geopandas will suffix them.

Hmm, actually let me think about this differently. The union overlay in geopandas creates fragments where each fragment is covered by a unique combination of original features. The output has _1 and _2 suffixes for columns from each input.

After overlay with how='union':
- Some fragments have park_id_1 (from official) and park_id_2 (from OSM) both non-null → "both"
- Some have only park_id_1 non-null → "official_only" 
- Some have only park_id_2 non-null → "osm_only"

Then I can dissolve/merge adjacent fragments with the same classification, and then merge adjacent official_only and osm_only fragments.

Let me write the script.

Wait, geopandas overlay can be slow for large datasets. Let me check: 1380 parks official, 1372 OSM - that's manageable.

Let me write the script.

Actually, I wonder if there's a simpler approach using just shapely operations on dissolved layers.

```python
off_dissolved = off.geometry.unary_union
osm_dissolved = osm.geometry.unary_union

both = off_dissolved.intersection(osm_dissolved)
official_only = off_dissolved.difference(osm_dissolved)
osm_only = osm_dissolved.difference(off_dissolved)
```

But then `both`, `official_only`, `osm_only` are just single (maybe multi) geometries. The "adjacent disagreement patches merged" is automatically handled because within each, all parts are merged by `unary_union`.

Then:
- `official_only` → source = "parks_official"
- `osm_only` → source = "parks_osm"
- `both` → source = "both"

But wait, the user wants "adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature". If I dissolve official_only and osm_only separately, they don't get merged together.

With the dissolved approach:
- official_only = all official-only areas (losses for shifted parks)
- osm_only = all osm-only areas (gains for shifted parks)

These are separate features. Adjacent ones (gain and loss for the same shifted park) are NOT merged.

So I need to merge official_only and osm_only where they touch.

Algorithm:
1. Compute dissolved geometries for both layers
2. Compute both = intersection, official_only = diff(off, osm), osm_only = diff(osm, off)
3. Now, I need to identify which parts of official_only and osm_only are adjacent and merge them
4. For each pair of touching parts from official_only and osm_only, create a merged geometry with source = "both"
5. Leftover official_only parts → source = "parks_official"
6. Leftover osm_only parts → source = "parks_osm"

But this adjacency detection between parts of MultiPolygons is complex. Let me use a polygonization approach instead.

Actually, here's a cleaner approach:

1. Compute the full set of boundary lines of both layers
2. Union the boundaries
3. Polygonize to get all the atomic fragments
4. For each fragment, test against official and OSM layers
5. Classify
6. Merge adjacent same-class fragments
7. Then merge adjacent official-only and osm-only fragments

This is essentially what overlay does, but gives me full control.

Let me use `shapely.ops.unary_union` and `shapely.ops.polygonize` to do this.

Actually, let me just use geopandas overlay. It's the simplest correct approach.

```python
import geopandas as gpd
from shapely.geometry import MultiPolygon

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

# Add source column to identify each layer
off['src'] = 'official'
osm['src'] = 'osm'

# Union overlay - but we need to handle column conflicts
# Let's just use geometry for overlay
union = gpd.overlay(off[['geometry']], osm[['geometry']], how='union')

# Now classify each fragment
# If it came from both inputs, it'll have both index values non-null
# Actually, geopandas overlay with how='union' creates rows where
# the fragment appears in the first, second, or both dataframes

# Hmm, the overlay result columns will be index from left and right.
# Let me check what columns we get.

# Actually, I need to add identification columns first
off['in_off'] = True
osm['in_osm'] = True

# Keep only relevant columns
union = gpd.overlay(off[['in_off', 'geometry']], osm[['in_osm', 'geometry']], how='union')

# Classify
def classify(row):
    if row['in_off'] and row['in_osm']:
        return 'both'
    elif row['in_off']:
        return 'parks_official'
    else:
        return 'parks_osm'

union['source'] = union.apply(classify, axis=1)

# Now merge same-class adjacent fragments
dissolved = union.dissolve(by='source')

# Now merge adjacent official-only and osm-only
# This is the tricky part
```

Hmm, geopandas overlay with how='union' returns fragments. For a shifted park:
- Fragment in official_only → in_off=True, in_osm=False
- Fragment in osm_only → in_off=False, in_osm=True
- Fragment in both → in_off=True, in_osm=True

These fragments are what I need. Now I need to:
1. Dissolve same-source fragments → 3 groups: official_only, osm_only, both
2. Merge official_only with adjacent osm_only

For step 2, I could check if official_only and osm_only geometries are adjacent (touch) and if so, union them.

But this is lossy because dissolving first loses adjacency info at finer granularity.

Alternative: don't dissolve first. Instead:
1. From the overlay, we have many small fragments each with a classification
2. Merge adjacent fragments that have classification "official_only" OR "osm_only" (i.e., all disagreement fragments)
3. But what source do we give the merged groups? If a group contains only official_only fragments → "parks_official". If only osm_only → "parks_osm". If both → "both"?

Actually yes! That's exactly the interpretation:
- A group with only official-only fragments → source = "parks_official"
- A group with only osm-only fragments → source = "parks_osm"
- A group with BOTH types of fragments → source = "both" (this is the "shifted park" case where gain and loss are adjacent)

And for areas where both agree, they're already classified as "both".

Let me implement this with a graph-based merging approach.

Actually, let me think about a simpler approach using shapely ops.

1. Dissolve both layers to single multipolygons
2. Compute the three sets: both, official_only, osm_only
3. For each geometry, get individual polygons
4. Build an adjacency graph between official_only polygons and osm_only polygons
5. Merge connected components

But this is quite complex. Let me try the overlay approach.

Actually, let me try a different, cleaner approach:

```python
from shapely.ops import unary_union, polygonize
from shapely.geometry import MultiPolygon, Polygon

off_dissolved = unary_union(off.geometry)
osm_dissolved = unary_union(osm.geometry)

# Get the boundary lines
boundaries = unary_union([
    off_dissolved.boundary,
    osm_dissolved.boundary
])

# Polygonize
polygons = list(polygonize(boundaries))

# For each polygon, classify
results = []
for poly in polygons:
    in_off = poly.representative_point().within(off_dissolved) or off_dissolved.contains(poly)
    in_osm = poly.representative_point().within(osm_dissolved) or osm_dissolved.contains(poly)
    if in_off and in_osm:
        src = 'both'
    elif in_off:
        src = 'parks_official'
    elif in_osm:
        src = 'parks_osm'
    else:
        continue  # shouldn't happen but skip
    results.append({'source': src, 'geometry': poly})
```

This gives me atomic fragments. Then I need to merge adjacent fragments that should be merged.

For merging, I need adjacency groups. Within same source, adjacent fragments should merge. Then for disagreement (official_only + osm_only), adjacent fragments regardless of source should also merge.

Let me think about this more carefully with the actual geometry.

Actually, let me use a simpler, more robust approach with geopandas unary_union and shapely:

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

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

# Dissolve each layer
off_union = unary_union(off.geometry)
osm_union = unary_union(osm.geometry)

# Compute the three regions
both = off_union.intersection(osm_union)
official_only = off_union.difference(osm_union)
osm_only = osm_union.difference(off_union)

# Now, official_only and osm_only might have multiple parts
# Merge adjacent parts from official_only and osm_only

# Convert to individual polygons
def get_polygons(geom):
    if geom.is_empty:
        return []
    if geom.geom_type == 'Polygon':
        return [geom]
    elif geom.geom_type == 'MultiPolygon':
        return list(geom.geoms)
    else:
        return []

off_polys = get_polygons(official_only)
osm_polys = get_polygons(osm_only)
both_polys = get_polygons(both)

# Classify each polygon
all_polys = []
for p in off_polys:
    all_polys.append({'source': 'parks_official', 'geometry': p})
for p in osm_polys:
    all_polys.append({'source': 'parks_osm', 'geometry': p})
for p in both_polys:
    all_polys.append({'source': 'both', 'geometry': p})

# Now merge adjacent disagreement patches
# Build adjacency graph
# ...
```

This approach is cleaner. The adjacency graph approach:

1. For each polygon, note its source
2. Build adjacency: two polygons are adjacent if they touch (share boundary)
3. Merge all polygons in the same connected component IF they are all disagreement (official_only or osm_only) or all same source

Wait, the rule should be:
- "both" polygons stay as is (or merge with adjacent "both" polygons)
- official_only and osm_only polygons that are adjacent should merge into one group

Let me think about the connected components:

Build a graph where nodes are polygons, edges connect adjacent polygons. Then find connected components. For each component:
- If all nodes are "both" → merge into one, source = "both"
- If all nodes are "parks_official" → merge into one, source = "parks_official"
- If all nodes are "parks_osm" → merge into one, source = "parks_osm"
- If mix of "parks_official" and "parks_osm" (and possibly "both") → merge into one, source = "both"

Wait, but "both" polygons shouldn't be merged with disagreement polygons unless they're adjacent. Actually, a "both" polygon adjacent to a disagreement polygon... hmm. Let me re-think.

The user says: "every patch claimed by exactly one side, with adjacent disagreement patches merged together"

So the merging is ONLY for disagreement patches (patches claimed by exactly one side). Both areas are separate.

But what if a "both" area is between an official_only and osm_only area? That seems unlikely from the geometry but possible.

Actually, I think the simplest correct approach is:

1. Calculate the overlay union → fragments with classification
2. For fragments classified as "both", dissolve them (source = "both")
3. For fragments classified as disagreement (official_only or osm_only), merge ALL adjacent disagreement fragments regardless of their specific source
4. For each merged disagreement group, determine its source:
   - If it contains only official_only fragments → "parks_official"
   - If it contains only osm_only fragments → "parks_osm"
   - If it contains both types → "both"

Let me implement this.

For the adjacency/merging approach, I'll use a simple graph-based method with shapely.

Let me write the full script now.

Actually, the polygonize approach might be cleaner since it gives atomic pieces and avoids issues with the overlay approach.

Let me write the complete script:

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

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

# CRS is EPSG:6539 (projected, meters)
crs = off.crs

# Dissolve each layer to a single (multi)polygon
off_union = unary_union(off.geometry)
osm_union = unary_union(osm.geometry)

# Get boundary lines and polygonize to get atomic fragments
boundaries = unary_union([
    off_union.boundary,
    osm_union.boundary
])
atomic_polys = list(polygonize(boundaries))

# Classify each atomic polygon
fragments = []
for poly in atomic_polys:
    centroid = poly.representative_point()
    
    in_off = off_union.contains(centroid) or (off_union.touches(poly) and off_union.contains(poly.representative_point()))
    in_osm = osm_union.contains(centroid) or (osm_union.touches(poly) and osm_union.contains(poly.representative_point()))
    
    # Actually, better approach: check if the polygon is mostly within the layer
    # Use a point that's definitely inside
    rep_point = poly.representative_point()
    in_off = off_union.contains(rep_point) or poly.within(off_union)
    in_osm = osm_union.contains(rep_point) or poly.within(osm_union)
    
    if in_off and in_osm:
        src = 'both'
    elif in_off:
        src = 'parks_official'
    elif in_osm:
        src = 'parks_osm'
    else:
        continue  # shouldn't happen
    
    fragments.append({'source': src, 'geometry': poly})

# Now merge adjacent fragments
# Build adjacency graph
G = nx.Graph()
for i, f in enumerate(fragments):
    G.add_node(i, **f)

# Add edges for touching polygons (only among disagreement)
for i in range(len(fragments)):
    for j in range(i+1, len(fragments)):
        if fragments[i]['geometry'].touches(fragments[j]['geometry']):
            G.add_edge(i, j)

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

# For each component, merge geometries
results = []
for comp in components:
    polys = [fragments[i]['geometry'] for i in comp]
    sources = set(fragments[i]['source'] for i in comp)
    
    merged = unary_union(polys)
    
    if 'both' in sources:
        # If any fragment is 'both', the whole group is 'both'
        src = 'both'
    elif len(sources) == 1:
        src = sources.pop()
    elif sources == {'parks_official', 'parks_osm'}:
        # Mixture of official-only and osm-only → 'both'
        src = 'both'
    else:
        src = 'both'  # fallback
    
    results.append({'source': src, 'geometry': merged})
```

Wait, this merges ALL adjacent fragments, including "both" with disagreement. That might not be desired. Let me separate "both" fragments:

The user wants:
1. Areas where both agree → "both" (these are NOT disagreement)
2. Areas of disagreement → patches claimed by exactly one side
3. Adjacent disagreement patches merged together

So "both" areas should NOT be merged with disagreement areas. They should stay as separate "both" groups.

Let me revise: 
1. Separate fragments into "both" and "disagreement" (official_only + osm_only)
2. Within "both", merge adjacent same-source fragments (all "both")
3. Within disagreement, merge ALL adjacent fragments (regardless of official_only/osm_only)
4. Assess source for each merged disagreement group

Let me rewrite:

```python
# Separate
both_frags = [f for f in fragments if f['source'] == 'both']
disc_frags = [f for f in fragments if f['source'] != 'both']

# Build graph for both fragments
G_both = nx.Graph()
for i, f in enumerate(both_frags):
    G_both.add_node(i, **f)
for i in range(len(both_frags)):
    for j in range(i+1, len(both_frags)):
        if both_frags[i]['geometry'].touches(both_frags[j]['geometry']):
            G_both.add_edge(i, j)

# Build graph for disagreement fragments
G_disc = nx.Graph()
for i, f in enumerate(disc_frags):
    G_disc.add_node(i, **f)
for i in range(len(disc_frags)):
    for j in range(i+1, len(disc_frags)):
        if disc_frags[i]['geometry'].touches(disc_frags[j]['geometry']):
            G_disc.add_edge(i, j)

# Process both groups
results = []
for comp in nx.connected_components(G_both):
    polys = [both_frags[i]['geometry'] for i in comp]
    merged = unary_union(polys)
    results.append({'source': 'both', 'geometry': merged})

# Process disagreement groups
for comp in nx.connected_components(G_disc):
    polys = [disc_frags[i]['geometry'] for i in comp]
    sources = set(disc_frags[i]['source'] for i in comp)
    merged = unary_union(polys)
    
    if len(sources) == 1:
        src = sources.pop()
    else:
        src = 'both'
    
    results.append({'source': src, 'geometry': merged})
```

This looks good. Now for label anchors, I need a point strictly inside each geometry. I can use `representative_point()` which guarantees a point within the geometry.

Let me write the full script now.

One concern: the polygonization might produce very many small fragments, and the O(n²) adjacency check could be slow. Let me use a spatial index to speed it up.

Let me also consider: the `polygonize` step requires the boundary lines to form a valid planar graph. With complex geometries, there might be precision issues. I should use a small buffer on the boundary before polygonizing, or use `snap` to snap geometries together.

Actually, for robustness, let me buffer slightly before polygonizing, or use a different approach.

Let me think of an alternative approach that avoids polygonize entirely:

1. Compute `both = intersection(off_union, osm_union)`
2. Compute `official_only = difference(off_union, osm_union)`
3. Compute `osm_only = difference(osm_union, off_union)`
4. Each of these is a (multi)polygon
5. Extract individual polygons from each
6. Build adjacency graph among all disagreement polygons (official_only + osm_only)
7. Merge connected components
8. Assess source for each merged group

This avoids polygonize entirely and uses shapely operations directly, which are more robust. The only tricky part is figuring out adjacency between polygons that might share very small boundaries.

Let me go with this approach.

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

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

# Dissolve each layer
off_union = unary_union(off.geometry)
osm_union = unary_union(osm.geometry)

# Make valid (fix any geometry issues)
off_union = make_valid(off_union) if hasattr(off_union, 'is_valid') and not off_union.is_valid else off_union
osm_union = make_valid(osm_union) if hasattr(osm_union, 'is_valid') and not osm_union.is_valid else osm_union

# Compute three regions
both_region = off_union.intersection(osm_union)
off_only_region = off_union.difference(osm_union)
osm_only_region = osm_union.difference(off_union)

# Extract individual polygons
def extract_polygons(geom):
    """Extract individual Polygon objects from a geometry"""
    if geom is None or geom.is_empty:
        return []
    geom = make_valid(geom)
    if geom.geom_type == 'Polygon':
        return [geom]
    elif geom.geom_type == 'MultiPolygon':
        return list(geom.geoms)
    elif geom.geom_type == 'GeometryCollection':
        result = []
        for g in geom.geoms:
            result.extend(extract_polygons(g))
        return result
    return []

both_polys = extract_polygons(both_region)
off_only_polys = extract_polygons(off_only_region)
osm_only_polys = extract_polygons(osm_only_region)

# Create fragment list
fragments = []
for p in both_polys:
    fragments.append({'source': 'both', 'geometry': p})
for p in off_only_polys:
    fragments.append({'source': 'parks_official', 'geometry': p})
for p in osm_only_polys:
    fragments.append({'source': 'parks_osm', 'geometry': p})
```

Now for adjacency: `touches()` in shapely checks if geometries share boundary but don't overlap. For the merged disagreement case, we want to know if official_only and osm_only polygons are adjacent (share a boundary).

However, there might be a tiny gap between them due to precision issues. I should use a small tolerance (e.g., buffer slightly).

Let me use an alternative adjacency check: two polygons are adjacent if their distance is 0 and they don't overlap (or the intersection is a line/point, not an area).

Actually, `touches` in shapely should work for this scenario. Since these geometries all came from the same dissolved union, their boundaries should line up exactly.

Let me proceed with the implementation and test it.

Let me also think about performance. With maybe hundreds of polygons after dissolving, the O(n²) check should be fine.

Let me write the complete script now.

One more thought: for the "both" areas, do we also need to consider that they might be formed by merging adjacent "both" fragments? Yes, I think so.

Let me also make sure `networkx` is installed.

Let me think about the task:

1. We have two layers: `parks_official` and `parks_osm` in EPSG:6539 (NY state plane, metric).
2. We need to find areas claimed by exactly one side - i.e., symmetric difference between the two layers.
3. But we need adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature. This means: a park that shifted position would have areas where Official has it but OSM doesn't (loss), and areas where OSM has it but Official doesn't (gain). These two areas are adjacent but belong to different sides. The user wants them merged into one feature.

Wait, let me re-read: "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 in Official but not OSM (official-only)
- Areas in OSM but not Official (osm-only)

If a park boundary shifted, the old area is "official-only" and the new area is "osm-only". These two might be adjacent (touching). The user wants them merged into a single MultiPolygon with `source = 'both'`? No wait...

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

Hmm, so the source attribute tells what was the original source of these disagreement patches? Let me think again.

The disagreement regions are:
1. Areas in `parks_official` but not in `parks_osm` → these are `parks_official` source
2. Areas in `parks_osm` but not in `parks_official` → these are `parks_osm` source
3. Areas in both → `both` source (but this doesn't show disagreement)

Wait, the user says "every patch claimed by exactly one side". So we're only looking at:
- Official \ OSM (official-only)
- OSM \ Official (osm-only)

And 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 an official-only patch and an osm-only patch are adjacent (touching), they should be merged into one feature. What would the source be for that merged feature? The merged feature contains both official-only and osm-only areas, so it's disputed by both sides. The source would be... `both`?

Actually, I think the logic is:
- Each merged group has a source. If the group contains only official-only areas, source = "parks_official". If only osm-only areas, source = "parks_osm". If both types (i.e., they touched and got merged), source = "both".

Let me plan the algorithm:

1. Compute `official_only = official.difference(osm_union)` - areas in official but not in osm
2. Compute `osm_only = osm.difference(official_union)` - areas in osm but not in official
3. Collect all these polygons (from both sets), assign them a label "official" or "osm"
4. Dissolve/merge polygons that are adjacent (touch) into groups, regardless of their label
5. For each merged group, determine source:
   - If all constituent parts are from official_only → source = "parks_official"
   - If all constituent parts are from osm_only → source = "parks_osm"
   - If mixed → source = "both"
6. Generate label anchors: point inside each merged group's geometry

Let me implement this properly.

The key challenge is Step 4: merging adjacent disagreement patches. I'll use a graph-based approach where nodes are individual polygons and edges connect touching polygons. Then each connected component becomes one merged MultiPolygon.

Let me write the solution.

Wait, I have an important concern. The problem says "adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two." This means I need to merge disagreement patches that are adjacent, but also a shifted park means the official-only and OSM-only areas may be adjacent to each other (touching). 

But also, what about areas that are in both layers? I need to think about this more carefully.

"every patch claimed by exactly one side" - this means the symmetric difference, i.e., areas in official but not OSM (official-only) plus areas in OSM but not official (OSM-only). These are the disagreement patches.

"with adjacent disagreement patches merged together" - if an official-only polygon touches an OSM-only polygon (because the park shifted), they should be merged.

But what about the "both" source? The spec says `source` attribute with value one of `parks_official`, `parks_osm`, or `both`. So if a merged group contains only official-only pieces → source = "parks_official". Only OSM-only pieces → "parks_osm". Mixed → "both".

Now, one thing I should also consider: there might be areas that are claimed by BOTH sides (intersection of the two layers). The problem says "every patch claimed by exactly one side" - so the intersection (both) patches are NOT included in the disagreement groups. The `both` source label is only for merged groups that contain a mix of official-only and OSM-only pieces that are adjacent.

Wait, let me re-read: "For each group I also need a label-anchor point that falls strictly inside the group's geometry."

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

So source "both" means a merged group has both official-only and osm-only patches.

Let me also think about whether I should include areas that are in both layers but are surrounded by or adjacent to disagreement areas. The spec says "every patch claimed by exactly one side" — so areas in both are not disagreement patches. Those areas should simply not appear. Only disagreement patches (symmetric difference) should be part of the groups.

This looks right. Let me run it.

One more thought: `touches` may not capture all cases if boundaries are slightly misaligned (topological errors). Using a tiny buffer could help. Let me add that.

Actually, `touches` in Shapely checks if the geometries share at least one boundary point but don't have interior intersections. Since these are from two different sources with potentially slightly different boundaries, `touches` might miss some adjacencies that are "almost touching." A buffer of e.g. 1 cm could help.

Let me improve the script.

Let me also verify the GeoJSON files are valid by checking their raw content structure, and ensure the `source` values exactly match what was asked for: `parks_official`, `parks_osm`, or `both`. Let me verify those strings are correct.

The user asked for:
- `source` attribute whose value is exactly one of `parks_official`, `parks_osm`, or `both`
- These should match the input GPKG layer names

The input layers are `parks_official` and `parks_osm`. Let me verify the values match exactly.

Both GeoJSON files have been created and verified. Let me confirm all the requirements are met:

1. ✅ `parks_disagreement.geojson` - FeatureCollection with `source` attribute having exactly `parks_official`, `parks_osm`, or `both`
2. ✅ `park_label_anchors.geojson` - FeatureCollection with same `source` attribute, one anchor per disagreement group in matching order
3. ✅ Disagreement patches merged by adjacency (connected components graph)
4. ✅ Each merged group is a MultiPolygon
5. ✅ Each anchor point falls strictly inside its group's geometry
6. ✅ Both files have 46 features with matching order and sources

Let me verify the files are on disk.

Everything looks good. Let me also verify the source attribute values are exactly the required strings.

The user requested:
- `source` attribute with values exactly one of `parks_official`, `parks_osm`, or `both`
- These match the input GPKG layer names exactly

Let me check that the output doesn't have a `name` or `crs` field in the GeoJSON (that's a pyogrio thing). Actually, looking at the output, it has `"name": "parks_disagreement"` and a `crs` object. While these are valid GeoJSON, they're not standard - but the GeoJSON specification actually says the only required members are "type" and "features". The CRS and name are allowed extensions. Let me check if these cause issues.

Actually, looking at the output more carefully, pyogrio adds `"name"` and `"crs"` at the top-level of the FeatureCollection. The standard GeoJSON spec (RFC 7946) says the CRS is always WGS84 and should not be explicitly given, but older implementations include it. Since we're working in EPSG:6539 (a projected CRS, not WGS84), pyogrio is correctly including the CRS info.

However, the user didn't ask for WGS84 output. They asked for GeoJSON. The data is in EPSG:6539 which is a projected CRS. pyogrio writes the CRS in the old GeoJSON crs format. Some consumers might expect WGS84 only. But since the task says to write GeoJSON with the source column, and the input data is in EPSG:6539, I think keeping the native CRS is correct.

Let me also verify the file contents are valid by loading them.

All checks passed. Let me do a final summary.