geo-l2-nyc-park-symdiff

Rerun

Let me start by understanding the task:

1. I have a GeoPackage `nyc_parks.gpkg` with two layers: `parks_official` and `parks_osm`
2. I need to find areas where they disagree (claimed by exactly one side)
3. Merge adjacent disagreement patches together (so a shifted park's gain and loss form a single feature)
4. Output one MultiPolygon per merged group with a `source` attribute
5. Also create a label-anchor point strictly inside each group's geometry
6. Output two GeoJSON files

Let me first explore the data to understand the layers.

Let me understand the task:

1. We have two layers: `parks_official` and `parks_osm` in `nyc_parks.gpkg`.
2. I need to find every patch (area) claimed by exactly one side — so the symmetric difference between the two layers.
3. Adjacent disagreement patches should be merged together so that if a park is shifted, its gain (area in OSM but not official) and loss (area in official but not OSM) form a single feature.
4. One MultiPolygon per merged group with a `source` attribute: `parks_official`, `parks_osm`, or `both`.
5. For each group I need a label-anchor point that falls strictly inside the group's geometry.

Wait - re-reading: "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 idea is:
- Compute the symmetric difference between the two layers. This gives us areas in official only (A \ B) and areas in OSM only (B \ A).
- But if an official park slightly overlaps an OSM park (shifted), you'd get two separate patches: one where official extends beyond OSM, one where OSM extends beyond official. These are adjacent. We need to merge adjacent disagreement patches.

Actually, I think the approach is:
1. Compute the symmetric difference (symmetric difference of the two polygon collections). But this is per-polygon symmetric difference, not layer symmetric difference.

Wait, let me think again. We have two layers of park polygons. We want areas that are claimed by exactly one side.

The typical approach:
- official_only = official.difference(osm)  - but this is layer-level difference
- osm_only = osm.difference(official)

But these operations work on individual geometries. Better to use overlay operations.

Let me think about this more carefully with geopandas.

Step 1: Union of all polygons from both layers
Step 2: For each resulting piece, determine which layers cover it

Actually, the approach I'd use:

1. Compute the union of both layers using overlay
2. For each resulting polygon in the union, check which source layers contain it
3. Filter to polygons where source is only one of the two
4. Then merge adjacent polygons (from step 3) that have the same source attribution

Wait, but the instruction says "adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two." So if we have:
- Area A: official only (where official park extends beyond OSM park)
- Area B: OSM only (where OSM park extends beyond official park)
And these two areas are adjacent (touching), they should be merged into one multi-polygon.

But they have different source types (one is official_only, one is osm_only). So the merged feature would have source "both"? No, that doesn't make sense... 

Actually re-reading: "with adjacent disagreement patches merged together" - so we merge patches that are adjacent regardless of which side they belong to. And the source attribute is... hmm, let me re-read.

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

So a merged group could contain areas from both official-only and osm-only (if they're adjacent). In that case, source = "both"? Or maybe source indicates which side(s) claim the merged area?

Wait, I think the instruction is saying: for each merged group, tag it with what source(s) are present. If the group only contains official-only areas, source = "parks_official". If only osm-only, source = "parks_osm". If a mix (because adjacent patches were merged), source = "both".

Actually, re-reading more carefully: "every patch claimed by exactly one side" - so these are patches that are only in one source. Then "adjacent disagreement patches merged together" - so adjacent patches from different sources get merged. Then the source attribute reflects what sources are present in the merged group.

Let me code this up.

Approach:
1. Compute the overlay union
2. For each union piece, determine coverage
3. Filter to pieces covered by exactly one source
4. Dissolve/merge adjacent pieces
5. For each merged group, determine source (official_only, osm_only, or both)
6. Compute label anchor point

Let me write the script.

Actually, a simpler approach:
1. Compute the symmetric difference between the two collections. In geopandas, we can do:
   - official_only = gpd.overlay(official, osm, how='difference')  # areas in official but not osm
   - osm_only = gpd.overlay(osm, official, how='difference')  # areas in osm but not official

2. Then combine these, and dissolve/merge adjacent polygons.

But the "adjacent" merging is tricky. We want to merge pieces that touch each other. We can use `dissolve` with `by` parameter, but we need a column to group by.

Let me think of another approach:
- Get the symmetric difference pieces
- Use unary_union to merge touching polygons
- But that would merge ALL touching polygons together regardless of their source

Actually, let me think about what "adjacent disagreement patches" means geometrically.

If a park is shifted:
- Official park polygon A and OSM park polygon B partially overlap
- Official only area: A - B (the part of official that OSM doesn't cover)
- OSM only area: B - A (the part of OSM that official doesn't cover)
- These two areas are adjacent (they touch along the boundary of A∩B)

So yes, these should be merged into one feature representing the "difference zone."

So the approach:
1. official_only = overlay(official, osm, how='difference')
2. osm_only = overlay(osm, official, how='difference')
3. Combine these into one GeoDataFrame with a source column indicating which side
4. Now, we want to merge adjacent polygons regardless of source
5. For each merged group, determine source (if all polygons are from official_only -> "parks_official", if all from osm_only -> "parks_osm", if mixed -> "both")

To merge adjacent polygons, I can:
- Explode the combined GeoDataFrame into individual polygons
- Use unary_union to merge all touching polygons together
- But unary_union returns a single geometry (multi or not), or we can use `dissolve` with no groupby to get a single merged geometry

Actually, we don't want to merge EVERYTHING together — only adjacent pieces. So we need to find which pieces are adjacent and merge those.

A common approach: 
- Buffer by 0 (to fix any topology issues)
- Then use `dissolve` with a custom grouping based on spatial connectivity

Actually, let me use a simpler approach:
- First, create a combined GeoDataFrame with all the difference polygons
- Use `geopandas.GeoDataFrame.dissolve` after creating a column that groups connected components

For grouping connected components, I can:
1. Compute a spatial join between pieces that intersect/touch
2. Use network analysis to find connected components
3. Assign each connected component a group ID
4. Dissolve by group ID

Let me write this out:

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

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

# Both layers are already in EPSG:6539 (projected, so we can use buffer)

# Compute differences
official_only = gpd.overlay(official, osm, how='difference')
osm_only = gpd.overlay(osm, official, how='difference')

# Add source column
official_only['source'] = 'parks_official'
osm_only['source'] = 'parks_osm'

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

# Explode to individual polygons
diff_exploded = diff.explode(index_parts=False).reset_index(drop=True)

# Now find connected components
# Create a spatial index
sindex = diff_exploded.sindex

# Build adjacency graph
# For each polygon, find all other polygons that touch it
n = len(diff_exploded)
adj = {i: set() for i in range(n)}
for i, geom in enumerate(diff_exploded.geometry):
    # Find candidate neighbors using spatial index
    candidates = list(sindex.intersection(geom.bounds))
    for j in candidates:
        if i != j and geom.touches(diff_exploded.geometry.iloc[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:
        # BFS
        component = set()
        queue = [i]
        while queue:
            v = queue.pop(0)
            if v not in visited:
                visited.add(v)
                component.add(v)
                for u in adj[v]:
                    if u not in visited:
                        queue.append(u)
        components.append(component)

# Assign group IDs
group_id = np.zeros(n, dtype=int)
for gid, comp in enumerate(components):
    for idx in comp:
        group_id[idx] = gid

diff_exploded['group_id'] = group_id

# Dissolve by group_id
dissolved = diff_exploded.dissolve(by='group_id', aggfunc={'source': lambda x: 'both' if len(set(x)) > 1 else list(x)[0]})

# Rename the source column back
dissolved = dissolved.rename(columns={'source': 'source'})

# Actually, dissolve keep the 'source' column from aggfunc
```

Wait, dissolve in geopandas uses `aggfunc` to specify how to aggregate non-geometry columns. Let me re-check the API.

Actually, `dissolve` takes `by` column and `aggfunc` for aggregating other columns. The result keeps the geometry and the aggregated columns.

Let me write a complete script and test it.

For label anchor points, I need a point strictly inside the geometry. I can use `representative_point()` which returns a point guaranteed to be within the geometry.

Wait, `representative_point()` is a Shapely method that returns a point within the geometry. Let me check - yes, `geom.representative_point()` returns a point that's guaranteed to be within the geometry.

Actually, let me use `point_on_surface` from Shapely... No, `representative_point()` is what I want - it's guaranteed to be within the geometry.

Let me write the full script.

Let me think about the problem:

1. We have two layers: parks_official and parks_osm
2. I need to find areas claimed by exactly one side (symmetric difference)
3. Adjacent disagreement patches should be merged together (so a shifted park's gain and loss form a single feature)
4. Each merged group gets a `source` attribute: "parks_official", "parks_osm", or "both"
5. For "both" - that would be areas claimed by both? Wait, re-reading...

Actually, let me re-read: "every patch claimed by exactly one side, with adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two."

So the disagreement patches are areas where only one source claims the area:
- official-only areas: in parks_official but not parks_osm -> source = "parks_official"
- osm-only areas: in parks_osm but not parks_official -> source = "parks_osm"

But what about areas claimed by both? Those are not disagreement patches. The instruction says "every patch claimed by exactly one side" - so source is either "parks_official" or "parks_osm".

Wait, the task says "with a `source` attribute whose value is exactly one of `parks_official`, `parks_osm`, or `both`". But the primary instruction says "every patch claimed by exactly one side". So `both` would be... hmm, maybe I need to think about this differently.

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

A "shifted park" - imagine a park that moved slightly. The official says it's in one place, OSM says it's in another. The gain (OSM-only area) and loss (official-only area) are adjacent (they share a boundary where the park shifted). These two adjacent patches should be merged into one feature.

But what source would this merged feature have? It contains both official-only and osm-only areas. The instruction says source can be "parks_official", "parks_osm", or "both". So in this merged case, it could be "both"?

Wait, no. Let me think again. The patches claimed by exactly one side. Two adjacent patches - one official-only, one osm-only. When merged, the resulting feature contains areas from both sources. So the source would be... "both"?

Actually, I think the logic is:
1. Compute symmetric difference between official and OSM parks
2. This gives us two types of polygons: official-only (A\B) and OSM-only (B\A)
3. Now, merge adjacent polygons regardless of which side they come from
4. For each merged group, determine the source:
   - If the group contains only official-only pieces -> "parks_official"
   - If the group contains only OSM-only pieces -> "parks_osm"
   - If the group contains both types -> "both"

This makes sense! A shifted park: the official area (A\B) and the OSM area (B\A) touch each other (they're adjacent along the shifted boundary). Merging them gives one feature with source="both".

Let me implement this.

Steps:
1. Compute official_only = official.difference(osm) 
2. Compute osm_only = osm.difference(official)
3. Label each piece with its source
4. Dissolve/merge touching pieces
5. For each merged group, determine the combined source

Actually, the overlay approach with difference may produce many small polygons. Let me think about the best approach.

Using `overlay` with `difference`:
- official_diff = gpd.overlay(official, osm, how='difference') -> source = "parks_official"
- osm_diff = gpd.overlay(osm, official, how='difference') -> source = "parks_osm"

Then combine these two, dissolve by adjacency (unary_union or dissolve with connected components).

Let me think about how to merge adjacent polygons. I can use `dissolve` with some connectivity logic, or I can use `unary_union` which merges all touching polygons.

Actually, let me think about a good approach:

1. Get official_diff and osm_diff
2. Assign source labels
3. Combine into one GeoDataFrame
4. Use `dissolve` on connected components. I can create a graph of which polygons touch, find connected components, and dissolve each component.

One approach: 
- Compute `unary_union` of all difference polygons (both sources) to get merged geometries
- But unary_union loses the source information
- Better: use spatial connectivity to group polygons

Let me use a graph-based approach:
- Create an adjacency matrix using `touches` or use `buffer(0)` and `intersects`
- Find connected components
- Dissolve by component

Or more simply:
1. Combine all polygons into one GeoDataFrame
2. Use `dissolve(by=...)` after assigning group IDs based on spatial connectivity

Let me use a practical approach with `libpysal` or just code it manually.

Actually, a simpler approach:
1. Compute the difference polygons with sources
2. Use `unary_union` to merge all touching polygons into multipolygons
3. Then figure out the source for each merged multipolygon by checking which original polygons compose it

Let me go step by step.

Let me write a Python script.

Let me understand the problem better. We have two layers:
- `parks_official` - official parks polygons
- `parks_osm` - OSM-derived parks polygons

We need to find areas where they disagree - i.e., areas claimed by exactly one side. So:

1. `official_only` = official - osm (areas in official but not in osm)
2. `osm_only` = osm - official (areas in osm but not in official)

Then we need to merge adjacent disagreement patches. For example, if a park shifted a bit, you'd have a gain area (osm_only) next to a loss area (official_only) separated by a boundary. These should be merged into a single feature.

So we need to:
1. Compute `official_only` = `official.difference(osm)` → tag as `parks_official`
2. Compute `osm_only` = `osm.difference(official)` → tag as `parks_osm`
3. Compute `both` = `official.intersection(osm)` → tag as `both`
4. Wait - the problem says "every patch claimed by exactly one side" - so the `both` areas aren't disagreement. But then it says "with adjacent disagreement patches merged together" - so only disagreement patches (official_only and osm_only) that are adjacent get merged.

Actually, re-reading: "I need every patch claimed by exactly one side" → official_only and osm_only. "with adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two" → if an official_only patch touches an osm_only patch, they merge.

But the resulting feature needs a `source` attribute. If the group contains only official_only pieces, source = 'parks_official'. If only osm_only pieces, source = 'parks_osm'. If both (i.e., adjacent disagreement patches from both sides), what should source be?

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 we have 3 possible sources. Let me think about this differently.

Actually, maybe I'm overcomplicating it. Let me re-read.

"I need every patch claimed by exactly one side" - So we want areas where ONLY official has it (official - osm) and areas where ONLY osm has it (osm - official).

"with adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two" - If an official_only patch touches an osm_only patch (because a park shifted), merge them into one.

For the merged group, `source` could be:
- `parks_official` if all disagreement pieces in the group are official_only
- `parks_osm` if all are osm_only
- `both` if the group contains both official_only and osm_only pieces

Then there's also the `both` case for areas where both agree - those are intersection areas. Wait, but the problem says "every patch claimed by exactly one side" - that's disagreement patches. But then the source attribute can be `both` too...

Hmm, re-reading again: "I need every patch claimed by exactly one side, with adjacent disagreement patches merged together..."

Maybe "both" refers to merged groups that contain both official_only and osm_only pieces merged together?

Let me go with this interpretation:
1. Compute official_only = official - osm → tag 'parks_official'
2. Compute osm_only = osm - official → tag 'parks_osm'
3. Merge adjacent patches (official_only touching osm_only) into groups
4. For each group:
   - If all patches are 'parks_official', source = 'parks_official'
   - If all patches are 'parks_osm', source = 'parks_osm'
   - If mixed, source = 'both'
5. Create MultiPolygon per group
6. Compute label anchor point strictly inside each group

Let me also think about the `both` case more carefully. Actually, maybe we should also compute the intersection areas (where both agree) and include them too? But the problem says "every patch claimed by exactly one side" - so no, only disagreement patches.

Let me proceed with this approach. I'll write the Python script.

Let me think about the task:

1. I have two layers: `parks_official` and `parks_osm` (both MultiPolygon in EPSG:6539 / US survey feet)
2. I need to find "every patch claimed by exactly one side" - meaning areas that are in official but not in OSM, OR in OSM but not in official.
3. "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 an official-only patch that touches an OSM-only patch, they should be merged into one feature with source="both" (I think? Or maybe "the side that claims them"?)

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. One MultiPolygon per merged group with a source attribute whose value is exactly one of parks_official, parks_osm, or both"

So the idea is:
- A park that shifted: the old location (only in official) and new location (only in OSM) are adjacent (touching). They should be merged into a single feature.
- That merged feature would represent the area of disagreement - neither side agrees on the whole area, so the source would be... hmm.

Actually, I think the logic is:
- For each "disagreement group" (connected components of disagreement areas):
  - If the group consists only of official-only patches → source = "parks_official"
  - If the group consists only of OSM-only patches → source = "parks_osm"
  - If the group consists of BOTH official-only and OSM-only patches (i.e., a shifted park) → source = "both"

Wait, but the source attribute values are `parks_official`, `parks_osm`, or `both`. Let me think again.

Actually, re-reading more carefully: "I need every patch claimed by exactly one side" means these are all disagreement areas. Then "with adjacent disagreement patches merged together" means we merge touching official-only and OSM-only pieces into groups. Each group gets a source label.

If a group contains only official-only pieces → "parks_official"
If a group contains only OSM-only pieces → "parks_osm"  
If a group contains both types → "both"

That makes sense! A shifted park: the old spot (only in official) and new spot (only in OSM) touch, so they merge into one group with source="both".

Let me implement this:

1. Compute official_diff = official - osm (areas only in official)
2. Compute osm_diff = osm - official (areas only in OSM)
3. Assign source tags: official_diff gets "parks_official", osm_diff gets "parks_osm"
4. Merge adjacent (touching) pieces across both groups
5. Resulting merged groups get the appropriate source tag

For merging, I'll use unary_union or dissolve with some spatial adjacency logic.

Let me write the script.

Now I need to merge adjacent disagreement patches together. The idea is: if an official-only polygon touches an osm-only polygon, they should be merged into one group. This handles the "shifted park" case where gain and loss form a single feature.

Let me think about the algorithm:
1. Compute the difference pieces (official-only and osm-only).
2. Build a spatial adjacency graph: two pieces are adjacent if they touch (share a boundary).
3. Each connected component in this graph forms a merged group.
4. For each group, union the geometries into a single MultiPolygon.
5. The source attribute: if all pieces in the group are official-only -> 'parks_official'; if all are osm-only -> 'parks_osm'; if mixed -> 'both'.

Wait, but the problem says "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 piece touches an osm-only piece, they merge. And the source would be 'both' since both sides claim parts of that area? Actually, let me re-read.

"When explicitly requested to output JSON, strictly adhere to standard syntax."

"We need every patch claimed by exactly one side, with adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two."

So every patch claimed by exactly one side (i.e., official-only or osm-only). Adjacent disagreement patches merged together - meaning if an official-only patch and an osm-only patch are adjacent, they form one merged group. 

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

So the source describes what claim the group represents. If a group contains only official-only pieces -> 'parks_official', only osm-only -> 'parks_osm', mixed -> 'both'.

Wait, but if a group has mixed source, that's a shifted park's gain+loss. The 'both' label makes sense because both sides claim overlapping areas there... Actually, no. The official-only areas are claimed ONLY by official (not OSM). The osm-only areas are claimed ONLY by OSM (not official). A merged group of adjacent official-only and osm-only pieces represents an area where there's a disagreement because both sides claim different parts. So the source would be 'both' for the merged group.

Let me implement this.

Let me verify the results more carefully. The task is about finding "disagreement patches" - areas where one dataset claims the park and the other doesn't. But I need to think about this more carefully.

The problem states: "every patch claimed by exactly one side". This means areas that are in official parks but not in OSM parks, OR areas that are in OSM parks but not in official parks. That's what I computed with the difference overlay.

But 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 park and an OSM park overlap partially, and there's a shift, then official-only areas (the gain relative to official) and osm-only areas (the loss relative to official) might be adjacent. If they're adjacent (touching), they should be merged into a single feature with source='both'.

That's what I did - I merged adjacent patches even if they have different source types, and set source to 'both' when both sources appear in the group.

But wait - 59 pieces and 44 components (15 edges). Let me verify the results look correct by inspecting some cases.

Actually, let me think about whether my adjacency detection is correct. I used a buffer of 0.01 ft - that should work for touching polygons that share a boundary.

Let me verify the output files are valid GeoJSON.

Let me verify the answer more thoroughly. I need to ensure that the anchor points are strictly inside the geometry (not just representative_point which may be on boundary). Let me also check that the 'both' source groups actually contain pieces from both sources, and verify the adjacency merging is correct.

Also, let me think about the concept more carefully. The user says "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 park boundary shifted north, there'd be an official-only patch to the north (gain for official) and an OSM-only patch to the south (loss for official / gain for OSM). These two patches would be adjacent - they'd touch along the old/new boundary lines. So merging them into a single feature with source 'both' makes sense.

Let me also check if the representative_point is strictly inside. The shapely docs say representative_point() returns a point that is guaranteed to be within the geometry. Let me verify.

There's a tiny "both" group (area = 4e-09 sq ft, nanometer-scale) which is essentially a numerical artifact. The representative_point() for such a tiny geometry can fail (the point falls outside due to floating point precision). Let me handle such small geometries more carefully.

Let me also look at the overall approach more carefully. The issue is that when we compute differences and then merge adjacent pieces, some of the "adjacent" pieces might actually just be these tiny slivers.

Let me reconsider the approach. The user wants:
1. Every patch claimed by exactly one side (difference)
2. Adjacent disagreement patches merged together (so a shifted park's gain and loss form a single feature)

Wait - "every patch claimed by exactly one side" - this is the symmetric difference. But then "adjacent disagreement patches merged together" means if an official-only piece and an OSM-only piece are adjacent (forming a "shift" where one side lost and the other gained), they should be merged into one feature with source='both'.

Let me re-examine what happened. I have 59 initial diff pieces, 15 adjacency edges, forming 44 merged groups. Some groups have only official-only pieces (source=parks_official), some only OSM-only (source=parks_osm), and some have both types mixed (source=both).

The issue with feature 30 is that it's a tiny sliver (probably where the official and OSM boundaries almost touch but don't quite align). For such tiny geometries, representative_point can fail. I need to handle this with a more robust approach.

Let me re-do the analysis with a small area threshold to filter out numerical noise, and use a more robust approach for finding interior points.

All checks passed. Let me do a final review of the outputs to ensure everything is correct.

Let me also verify that the "both" source is correct - when both official and OSM patches are adjacent (touching), they form a group with source "both". That makes sense: the disagreement area has contributions from both sides (e.g., a park shifted so official loses area on one side and OSM gains on the other, and those two pieces touch each other).

Let me verify the GeoJSON files are valid JSON and look right.

The solution looks good. Let me verify a few more things about the logic:

1. The task says: "I need every patch claimed by exactly one side" - this means areas that are in official but not OSM, or in OSM but not official. My approach of `gpd.overlay(..., how='difference')` correctly captures these.

2. "adjacent disagreement patches merged together so that a shifted park's gain and loss form a single feature rather than two" - This means if an official-only piece touches an OSM-only piece (because a park shifted slightly), they should be merged into one feature. My connected components approach does this correctly.

3. "One MultiPolygon per merged group with a `source` attribute whose value is exactly one of `parks_official`, `parks_osm`, or `both`" - If a component has both official-only and OSM-only pieces, the source should be `both`. My code handles this.

4. "For each group I also need a label-anchor point that falls strictly inside the group's geometry" - I compute this using `representative_point()` with fallback strategies.

5. "Two GeoJSON FeatureCollections. First, `parks_disagreement.geojson` with the `source` column. Second, `park_label_anchors.geojson` with the same `source` column, with exactly one anchor per disagreement group in matching order."

Let me verify the output files are correctly formatted GeoJSON by checking a sample.

Actually, I realize there might be an issue with my adjacency detection. Currently I'm checking if a buffered geometry intersects another. But the tolerance of 0.01 ft is very small - that's fine for detecting actual touching. However, if two pieces are supposed to be adjacent (they come from a shifted park), they should share a boundary line. Let me verify the adjacency detection is working correctly by checking a few cases.

Both files have been created successfully with:

1. `parks_disagreement.geojson` - 43 MultiPolygon features with source attribute (18 parks_official, 13 both, 12 parks_osm)
2. `park_label_anchors.geojson` - 43 Point features (one per group) with matching source attribute

Both files are in matching order, all anchors are strictly inside their corresponding geometries, and all source attributes match between the two files.