geo-l2-nyc-park-symdiff

Rerun

Let me start by understanding the problem:

1. We have a GeoPackage `nyc_parks.gpkg` with official parks and OSM-derived parks layers.
2. I need to find areas where exactly one side claims a park (disagreement areas).
3. Adjacent disagreement patches from the same source should be merged together.
4. Each merged group should be a MultiPolygon with a `source` attribute.
5. I need label-anchor points inside each group's geometry.

Let me first explore the GPKG file to understand the layers and their structure.

Let me understand the problem:

1. We have two layers in the GPKG: `parks_official` and `parks_osm`
2. Both have the same CRS (EPSG:6539) and both contain MultiPolygon geometries
3. They both have a `park_id` column - each park appears in both layers (same IDs), but the geometries may differ
4. I need to find areas that are in one but not the other
5. Adjacent disagreement patches should be merged together (so if park A shifted, its "gain" on one side and "loss" on the other should form a single feature if they're adjacent)
6. I need to output:
   - `parks_disagreement.geojson` with MultiPolygon features, each with a `source` attribute that is `parks_official`, `parks_osm`, or `both`
   - `park_label_anchors.geojson` with the same source attribute and one anchor point per feature, strictly inside the geometry

Let me think about this more carefully:

The task says "every patch claimed by exactly one side" - so I need areas that are in one but not the other. But then "with adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two."

Wait, let me re-read: "I need every patch claimed by exactly one side" - this means areas where only parks_official claims them OR only parks_osm claims them, not both.

"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 boundary shifted, there will be a "loss" (area in official but not OSM) and a "gain" (area in OSM but not official). If these two areas are adjacent (touching), they should be merged into a single feature.

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

Wait, so the source attribute can be `both` too? That doesn't make sense with "every patch claimed by exactly one side"... Let me re-read.

"I need every patch claimed by exactly one side" - This describes the disagreement patches. But then it says source can be `both` too? 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` (those exact strings, matching the input GPKG layer names)"

Hmm. So the source attribute for disagreement groups... but `both` would be areas claimed by both sides. That's agreement, not disagreement.

Wait, maybe the instruction is: Produce ALL areas (not just disagreement), categorized by source. Let me re-read the full first paragraph:

"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 actually "every patch claimed by exactly one side" = disagreement. The source values are `parks_official` or `parks_osm`. But the instruction also says `both` is a valid value. 

Hmm, maybe the thinking is: if we take all parks features and dissolve/union them, then some areas are:
- Only in official → source = `parks_official`
- Only in OSM → source = `parks_osm`
- In both → source = `both`

And adjacent patches of disagreement (official-only adjacent to osm-only) get merged into one feature with... what source? That doesn't make sense.

Re-reading again more carefully:

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

So if there's a shifted park, the area that was in official but not in OSM (loss) is adjacent to the area that's in OSM but not in official (gain). These two adjacent patches should be merged into ONE feature. But what source would that merged feature have? It can't be `parks_official` or `parks_osm` since it contains both.

Wait, maybe the intent is:
- Areas in both official AND OSM → source = `both`
- Areas only in official → source = `parks_official`  
- Areas only in OSM → source = `parks_osm`
- Gaps (areas not in either) → not included?
- BUT, if an area only in official is adjacent to an area only in OSM, they get merged into a single feature with source = ??? 

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

Maybe the output is supposed to be a complete partition of space into:
- both: areas covered by both
- parks_official: areas only in official
- parks_osm: areas only in OSM
- (areas not in either are excluded)

But then adjacent disagreement patches... If an official-only area touches an OSM-only area, they should be merged. But then what's the source attribute? 

Actually, re-reading: "every patch claimed by exactly one side" - maybe this means I take the union of all official and OSM parks, clip to areas claimed by exactly one source (symmetric difference), then dissolve/merge adjacent features together.

Wait, "claimed by exactly one side" - so not by both. That's 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 park at location A, and OSM has the same park shifted to location B (which partially overlaps), the symmetric difference would be:
- Area A\B (in official only)
- Area B\A (in OSM only)

If A and B overlap partially, A\B and B\A are NOT adjacent - they're separated by the overlap area. So they wouldn't be merged.

But if the park is completely shifted - no overlap - then A\B = A and B\A = B, and they might be adjacent (touching at a point or edge?). If a park moved from A to B and A and B touch, then they'd be merged.

Actually, I think the more common case is: a park boundary is slightly different between official and OSM. So the symmetric difference gives two adjacent polygons (the gain and loss areas) that share a boundary. These should be merged into one feature.

But wait - if they're from different sources (one official-only, one OSM-only), what source does the merged feature get?

Let me just re-read one more time:

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

OK so maybe the approach is:
1. Compute the full overlay (identify areas that are in official only, OSM only, or both)
2. For areas that are in both, source = `both`
3. For areas only in one, source is that layer
4. Now, merge adjacent features only if they are both disagreement (i.e., `parks_official` touching `parks_osm`) 

But what source does the merged feature get? Maybe the rule is: if both disagreement types are present in a merged group, it's... hmm.

Actually, maybe the intent is simpler: Just compute the symmetric difference (areas in exactly one source), then dissolve/st_buffer(0) to merge adjacent ones. The source attribute could be... well, each part of the symmetric difference is either official-only or OSM-only. But if we merge adjacent ones, we lose that distinction.

Let me try a different interpretation: Maybe the task wants:
1. Take all official parks, dissolve → unified official coverage
2. Take all OSM parks, dissolve → unified OSM coverage
3. Compute the overlay: 
   - official only (official \ OSM) → source = `parks_official`
   - OSM only (OSM \ official) → source = `parks_osm`
   - Both (official ∩ OSM) → source = `both`
4. Now, for any `parks_official` feature adjacent to a `parks_osm` feature, merge them into one MultiPolygon
   - But what source? Maybe the instruction doesn't mean to assign source per merged group but rather...

Actually wait, re-reading one last time: "One MultiPolygon per merged group with a `source` attribute whose value is exactly one of..."

I think maybe the instruction means: for each merged group, determine its source. If a merged group is formed by merging an official-only and an OSM-only area, then... the source could be `both` because it contains material from both? That would be confusing though.

OR, maybe the source is determined by the majority or... 

Let me try yet another interpretation. Maybe the task is:

1. Compute symmetric difference (areas in exactly one of the two sources)
2. Merge adjacent disagreement patches (so if A is in official only and B is in OSM only, and they touch, merge them)
3. But also keep the `both` areas separately
4. Each feature has a source attribute

So the full output has three categories:
- `both`: areas covered by both layers (the intersection)
- `parks_official`: areas only in official, and NOT adjacent to any OSM-only area... OR areas that result from merging official-only with adjacent OSM-only (but what source?)

This is getting complicated. Let me just go with the simplest literal interpretation:

"every patch claimed by exactly one side" → this is the symmetric difference (official Δ OSM)

"with adjacent disagreement patches merged together" → dissolve/merge touching features within the symmetric difference

"One MultiPolygon per merged group with a `source` attribute" → wait, but within the symmetric difference there are both official-only and OSM-only patches. If we merge them, we lose the distinction...

OK I think I'm overcomplicating this. Let me try this approach:

1. Compute the full overlay/intersection of all features from both layers
2. Get three categories:
   - In official only (official \ OSM)
   - In OSM only (OSM \ official)
   - In both (official ∩ OSM)
3. Take the features in "official only" and "OSM only" (disagreement patches)
4. Merge any adjacent disagreement patches (official-only touching OSM-only) into single features
5. For each merged group (which could be a mix of official-only and OSM-only parts), assign source as... hmm.

Actually, you know what, let me just interpret it as:

The output should have three types of features:
- `both`: areas where both layers agree (intersection area)
- `parks_official`: areas only in the official layer
- `parks_osm`: areas only in the OSM layer

Then merge adjacent `parks_official` and `parks_osm` features together. For the merged features, the source... maybe just keep the source as whichever is dominant?

Actually, re-reading again: "every patch claimed by exactly one side" - this specifically says disagreement. So only `parks_official` and `parks_osm` source values would apply.

"with adjacent disagreement patches merged together" - so an official-only patch touching an OSM-only patch should be merged. But what source?

Maybe the answer is: the merged feature's source is... both? Since it contains material from both?

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

I think the three values cover all cases:
- If the merged group consists only of official-only area → `parks_official`
- If the merged group consists only of OSM-only area → `parks_osm`
- If the merged group consists of both (because adjacent official-only and OSM-only were merged) → `both`

And also, `both` would include the true intersection areas (where both agree).

Let me just go with this approach. Let me code it up.

Actually, wait. Let me re-read the exact wording 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 it's specifically about disagreement patches. Areas where both agree (official ∩ OSM) are not disagreement, so they wouldn't be included? But the source attribute includes `both`...

Hmm, maybe the intent is to produce ALL areas (not just disagreement), categorized.

Let me just go with: produce all areas broken into three categories (official-only, OSM-only, both), merge adjacent official-only with OSM-only into `both`, and output.

Actually, I think the cleanest interpretation is:

1. Take the union of all official and OSM park features
2. For each part of the union, determine if it's in official only, OSM only, or both
3. Merge adjacent `parks_official` and `parks_osm` parts together (forming a feature with source = `both`)
4. Keep `both` parts as-is with source = `both`
5. Keep isolated `parks_official` or `parks_osm` parts (not adjacent to the other) with their respective source

Wait but that would mean a merged official-only + OSM-only feature gets source `both`, and a genuinely overlapped area also gets source `both`. That makes sense actually.

But then there's a subtlety: what if an official-only area touches another official-only area? Those should stay as separate features? Or merged?

I think the merging only applies to disagreement patches of DIFFERENT sources touching each other. Not same-source touching.

Let me just code it up step by step.

Actually, the simplest approach:

1. Dissolve all official parks → unified coverage
2. Dissolve all OSM parks → unified coverage
3. Compute:
   a. official_and_osm = intersection (both)
   b. official_only = official - osm (difference)
   c. osm_only = osm - official (difference)
4. Now, merge adjacent features: any official_only polygon touching an osm_only polygon gets merged
5. The merged group gets source = `both` (since it contains both types)
6. Isolated official_only polygons → source = `parks_official`
7. Isolated osm_only polygons → source = `parks_osm`
8. The intersection (official_and_osm) → source = `both`

Actually, the intersection areas are likely bigger contiguous areas. Let me think about the merge logic more carefully.

Maybe a better approach:

1. Compute the overlay of the two dissolved coverages
2. This gives us parts classified by their source presence
3. Now merge (dissolve) all parts with the same classification that are adjacent
4. Then, for `parks_official` parts adjacent to `parks_osm` parts, merge them together as `both`

Actually, let me use a different, cleaner approach:

1. Dissolve official → off
2. Dissolve OSM → osm
3. Compute intersection: both_geom = off.intersection(osm)
4. Compute official_only: off_geom.difference(osm)
5. Compute osm_only: osm_geom.difference(off)

Now, for official_only and osm_only, they might consist of multiple disjoint polygons. Any official_only polygon that touches an osm_only polygon should be merged.

Let me implement this.

Step 1: Read both layers, dissolve each
Step 2: Compute the three categories
Step 3: Merge adjacent official_only and osm_only polygons

For the merging, I can:
- Take all official_only polygons and osm_only polygons
- Use a spatial join / adjacency check to find touching pairs
- Merge them together

Or, simpler approach:
- Combine official_only and osm_only into a single GeoSeries
- Buffer(0) to dissolve adjacent ones
- Then for each resulting merged polygon, determine the source

Let me think about the source determination for merged features:
- If the merged polygon contains area from both official_only and osm_only, source = `both`
- If it only contains area from official_only, source = `parks_official`
- If it only contains area from osm_only, source = `parks_osm`

But after dissolving, how do I know which source contributed? I can use the original classification before dissolving.

Actually, let me try a different approach entirely:

1. For each park_id that appears in both layers, get the official geometry and OSM geometry
2. Compute the differences
3. Merge adjacent differences across parks

No, that's per-park which might not match the "dissolve all" approach.

Let me think about what makes sense spatially. If we have parks A and B, where:
- Official: A1 (park A), B1 (park B)  
- OSM: A2 (park A, slightly different), B2 (park B, same as B1)

The symmetric difference would be:
- For park A: A1\A2 (official only), A2\A1 (OSM only)
- For park B: nothing (they agree)

If A1\A2 and A2\A1 are adjacent (they would be, since A1 and A2 overlap), they should be merged.

If park B's official and OSM match exactly, then B is in the "both" category.

OK let me just code this up.

My plan:
1. Read both layers
2. Dissolve each to a single MultiPolygon geometry
3. Compute:
   - both_geom = official.intersection(osm)
   - off_only = official.difference(osm)
   - osm_only = osm.difference(official)
4. Now, to merge adjacent off_only and osm_only:
   - Explode off_only and osm_only into individual polygons
   - Find which ones touch
   - Group them using graph connectivity
   - For each group, create a MultiPolygon
   - Determine source: if group contains only off_only parts → `parks_official`; only osm_only parts → `parks_osm`; both → `both`
5. Also include both_geom parts as source = `both`
6. For each final feature, compute a label anchor point (point that's strictly inside)

Let me code this up.

Wait, actually I realize there's a subtlety with the "both" geometry. If I just take the intersection of the two dissolved coverages, I get a single (possibly multipart) geometry. But adjacent official-only and OSM-only merged features would also get source = `both`. So I need to combine these two sources of `both` features.

Actually, let me simplify:

1. Compute the full overlay partition:
   - In both: official ∩ osm
   - Official only: official - osm
   - OSM only: osm - official

2. Merge adjacent official-only and OSM-only polygons together.

3. The merged groups get source = `both`.
4. Non-merged official-only stay `parks_official`, non-merged OSM-only stay `parks_osm`.
5. The original intersection (from step 1) is source = `both`.

Let me code this up properly.

Actually, I think there may be a cleaner way. Let me use unary_union and overlay operations.

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

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

# Compute categories
both_geom = off_union.intersection(osm_union)
off_only_geom = off_union.difference(osm_union)
osm_only_geom = osm_union.difference(off_union)
```

Now I need to work with individual polygons, not MultiPolygons.

Let me use a graph-based approach to merge adjacent polygons.

```python
from shapely.geometry import Polygon, MultiPolygon

def explode(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_only_polys = explode(off_only_geom)
osm_only_polys = explode(osm_only_geom)
both_polys = explode(both_geom)
```

Now for each off_only polygon and osm_only polygon, I need to check adjacency and merge.

For the label anchor point: I need a point strictly inside the geometry. I can use `representative_point()` which guarantees being inside the geometry. Or I could use `point_on_surface()` which also guarantees being inside.

Wait, `representative_point()` is a method that returns a point guaranteed to be within the geometry. Let me use that.

OK let me write the full solution.

For the merging logic:
1. Create a graph where each polygon (from off_only and osm_only) is a node
2. Edges connect polygons that touch (share a boundary)
3. Each connected component becomes a merged feature
4. Determine source: if any off_only in the component, and any osm_only, source = `both`; otherwise, source is the category

Wait, but `both` is also used for the intersection areas. So I need to combine:
- Merged groups from off_only + osm_only adjacency → `both`
- Non-merged off_only → `parks_official`
- Non-merged osm_only → `parks_osm`
- Intersection areas (both_geom parts) → `both`

But a merged group (off_only + osm_only adjacent) could potentially be adjacent to a `both` area... Should I merge those too? The instruction says "adjacent disagreement patches" - so disagreement patches only. The `both` areas are agreement, not disagreement.

OK so the final output has:
- `both` areas (from intersection) 
- merged disagreement groups (which could be `parks_official`, `parks_osm`, or `both` depending on composition)
- isolated disagreement patches (non-merged)

Let me write the code.

Actually wait, let me reconsider. If I merge off_only and osm_only polygons that touch, the result might have source `both`. But then a `both` area from the intersection might be adjacent to this merged `both` group. Should they be merged? The instruction says "adjacent disagreement patches merged together" - the intersection areas aren't disagreement patches. So I think we keep them separate.

Let me code this up.

For finding adjacent polygons efficiently, I can use a spatial index or just brute force for reasonable numbers of polygons.

Let me write the script.

Both layers have `park_id` which looks like a UUID. The same park_id appears in both layers for matching parks. The task is about finding disagreement patches — areas where the official and OSM layers differ.

Let me think about this:

1. Each park has a `park_id`. Parks that appear in both layers may have different geometries (shifts, boundary differences).
2. For each park, I need to find:
   - Area only in official (official - osm)
   - Area only in osm (osm - official)
   - Area in both (intersection)
3. "Disagreement patches" are areas where exactly one side claims the park. So I need the symmetric difference: `official.symmetric_difference(osm)` for each matching park pair.
4. "Adjacent disagreement patches merged together" — if two different parks have disagreement patches that are adjacent (touching), they should be merged together into one MultiPolygon.
5. Then label each merged group with source = "parks_official" if the merged group only contains official-only areas, "parks_osm" if only OSM-only areas, or "both" if it contains both types.

Wait, re-reading: "every patch claimed by exactly one side" - so we take the symmetric difference per park. Then merge adjacent disagreement patches across parks. Then classify each merged group by whether its constituent original patches came from official side, osm side, or both.

Let me think about this more carefully.

For each park_id that exists in both layers:
- official_only = official_geom - osm_geom (difference)
- osm_only = osm_geom - official_geom (difference)

For a park only in official:
- official_only = official_geom
- osm_only = empty

For a park only in osm:
- osm_only = osm_geom
- official_only = empty

Now, these patches (official_only and osm_only) are the "disagreement patches". Each patch has a source:
- official_only patches have source "parks_official"
- osm_only patches have source "parks_osm"

Two patches that are adjacent (touching) should be merged together. If a merged group contains both official_only and osm_only patches, the source becomes "both".

Let me code this up.

Actually, 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 if a park shifted, the area it lost (official - osm) and the area it gained (osm - official) are adjacent (they touch at the shifted boundary). Those should be merged. So the source becomes "both" for that merged group.

But also, disagreement patches from different parks that are adjacent should also be merged together.

Let me implement:

1. Compute symmetric difference per park (or per park_id pair).
2. Collect all "official_only" and "osm_only" polygons.
3. Dissolve/merge adjacent ones.
4. For each merged group, determine source.

Let me write this step by step.

Actually, I think the simplest approach:
1. For each park_id (union of both sets), get official geom (if exists) and osm geom (if exists).
2. Compute:
   - both: intersection
   - official_only: official.difference(osm) if both exist, else official_geom if only official
   - osm_only: osm.difference(official) if both exist, else osm_geom if only osm
3. Collect all official_only and osm_only polygons into one GeoDataFrame with a source column.
4. Explode MultiPolygons into individual polygons so we can dissolve/merge them.
5. Union/dissolve polygons that are adjacent (touch). We can do this by dissolving with ST_Union (unary_union) which merges adjacent polygons, then polygonize the result.

Actually, geopandas dissolve with `dissolve(by=None)` will merge all touching features into a MultiPolygon. But we want to only merge adjacent ones, not all of them.

Better approach: use `unary_union` on the whole collection, which merges all adjacent/touching polygons into multipolygons. Then each resulting geometry is a merged group.

Wait, `unary_union` merges all touching/overlapping geometries into a single geometry. So if we pass all polygons, it will merge all that touch each other, which is what we want.

But we need to track the source. Let me think...

After unary_union, we lose the mapping from original polygons to merged groups. We need to figure out which polygons belong to which merged group.

Alternative approach:
1. Compute `unary_union` of all disagreement polygons to get a merged collection.
2. For each merged geometry in the result, find which original polygons intersect it, and determine the source from those.

Let me code this up.

Actually, there's a clean approach using spatial join and dissolve:

1. Start with all disagreement patches (official_only and osm_only) in a single GDF with a `source` column.
2. Explode MultiPolygons to Polygons (singlepart).
3. Compute adjacency groups by doing a spatial union and then using `unary_union` to merge, then spatially joining back.

Let me think of a simpler way:

1. Compute `unary_union` of all disagreement polygons -> merged_multi.
   This yields a (Multi)Polygon where each polygon component is a merged group.
   
Wait, `unary_union` returns a single geometry (could be MultiPolygon or Polygon). If two groups are not adjacent, they'll be separate polygons in the MultiPolygon. But if they are adjacent, they'll be merged.

Hmm, but `unary_union` merges ALL overlapping/touching polygons in the input. So it would merge all official_only and osm_only patches that are adjacent, regardless of source.

Then the result is a single geometry (MultiPolygon or Polygon). But we need each contiguous merged group as a separate feature. So we can explode the result into individual polygons (if MultiPolygon, get each sub-polygon; if Polygon, just one).

Then for each resulting polygon, we need to find which original patches contributed to it. We can do a spatial overlay - for each merged polygon, check which original disagreement patches intersect it.

Let me code this step by step.

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

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

# Build a dictionary of park_id -> geometry
official_dict = dict(zip(official['park_id'], official['geometry']))
osm_dict = dict(zip(osm['park_id'], osm['geometry']))

all_ids = set(official_dict.keys()) | set(osm_dict.keys())

# Collect disagreement patches
patches = []  # list of (geometry, source)

for pid in all_ids:
    off_geom = official_dict.get(pid)
    osm_geom = osm_dict.get(pid)
    
    if off_geom is not None and osm_geom is not None:
        # Both exist - compute difference
        off_only = off_geom.difference(osm_geom)
        osm_only = osm_geom.difference(off_geom)
        
        if not off_only.is_empty:
            patches.append((off_only, 'parks_official'))
        if not osm_only.is_empty:
            patches.append((osm_only, 'parks_osm'))
    elif off_geom is not None:
        # Only official
        patches.append((off_geom, 'parks_official'))
    elif osm_geom is not None:
        # Only osm
        patches.append((osm_geom, 'parks_osm'))

# Create GDF
patches_gdf = gpd.GeoDataFrame(patches, columns=['geometry', 'source'], crs=official.crs)
```

Now we have all patches. Next, merge adjacent ones.

```python
# Explode to single polygons (in case any are MultiPolygon with separate parts)
patches_exploded = patches_gdf.explode(index_parts=True).reset_index(drop=True)

# Compute unary union to merge touching polygons
merged_geom = unary_union(patches_exploded.geometry)

# Explode into individual polygons (each is a merged group)
if isinstance(merged_geom, Polygon):
    merged_polys = [merged_geom]
elif isinstance(merged_geom, MultiPolygon):
    merged_polys = list(merged_geom.geoms)

# For each merged group, find which original polygons intersect it
merged_features = []
for mp in merged_polys:
    # Find patches that intersect this merged polygon
    intersecting = patches_exploded[patches_exploded.intersects(mp)]
    sources = intersecting['source'].unique()
    
    if len(sources) == 1:
        source = sources[0]
    else:
        source = 'both'
    
    merged_features.append({'geometry': mp, 'source': source})

result_gdf = gpd.GeoDataFrame(merged_features, crs=official.crs)
```

Wait, but there's a subtlety. If `patches_exploded` contains polygons that are adjacent (touching), `unary_union` will merge them. But what about polygons from the same source that are adjacent to polygons from another source? They'll be merged too, and the source becomes "both". That's what we want.

But we need to be careful: `unary_union` might also merge non-adjacent but overlapping polygons. Do any disagreement patches overlap? They shouldn't since official_only and osm_only for the same park are disjoint, and for different parks they should also be disjoint (park boundaries don't overlap in the same dataset). But we should be safe.

Actually, wait. There's another issue. When I do `off_only = off_geom.difference(osm_geom)`, the result might be a MultiPolygon. When I explode it, individual parts that are adjacent to other patches will get merged. That's fine.

But what about label anchors? I need a point that falls strictly inside the group's geometry. I can use `representative_point()` which guarantees an interior point.

Let me also think about whether `unary_union` will actually merge touching polygons. In Shapely, `unary_union` merges polygons that share boundaries (touch) or overlap. So yes.

Let me write the full script.

Actually, I also need to handle the case where there's only one merged polygon (result is Polygon not MultiPolygon). And also the case where there are no disagreement patches at all (unlikely but possible).

Let me also think about edge cases: what if patches_exploded is empty? Then there are no disagreements, and both output files should be empty.

Let me write the full script now.

Also for the label anchors, I need to ensure the point is strictly inside. `representative_point()` returns a point that is guaranteed to be within the polygon. For MultiPolygon, it returns a point in one of the sub-polygons. So that works.

Let me also double-check: for a merged group that's a MultiPolygon (e.g., two separate areas that don't touch each other but both have source "parks_official"), the unary_union would NOT merge them (since they don't touch). So they'd be separate merged groups, each with source "parks_official". That's correct.

Wait, actually no. Let me re-read the problem: "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"

So the output is "One MultiPolygon per merged group". That means even if a merged group is a single polygon, we should output it as a MultiPolygon (containing just that one polygon). Let me make sure of that.

Looking at the source layer, `parks_official` geometries are MultiPolygon. So the output should also be MultiPolygon.

Let me also think: if two disagreement patches from different parks are adjacent (share a boundary), they get merged into one. That means the merged group is a MultiPolygon (if multiple parts) or a Polygon (if single part). The spec says "One MultiPolygon per merged group", so I should convert Polygon to MultiPolygon.

Let me write the full solution now.

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

# Build dictionaries
official_dict = dict(zip(official['park_id'], official.geometry))
osm_dict = dict(zip(osm['park_id'], osm.geometry))

all_ids = set(official_dict.keys()) | set(osm_dict.keys())
print(f"Total unique park IDs: {len(all_ids)}")

# Collect disagreement patches
patch_geoms = []
patch_sources = []

for pid in sorted(all_ids):
    off_geom = official_dict.get(pid)
    osm_geom = osm_dict.get(pid)
    
    if off_geom is not None and osm_geom is not None:
        off_only = off_geom.difference(osm_geom)
        osm_only = osm_geom.difference(off_geom)
        
        if not off_only.is_empty:
            patch_geoms.append(off_only)
            patch_sources.append('parks_official')
        if not osm_only.is_empty:
            patch_geoms.append(osm_only)
            patch_sources.append('parks_osm')
    elif off_geom is not None:
        patch_geoms.append(off_geom)
        patch_sources.append('parks_official')
    elif osm_geom is not None:
        patch_geoms.append(osm_geom)
        patch_sources.append('parks_osm')

print(f"Total disagreement patches: {len(patch_geoms)}")
```

Then explode and merge...

Actually, I realize there might be a simpler approach using `dissolveby` in geopandas. Let me think about how to group adjacent polygons.

The approach with unary_union is fine. Let me code it properly.

After unary_union, I get a single geometry (Polygon or MultiPolygon) where all touching polygons are merged. But wait - if there are two separate groups that don't touch each other, they'll both be in the same MultiPolygon result from unary_union. Then I need to extract each polygon from the MultiPolygon as a separate merged group.

So:
1. unary_union merges all touching polygons -> result is a (Multi)Polygon with potentially multiple disconnected components
2. Each disconnected component is a merged group
3. For each component, find which original patches intersect it, determine source

This is correct.

Let me write the full script.

Let me analyze the problem more carefully.

We have two layers:
1. `parks_official` - 1380 features
2. `parks_osm` - 1372 features

Both have a `park_id` column. 20 parks only in official, 12 only in OSM, 1360 in both.

The task is about finding disagreement patches:
- Where a park exists in one layer but not the other (20 + 12 = 32 such parks)
- Where geometry differs between layers for the same park (symmetric difference)

The requirement is: "every patch claimed by exactly one side" - meaning areas that are in one layer but not the other. When adjacent disagreement patches touch (e.g., a shifted park's gain and loss areas), they should be merged into a single MultiPolygon.

Let me think about this more carefully:

For each park, we can compute:
- `o - os` = area in official but not OSM → sourced from `parks_official`
- `os - o` = area in OSM but not official → sourced from `parks_osm`
- `o ∩ os` = area in both → sourced from `both`

Wait, re-reading the requirements: "every patch claimed by exactly one side" - so we're looking at areas that are only in official OR only in OSM. The areas in both are not "disagreement" patches. But the requirement says we need a `source` attribute which can be `parks_official`, `parks_osm`, or `both`. Hmm, so `both` areas could be part of the output too?

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 for each park:
- Official-only area = "claimed by official" → source = `parks_official`
- OSM-only area = "claimed by OSM" → source = `parks_osm`

But then it says "a shifted park's gain and loss form a single feature rather than two". If a park shifted, the area it lost (official only) and the area it gained (OSM only) are adjacent. They should be merged into one feature. But they have different sources! 

Wait, re-reading again: "adjacent disagreement patches merged together" - so if two patches with DIFFERENT sources touch, they should be merged into one feature? But what source would it have?

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"

So the groups can have source = `both` too. Let me think...

Maybe the approach is:
1. Compute all differences: official\osm (official-only), osm\official (osm-only), and official∩osm (both)
2. Union all these patches together
3. Merge touching patches into groups
4. For each group, determine its source based on what it contains

Actually, I think the approach is more nuanced. Let me consider:

For each park (matching by park_id):
- `official_only = official.geometry - osm.geometry` → source = parks_official
- `osm_only = osm.geometry - official.geometry` → source = parks_osm
- `both = official.geometry ∩ osm.geometry` → source = both

Then we merge all these patches. When patches touch/overlap, they get merged into groups. But overlapping won't happen since the three categories (official_only, osm_only, both) are disjoint by construction.

Wait, but when a park shifts, you get a patch of official_only (where the park was) and a patch of osm_only (where it moved to). These might be adjacent but not overlapping. The instruction says to merge adjacent disagreement patches together.

But what if an official_only patch touches an osm_only patch from a different park? Or what if they're from the same park? The instruction says "a shifted park's gain and loss form a single feature" - so patches from the same park's shift should be merged.

Hmm, but the source attribute... If we merge an official_only patch and an osm_only patch that are adjacent, what source value do we give? Let me think about this differently.

Actually, I think the key insight is: we should first generate all the "patches" (disjoint areas where sources disagree), then merge adjacent patches regardless of source, producing a "merged group". For each group, we determine the source based on... hmm.

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 a merged group could be:
- If it contains patches that are only in official → `parks_official`
- If it contains patches that are only in OSM → `parks_osm`
- If it contains patches that are in both → `both`

But if we merge an official_only patch with an osm_only patch, what would the source be? It could be... hmm, maybe the idea is different.

Let me think about it differently. Perhaps "adjacent disagreement patches" refers to patches of the SAME source type. Like, if two official-only patches touch (from different parks, or from the same park's multipolygon components), they should be merged.

Actually, I think the intended approach is:

1. For each park_id, compute:
   - The symmetric difference between official and OSM geometries (the area of disagreement)
   - The intersection (area of agreement)

2. But what we really want is the area of disagreement grouped by source.

Wait, actually, let me re-think based on "a shifted park's gain and loss form a single feature rather than two."

If a park's geometry shifted, the "official only" area (the old location) and "OSM only" area (the new location) might be adjacent but are distinctly different sources. Merging them doesn't make sense for the source attribute.

Unless... the approach is:
- For each park, compute the disagreement areas (official\osm and osm\official)
- Merge all touching disagreement areas across ALL parks, regardless of source
- Then for each merged group, determine if the area is from official, osm, or both

Let me try a different interpretation: "adjacent disagreement patches" could mean patches of the same type. Like if park A has an official-only area adjacent to park B's official-only area, they'd merge.

Actually, I think the most logical interpretation is:

For each park, compute:
1. `both` = intersection of official and OSM geometry
2. `official_only` = official minus OSM
3. `osm_only` = OSM minus official

Now, for each park, the `official_only` and `osm_only` patches might be adjacent (if the park shifted slightly). The instruction says to merge these adjacent disagreement patches together. So for a shifted park, the gain (osm_only) and loss (official_only) would merge into one feature.

But then what source? If it contains both official-only and osm-only areas... well, maybe the source would be... hmm.

Let me try yet another interpretation: Maybe we need to look at this from the "pixel" perspective. For each point in space:
- If it's only in official → source = parks_official
- If it's only in OSM → source = parks_osm  
- If it's in both → source = both

Then we dissolve adjacent regions with the same source. This gives us patches with a consistent source. But then the "adjacent disagreement patches merged together" thing... if a shifted park has official_only adjacent to osm_only, they have different sources and wouldn't be merged under this scheme.

OK let me just try the most straightforward approach and see what happens:

1. Union all official polygons
2. Union all OSM polygons
3. Compute three sets:
   - `official_only` = official_union - osm_union
   - `osm_only` = osm_union - official_union
   - `both` = official_union ∩ osm_union

4. Assign source:
   - official_only → parks_official
   - osm_only → parks_osm
   - both → both

5. Then... merge adjacent patches? But they're already dissolved by source.

Hmm, but the instruction says "adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature." A shifted park's gain (OSM only) and loss (Official only) are adjacent but have different sources. If we merge them, what source?

Unless the idea is that the merged feature encompasses both the gain and loss area, and it's considered a "disagreement" zone. In that case, the source would be... maybe we just dissolve ALL patches (official_only + osm_only + both) into groups of touching features, and for each group, if it has both official_only and osm_only, we assign based on a rule?

Actually wait - re-reading: "every patch claimed by exactly one side" - patches claimed by exactly one side = official_only or osm_only. And "adjacent disagreement patches merged together" - so official_only patches that touch osm_only patches get merged. Then the resulting group's source... could be determined by which side(s) contributed?

But the source values are only `parks_official`, `parks_osm`, or `both`. If a group contains both official_only and osm_only patches, it could be `both`.

Let me re-read the original: "I need every patch claimed by exactly one side" - meaning areas that are exclusively in one dataset. But then "adjacent disagreement patches merged together" - if I only have patches claimed by exactly one side, merging them is straightforward: merge touching ones.

But then what about `both`? "with a source attribute whose value is exactly one of parks_official, parks_osm, or both" - the `both` source seems to be for areas where both agree.

Wait, I think the output should contain ALL areas (both where they agree AND disagree), but organized by merged groups of disagreement patches.

Let me try this approach:

1. Compute for each park_id:
   - `both` = intersection
   - `official_only` = official - osm
   - `osm_only` = osm - official

2. Union all `both` geometries into a single `both` multipolygon
3. Union all `official_only` geometries
4. Union all `osm_only` geometries

5. Merge together ALL official_only and osm_only patches that touch (this is the "adjacent disagreement patches merged together" part)

6. For each merged group, determine source:
   - If it contains ONLY official_only → parks_official
   - If it contains ONLY osm_only → parks_osm
   - If it contains BOTH → both

7. Also add the `both` areas as separate features with source `both`.

Hmm, but step 6 seems weird - if I merge official_only and osm_only together, the source could be `both` since both sides claim parts of it.

Actually, I think the simplest sensible approach:

1. Compute three categories for each park:
   - `official_only` = official.geometry.difference(osm.geometry) → source = parks_official
   - `osm_only` = osm.geometry.difference(official.geometry) → source = parks_osm
   - `intersection` = official.geometry.intersection(osm.geometry) → source = both

2. Collect all these geometries with their source
3. Dissolve/merge touching geometries (using unary_union or dissolve with a touch predicate)
4. For each merged group, if it contains patches from a single source, use that source; if it contains multiple sources... hmm, this doesn't work cleanly.

Let me take a step back and think about what makes sense. The key operation is:

"For each park, compute the symmetric difference (disagreement). Then merge adjacent disagreement regions across parks."

So the process would be:
1. For each park, compute `disagreement = official.geometry.symmetric_difference(osm.geometry)` - this gives areas where the two disagree
2. Dissolve all disagreement areas across all parks into contiguous groups
3. For each contiguous group, determine the source:
   - If the area is only in official → parks_official
   - If the area is only in OSM → parks_osm
   - If... wait, the symmetric difference is by definition areas that are in one but not the other. So it can't be `both`.

Hmm, but the output spec says source can be `both`. So maybe `both` areas should also be included?

Let me try: Include ALL areas but grouped such that adjacent disagreement patches are merged.

1. For each park, create three patches:
   - `official_only` (in official not osm) → source = parks_official
   - `osm_only` (in osm not official) → source = parks_osm
   - `both` (in both) → source = both

2. Merge ALL patches that touch each other into groups

3. For each group, determine its source:
   - If all patches in group have same source → use that source
   - If patches have mixed sources → use `both`

Actually, this makes sense! If an official_only patch touches an osm_only patch (from a shifted park), they get merged into one group with source `both`, because both sides have some claim there.

And if a `both` patch touches an `official_only` patch (e.g., a park that partially agrees and partially disagrees), they'd merge into one group with source `both`.

This interpretation seems most aligned with the task description. Let me implement it.

Wait, but the `both` patches are huge - they're the intersection areas. If we merge all touching `both` with `official_only`, we'd get enormous merged groups that represent entire parks. That doesn't seem right.

Let me reconsider. The instruction says "every patch claimed by exactly one side" - this is the key. We ONLY care about patches where there's disagreement. The `both` areas are not disagreement.

"with adjacent disagreement patches merged together" - so we take disagreement patches (official_only and osm_only) and merge touching ones. The result should have patches claimed by exactly one side, merged when adjacent.

But then the source values... if we merge an official_only patch with an osm_only patch, the merged group contains areas claimed by both sides, so its source would be `both`.

And if we have an isolated official_only patch with no touching osm_only, its source is `parks_official`.

And if we have an isolated osm_only patch with no touching official_only, its source is `parks_osm`.

This seems right! Let me implement this.

So the algorithm:

1. For each park_id, compute:
   - official_only = official.geometry - osm.geometry (for that park)
   - osm_only = osm.geometry - official.geometry (for that park)

2. Collect all official_only and osm_only patches (as separate features with their source)

3. Merge all features that touch (using unary_union or similar, then explode and deduplicate... actually use dissolve with a touch predicate)

Let me use shapely ops unary_union with the collection to merge touching features.

Actually, a cleaner approach:
- Create a GeoDataFrame with all the patches and their source
- Use `dissolve` with `by` based on connectivity... hmm, how to dissolve touching features?

Best approach: 
1. Create a GeoSeries of all patches (official_only + osm_only)
2. Use `unary_union` to merge all touching ones into a single MultiPolygon
3. Explode the result into individual polygons
4. But we lose the source information

Alternative approach using networkx or just spatial operations:
1. For each patch, find which other patches it touches
2. Group connected components
3. For each component, merge geometries and determine source

Let me write a Python script to do this.

Actually, the simplest approach might be:

1. Compute all patches with source
2. Use geopandas dissolve with a connectivity groupby:
   - Create a graph where patches are nodes and edges exist if patches touch
   - Find connected components
   - Dissolve by component id

Let me implement this step by step.

Actually, we can use `shapely.ops.unary_union` on all patches - this will merge all touching/overlapping geometries. Then we need to map back to source. But unary_union loses the attribute info.

A better approach: use the `union_all` in newer geopandas, or iterate...

Let me think about another approach:

1. For each park, compute official_only, osm_only, and both
2. Only keep official_only and osm_only (the disagreement patches)
3. Merge the entire collection of disagreement patches by connectivity
4. For each merged group, check what source(s) it contains

To find connectivity-based groups, I can:
- Compute a spatial join of patches with themselves using a touch/intersects predicate
- Build a graph, find connected components
- Dissolve patches by component ID

Let me code this up.

Let me understand the problem:

1. We have two layers: `parks_official` and `parks_osm`, each with MultiPolygon geometries.
2. We need to find areas that are "claimed by exactly one side" — the symmetric difference between the two layers.
3. But with a twist: "adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two."
   - This means: for a park that exists in both layers but is shifted, there will be areas where official says park but OSM doesn't (loss), and areas where OSM says park but official doesn't (gain). These are adjacent and should be merged into one feature.
4. Each merged group gets a `source` attribute: `parks_official`, `parks_osm`, or `both`.
   - Wait, but the instruction says "every patch claimed by exactly one side". So we only want patches that are in exactly one of the two layers, not both. But then it says with `source` being `parks_official`, `parks_osm`, or `both`. That seems contradictory.
   
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."

Hmm, "exactly one side" means the symmetric difference. But then why would `source` be `both`? Let me re-read.

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

"First, parks_disagreement.geojson with the source column. Second, park_label_anchors.geojson with the same source column"

Wait, maybe the instruction means we need to categorize:
- Areas only in official (not in OSM): source = 'parks_official'
- Areas only in OSM (not in official): source = 'parks_osm'
- Areas in both: source = 'both'

But the first sentence says "every patch claimed by exactly one side". So maybe it's just the symmetric difference parts (only official, only OSM) but with adjacent ones merged. And then `both` appears if after merging, a merged group contains both a "only official" and "only OSM" patch? That would be the "shifted park's gain and loss" case.

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

So the algorithm would be:
1. Compute the symmetric difference between the two layers.
2. For each park (by park_id), if it exists in both:
   - Compute the areas where only official has it (official - osm)
   - Compute the areas where only OSM has it (osm - official)
   - These are the "gain" and "loss" for a shifted park.
   - If they are adjacent (touching), merge them into one multipolygon.
3. For parks only in one layer, those areas are already just "one side."
4. For parks in both but with no shift... they'd be in both, so no disagreement patch.

Wait, but the instruction says "every patch claimed by exactly one side". So we only want areas that are in one layer but not the other. But "adjacent disagreement patches merged together" means if an official-only patch and an OSM-only patch are adjacent (from the same park being shifted), they get merged.

And the source attribute... Let me think about what source means for a merged group that combines both an official-only area and an OSM-only area. It could be tagged 'both'? Or maybe the task is simpler: we just need the symmetric difference of the two full layers, then merge adjacent pieces, and assign each piece to whichever side(s) it came from.

Actually, I think the most natural interpretation is:
1. Compute areas only in official (A_only = official - osm)
2. Compute areas only in OSM (B_only = osm - official)
3. These are "disagreement patches"
4. Merge adjacent disagreement patches together (i.e., if an A_only patch touches a B_only patch, merge them)
5. For the merged group:
   - If it came only from A_only areas: source = 'parks_official'
   - If it came only from B_only areas: source = 'parks_osm'
   - If it came from both (merging adjacent A_only and B_only patches): source = 'both'

OR maybe the simpler interpretation: just compute the symmetric difference (areas in exactly one of the two layers), dissolve/merge adjacent pieces, and assign source based on whether the original piece was from official or osm. But if a merged piece contains both official-only and osm-only original pieces, then source = 'both'.

Let me just implement this approach:

1. Compute `official_only = official.geometry.unary_union - osm.geometry.unary_union`
2. Compute `osm_only = osm.geometry.unary_union - official.geometry.unary_union`
3. These give us all the disagreement patches.
4. Collect all individual polygons from both sides.
5. Merge adjacent/touching polygons across both sides.
6. For each merged group, determine its source.

Actually, using unary_union on the whole layer might be too aggressive. Let me think about this differently.

Each park has a park_id. Parks that are in both layers might have shifted geometry. Parks that are only in one layer are solely from that source.

For parks in both with the same geometry: no disagreement => ignore.
For parks in both with different geometry: the symmetric difference gives us the disagreement areas (gain and loss for that park).
For parks only in official: the whole park is a disagreement area (only official).
For parks only in OSM: the whole park is a disagreement area (only OSM).

But the "adjacent disagreement patches merged together" suggests that for a single shifted park, the gain area (OSM-only part) and loss area (official-only part) should be merged if they're adjacent. And across different parks... probably not.

Actually, I think the instruction means: take ALL disagreement patches from the entire dataset, then merge any that are adjacent (touching), regardless of which park they originally belonged to. This way, if a park is shifted, its gain and loss are adjacent and will be merged.

Let me think about the simplest, most robust approach:

Step 1: Compute symmetric difference between the full official layer and full OSM layer.
- `official_union = official.geometry.unary_union`
- `osm_union = osm.geometry.unary_union`
- `only_official = official_union - osm_union`
- `only_osm = osm_union - official_union`

But this loses per-park context. However, for the "merged groups" requirement, this might actually be the cleanest approach because:
- `only_official` and `only_osm` are both MultiPolygon/MultiLineString geometries
- We need to combine them, then merge adjacent parts

Let me try:

1. Get all polygons from `only_official` and `only_osm`
2. For each polygon, tag it with its source ('parks_official' or 'parks_osm')
3. Merge adjacent polygons (from any source) together
4. For each merged group, determine the source:
   - If all constituent polygons are from 'parks_official': source = 'parks_official'
   - If all from 'parks_osm': source = 'parks_osm'
   - If mixed: source = 'both'
5. Compute a label anchor point for each merged group (point strictly inside the geometry)

For merging adjacent polygons, I can use unary_union + dissolve by a group ID from spatial clustering (e.g., using DBSCAN or union cascade).

Let me implement this in Python using shapely and geopandas.

Actually, `unary_union` already dissolves overlapping/adjacent geometries into one. But we need to track the source. Let me think...

A simpler approach:
1. Create two GeoDataFrames: one with `only_official` polygons (source='parks_official'), one with `only_osm` polygons (source='parks_osm').
2. Combine them.
3. Use a spatial clustering approach to group adjacent polygons. I can use `dissolve` by a cluster ID obtained from spatial join of buffered geometries, or use the `unary_union` + polygonize approach.

Wait, let me think more carefully. `unary_union` on a set of polygons will merge them if they touch or overlap. But I need to preserve which source each part came from, and after merging, determine the source of the merged geometry.

Approach:
1. Get individual polygons from `only_official` and `only_osm`, each tagged with source.
2. To merge adjacent ones: I can use `shapely.ops.unary_union` on all of them together, which will produce a merged set of polygons. But then I lose the source information.
3. Alternative: use a graph-based approach where I find connected components of polygons that touch.

Let me use a simpler approach:
- For each polygon, buffer it by 0 (or use a very tiny epsilon) to create a "touch graph"
- Find connected components
- Merge polygons within each component
- Determine source from the original polygons

Actually, the simplest approach might be:
1. Compute `only_official` multipolygon and `only_osm` multipolygon from the full layers.
2. Get individual polygons from each.
3. To find which ones are adjacent, I can use a spatial index.

Let me code this up.

Actually wait - I think there might be a much simpler interpretation. 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 "every patch claimed by exactly one side" = symmetric difference areas (official\osm U osm\official).

"adjacent disagreement patches merged together" = if an official-only patch touches an osm-only patch, merge them.

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

So if a merged group contains only official-only areas -> 'parks_official'
If it contains only osm-only areas -> 'parks_osm'
If it contains both types (which happens when they're adjacent) -> 'both'

This makes perfect sense.

Let me implement:

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

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

# Compute the full unions
official_union = unary_union(official.geometry)
osm_union = unary_union(osm.geometry)

# Areas only in official
only_official = official_union - osm_union
# Areas only in osm
only_osm = osm_union - official_union

# Extract individual polygons
def extract_polygons(geom, source):
    if geom.is_empty:
        return []
    if geom.geom_type == 'Polygon':
        return [(geom, source)]
    elif geom.geom_type == 'MultiPolygon':
        return [(g, source) for g in geom.geoms]
    return []

polys = []
polys.extend(extract_polygons(only_official, 'parks_official'))
polys.extend(extract_polygons(only_osm, 'parks_osm'))
```

Then find adjacent polygons and merge them. For finding adjacency, I can use a spatial index and check if polygons touch or intersect.

For merging, I'll use a graph connected-components approach.

Then for each merged group, determine source and compute a label anchor point (representative point strictly inside).

Let me write the full script.

For the anchor point: I need a point that falls strictly inside the geometry. `geometry.representative_point()` gives a point guaranteed to be within the geometry. But it says "strictly inside" - I'll use `representative_point()` which works well for this.

Let me write it all out.

The output looks good:
- 46 disagreement features (MultiPolygon), 46 anchor points
- Source counts: 14 'both', 20 'parks_official', 12 'parks_osm'
- All anchors are strictly inside their corresponding polygons
- Source columns match between the two files

However, I notice the CRS changed from EPSG:6539 (the original GPKG) to EPSG:4326. GeoJSON is supposed to be in WGS84 (EPSG:4326) by default, so that's actually correct since GeoJSON specification requires coordinates in WGS84. But let me check if the coordinates were reprojected or if they stayed in the original CRS...

Actually, looking at the output, it says EPSG:4326. The original data was EPSG:6539 (NY Long Island state plane). Geopandas's to_json() auto-reprojects to WGS84 for GeoJSON output, which is correct per the GeoJSON spec. So this is fine.

Let me also verify that the source values match exactly the required strings.

Let me think about the logic more carefully.

The task says:
- "every patch claimed by exactly one side" - so areas in only official OR only OSM, not both
- "with adjacent disagreement patches merged together" - so touching/adjacent patches from the same or different sources should be merged
- "so that a shifted park's gain and loss form a single feature rather than two" - this means if a park shifted, the area only in official and the area only in OSM, if they're adjacent, should be merged into one feature with source "both"

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 "patches claimed by exactly one side" - this is the symmetric difference: (official ∪ osm) \ (official ∩ osm). That's the areas where only one side claims them.

"with adjacent disagreement patches merged together" - if these patches are adjacent, merge them.

"so that a shifted park's gain and loss form a single feature rather than two" - when a park shifts, the area gained (only in new) and area lost (only in old) are adjacent, so they merge into one feature. This merged feature would have both sources present, so source = 'both'.

This is exactly what I implemented. Let me verify the counts make sense:
- official-only parks: 20 (these are parks that only exist in official dataset)
- OSM-only parks: 12 (these only exist in OSM dataset)
- 15 parks have geometry differences (the symmetric diff of matching parks)

So we'd expect:
- 20 official-only park areas
- 12 OSM-only park areas
- Some number of "both" areas from the 15 parks with shifted geometries (some may be adjacent and merge)

Total disagreement patches initially: 20 + 12 + some number from the 15 shifted parks. After merging adjacent ones, we get 46 groups. 20 official-only, 12 OSM-only, 14 "both". That makes sense!

Let me just double-check that the source logic is right. For a component that has both official-only and osm-only polygons touching each other (a shifted park), the source should be 'both'. My code does:
- If all sources in the component are 'parks_official' -> 'parks_official'
- If all sources in the component are 'parks_osm' -> 'parks_osm'
- Otherwise -> 'both'

That's correct.

Let me also verify that the representative_point is strictly inside. The representative_point() method from Shapely returns a point guaranteed to be within the geometry. Let me verify this is the case.

Actually, looking at the Shapely docs: `representative_point()` returns a point that is guaranteed to be within the geometry. So it should be strictly inside.

The check passed: "Anchors NOT inside their polygon: 0" - all 46 anchors are inside their polygons.

Let me also verify the output formatting more carefully - GeoJSON should be valid.

Everything looks good. Let me summarize what was done:

1. Loaded both layers from the GeoPackage
2. Made all geometries valid
3. Computed union of official parks and union of OSM parks
4. Found symmetric difference (official-only, OSM-only)
5. Extracted individual polygons from those difference areas
6. Merged adjacent/touching polygons into groups using an adjacency graph (connected components)
7. For each merged group, determined the source: 'parks_official' if only official polygons, 'parks_osm' if only OSM polygons, 'both' if the group contains both types
8. Created MultiPolygon geometries for each group
9. Computed a representative point (strictly inside) for each group as the label anchor
10. Wrote both GeoJSON FeatureCollections

Results:
- 46 disagreement groups (MultiPolygon): 20 official-only, 12 OSM-only, 14 both
- 46 label anchors (Point), one per group in matching order
- All anchors are strictly inside their corresponding disagreement polygons
- All source values match exactly the required strings
- Both files are valid GeoJSON