run-20260616-001510Z / dc-l2-lagos-snap-normalize

0:00

I need to get our legacy zoning layer ready for the state portal. lagos_zoning_legacy.gpkg has inconsistent zoning_class values (mixed casing, spelling variations, and some blank entries) plus vertex precision issues that cause tiny gaps and overlaps between adjacent parcels. Can you clean up the class names to consistent canonical TitleCase values, drop any rows with blank or whitespace-only classes, fix the vertex precision problems, and produce one merged geometry per canonical zoning class? Please write the result to zoning_aggregated.gpkg, with each row carrying two columns: zoning_class and a numeric area_m2.

0:00

System info

{
  "adapter": "openrouter",
  "model": "deepseek/deepseek-v4-flash",
  "base_url": "https://openrouter.ai/api/v1"
}

0:03

Let me start by exploring the input file to understand its structure and contents.

0:03

Bash cd /work && python3 -c "

cd /work && python3 -c "
import geopandas as gpd
gdf = gpd.read_file('lagos_zoning_legacy.gpkg')
print('Columns:', list(gdf.columns))
print('Shape:', gdf.shape)
print(gdf.head(10).to_string())
print('---')
print('CRS:', gdf.crs)
print('---')
print('zoning_class values:')
print(gdf['zoning_class'].value_counts(dropna=False))
print('---')
print('Geometry types:', gdf.geometry.type.unique())
"

Columns: ['parcel_id', 'lga_source', 'zoning_class', 'area_m2', 'geometry']
Shape: (10080, 5)
   parcel_id lga_source zoning_class  area_m2                                                                               geometry
0          1   Surulere  RESIDENTIAL    100.0  POLYGON ((540000 720000, 540010 720000, 540010 720010, 540000 720010, 540000 720000))
1          2       Yaba  residential    100.0  POLYGON ((540010 720000, 540020 720000, 540020 720010, 540010 720010, 540010 720000))
2          3      Apapa        Resi.    100.0  POLYGON ((540020 720000, 540030 720000, 540030 720010, 540020 720010, 540020 720000))
3          4    Eti-Osa        resi.    100.0  POLYGON ((540030 720000, 540040 720000, 540040 720010, 540030 720010, 540030 720000))
4          5        Ojo        RESI.    100.0  POLYGON ((540040 720000, 540050 720000, 540050 720010, 540040 720010, 540040 720000))
5          6      Ikeja  Residential    100.0  POLYGON ((540050 720000, 540060 720000, 540060 720010, 540050 720010, 540050 720000))
6          7   Surulere  RESIDENTIAL    100.0  POLYGON ((540060 720000, 540070 720000, 540070 720010, 540060 720010, 540060 720000))
7          8       Yaba  residential    100.0  POLYGON ((540070 720000, 540080 720000, 540080 720010, 540070 720010, 540070 720000))
8          9      Apapa        Resi.    100.0  POLYGON ((540080 720000, 540090 720000, 540090 720010, 540080 720010, 540080 720000))
9         10    Eti-Osa        resi.    100.0  POLYGON ((540090 720000, 540100 720000, 540100 720010, 540090 720010, 540090 720000))
---
CRS: EPSG:26331
---
zoning_class values:
zoning_class
Residential     420
COMMERCIAL      420
industrial      420
Agri.           420
residential     419
resi.           419
Comm.           419
COMM.           419
indus.          419
AGRI.           419
Industrial      418
AGRICULTURAL    418
RESIDENTIAL     417
RESI.           417
comm.           417
commercial      417
Indus.          417
INDUSTRIAL      417
Agricultural    417
agri.           417
Resi.           416
Commercial      416
INDUS.          416
agricultural    416
                 13
                 13
NaN              12
\t               12
Name: count, dtype: int64
---
Geometry types: <ArrowStringArray>
['Polygon']
Length: 1, dtype: str

0:09

Let me analyze the data:

1. **zoning_class values** - There are many inconsistent values:
   - Residential variants: 'RESIDENTIAL', 'residential', 'Residential', 'RESI.', 'resi.', 'Resi.'
   - Commercial variants: 'COMMERCIAL', 'commercial', 'Commercial', 'Comm.', 'COMM.', 'comm.'
   - Industrial variants: 'industrial', 'Industrial', 'INDUSTRIAL', 'Indus.', 'indus.', 'INDUS.'
   - Agricultural variants: 'Agri.', 'AGRI.', 'AGRICULTURAL', 'Agricultural', 'agricultural', 'agri.'
   - Blank/whitespace: '' (13), '\t' (12), NaN (12)

2. Let me check the geometry to understand the vertex precision issues.

Let me look at the geometries more carefully.

0:09

Bash cd /work && python3 -c "

cd /work && python3 -c "
import geopandas as gpd
gdf = gpd.read_file('lagos_zoning_legacy.gpkg')
# Check some geometries for vertex precision
print('Sample geometries:')
for i in range(5):
    print(f'{i}: {gdf.geometry.iloc[i]}')
print('...')
# Check bounds
print('Total bounds:', gdf.total_bounds)
# Check if there are gaps by looking at unique coordinate values
import numpy as np
coords = np.array([list(g.exterior.coords) for g in gdf.geometry])
print('Coord sample:', coords[:3])
"

Sample geometries:
0: POLYGON ((540000 720000, 540010.00002 720000.000005, 540010.000005 720010.000015, 540000.000025 720010.00002, 540000 720000))
1: POLYGON ((540010.00001 720000.000025, 540020 720000.000005, 540020.00002 720010.00001, 540010.000005 720010.00002, 540010.00001 720000.000025))
2: POLYGON ((540020.000025 720000.00002, 540030.000015 720000, 540030 720010.00001, 540020.00002 720010.000015, 540020.000025 720000.00002))
3: POLYGON ((540030.000005 720000.00002, 540040.000025 720000.000025, 540040.000015 720010.000005, 540030 720010.000015, 540030.000005 720000.00002))
4: POLYGON ((540040.00002 720000.000015, 540050.000005 720000.000025, 540050.00003 720010, 540040.000015 720010.00001, 540040.00002 720000.000015))
...
Total bounds: [540000.       720000.       542547.       721000.000025]
Coord sample: [[[540000.       720000.      ]
  [540010.00002  720000.000005]
  [540010.000005 720010.000015]
  [540000.000025 720010.00002 ]
  [540000.       720000.      ]]

 [[540010.00001  720000.000025]
  [540020.       720000.000005]
  [540020.00002  720010.00001 ]
  [540010.000005 720010.00002 ]
  [540010.00001  720000.000025]]

 [[540020.000025 720000.00002 ]
  [540030.000015 720000.      ]
  [540030.       720010.00001 ]
  [540020.00002  720010.000015]
  [540020.000025 720000.00002 ]]]

0:23

Let me understand the data:

1. **zoning_class values** - There are many inconsistencies:
   - Full word forms: residential, RESIDENTIAL, Residential, resi., Resi., RESI.
   - Commercial: commercial, COMMERCIAL, Commercial, comm., Comm., COMM.
   - Industrial: industrial, INDUSTRIAL, Industrial, indus., Indus., INDUS.
   - Agricultural: agricultural, AGRICULTURAL, Agricultural, agri., Agri., AGRI.
   
   The canonical TitleCase forms should be: "Residential", "Commercial", "Industrial", "Agricultural"

2. **Blank/whitespace entries** - There are 13 blank, 13 with whitespace, 12 NaN, and 12 with tabs. Need to drop those.

3. **Vertex precision** - Coordinates have tiny variations like 540000 vs 540000.00002, 720000 vs 720000.000005. Need to snap/round to clean them up.

4. **Final output** - One merged geometry per canonical zoning class, with zoning_class and area_m2 columns.

Let me write a comprehensive Python script.

For the vertex precision: I'll round coordinates to a reasonable precision (e.g., 0.001 or maybe just to integers since these seem to be 10m grid based). Looking at the data, the base coordinates are 540000, 540010, 540020... with small perturbations. I should snap to a grid. Since they seem to be on a 10-meter grid, I'll round to 0 decimal places (nearest integer) or maybe to 3 decimal places. But actually the perturbations like 0.00002 suggest these are sub-mm. Let me round to nearest 0.001 or maybe just to integers since the pattern is clearly a 10m grid.

Actually, let me think about this more carefully. The coordinates are like:
540000, 540010, 540020... (10m spacing)
with perturbations like +/- 0.00002

So rounding to the nearest integer would clean them up perfectly. But let me check if there are any coordinates that aren't aligned to 10m.

Actually, looking at the bounds: 540000 to 542547, and 720000 to 721000.000025. The y range is roughly 1000 units. Let me just round to the nearest 1 meter (0 decimal places) to fix vertex precision.

Then after cleaning vertices, I'll:
1. Standardize zoning_class names
2. Drop blank/whitespace rows
3. Dissolve by zoning_class
4. Calculate area_m2
5. Write to zoning_aggregated.gpkg

Let me write the script.

0:23

Bash cd /work && python3 -c "

cd /work && python3 -c "
import geopandas as gpd
gdf = gpd.read_file('lagos_zoning_legacy.gpkg')
# Check all unique coordinate deltas to understand spacing
import numpy as np
all_coords = []
for g in gdf.geometry:
    all_coords.extend(list(g.exterior.coords))
all_coords = np.array(all_coords)
print('X min/max:', all_coords[:,0].min(), all_coords[:,0].max())
print('Y min/max:', all_coords[:,1].min(), all_coords[:,1].max())
print('Unique X values (first 20):', sorted(np.unique(all_coords[:,0].round(6)))[:20])
print('Unique X values (last 20):', sorted(np.unique(all_coords[:,0].round(6)))[-20:])
# Check if values are mostly on integer grid
x_mod = np.mod(all_coords[:,0], 1)
print('X fractional parts range:', x_mod.min(), x_mod.max())
print('X non-integer count:', np.sum(x_mod > 0.001))
y_mod = np.mod(all_coords[:,1], 1)
print('Y non-integer count:', np.sum(y_mod > 0.001))
"

X min/max: 540000.0 542547.0
Y min/max: 720000.0 721000.000025
Unique X values (first 20): [np.float64(540000.0), np.float64(540000.000005), np.float64(540000.00001), np.float64(540000.000015), np.float64(540000.00002), np.float64(540000.000025), np.float64(540000.00003), np.float64(540010.0), np.float64(540010.000005), np.float64(540010.00001), np.float64(540010.000015), np.float64(540010.00002), np.float64(540010.000025), np.float64(540010.00003), np.float64(540020.0), np.float64(540020.000005), np.float64(540020.00001), np.float64(540020.000015), np.float64(540020.00002), np.float64(540020.000025)]
Unique X values (last 20): [np.float64(542525.0), np.float64(542525.5), np.float64(542526.0), np.float64(542527.0), np.float64(542530.0), np.float64(542530.5), np.float64(542531.0), np.float64(542532.0), np.float64(542535.0), np.float64(542535.5), np.float64(542536.0), np.float64(542537.0), np.float64(542540.0), np.float64(542540.5), np.float64(542541.0), np.float64(542542.0), np.float64(542545.0), np.float64(542545.5), np.float64(542546.0), np.float64(542547.0)]
X fractional parts range: 0.0 0.5
X non-integer count: 30
Y non-integer count: 0

0:41

Let me analyze the data:

1. **zoning_class values**: There are multiple variations:
   - Residential: `RESIDENTIAL`, `Residential`, `residential`, `Resi.`, `resi.`, `RESI.`, `RESI.`
   - Commercial: `COMMERCIAL`, `Commercial`, `commercial`, `Comm.`, `comm.`, `COMM.`, `COMM.`
   - Industrial: `INDUSTRIAL`, `Industrial`, `industrial`, `Indus.`, `indus.`, `INDUS.`, `INDUS.`
   - Agricultural: `AGRICULTURAL`, `Agricultural`, `agricultural`, `Agri.`, `agri.`, `AGRI.`, `AGRI.`
   - Also blank/whitespace: `'' (13)`, `NaN (12)`, `\t (12)`

2. **Vertex precision issues**: Coordinates have tiny offsets like 0.000005, 0.00001, etc. These should be snapped to the nearest integer (grid).

Let me check more carefully what the canonical classes should be. Looking at the data:
- Residential, RESIDENTIAL, residential, Resi., resi., RESI. → "Residential"
- COMMERCIAL, Commercial, commercial, Comm., comm., COMM. → "Commercial"
- INDUSTRIAL, Industrial, industrial, Indus., indus., INDUS. → "Industrial"
- AGRICULTURAL, Agricultural, agricultural, Agri., agri., AGRI. → "Agricultural"

Now I need to:
1. Clean zoning_class to canonical TitleCase
2. Drop rows with blank/whitespace-only classes
3. Fix vertex precision (snap to integer grid since the grid seems to be based on 10m intervals)
4. Union/dissolve by zoning_class
5. Compute area_m2
6. Write to zoning_aggregated.gpkg

Let me write the script.

0:41

Write /work/solve.py

content · 79 lines · py

import geopandas as gpd
import numpy as np
from shapely import wkt
from shapely.geometry import Polygon
import re

# ── 1. Read ──────────────────────────────────────────────────────────────
gdf = gpd.read_file("/work/lagos_zoning_legacy.gpkg")
print(f"Read {len(gdf)} rows")

# ── 2. Drop rows with blank / whitespace-only / NaN zoning_class ─────────
gdf = gdf[gdf["zoning_class"].notna()].copy()
gdf = gdf[gdf["zoning_class"].str.strip().astype(bool)].copy()
print(f"After dropping blank/NaN: {len(gdf)} rows")

# ── 3. Normalise zoning_class to canonical TitleCase ────────────────────
# Build a mapping from the raw values to canonical names
def canonicalise(val):
    v = str(val).strip().lower()
    # Handle abbreviated versions
    if v in ("residential", "resi."):
        return "Residential"
    if v in ("commercial", "comm."):
        return "Commercial"
    if v in ("industrial", "indus."):
        return "Industrial"
    if v in ("agricultural", "agri."):
        return "Agricultural"
    # Full uppercase variants (already covered by lower() above)
    return val.strip().title()

gdf["zoning_class"] = gdf["zoning_class"].apply(canonicalise)

# Double-check all values
print("Canonicalised zoning_class distribution:")
print(gdf["zoning_class"].value_counts())

# ── 4. Snap vertex coordinates to the nearest integer metre ──────────────
# The grid looks like it should be integer metres (10m spacing)
# We'll round all coordinates to the nearest integer (0 decimal places)

def snap_polygon(geom):
    """Snap all exterior & interior ring coords to the nearest integer."""
    if geom is None or geom.is_empty:
        return geom

    def snap_ring(coords):
        return [(round(x), round(y)) for x, y in coords]

    exter = snap_ring(geom.exterior.coords)
    inters = [snap_ring(r.coords) for r in geom.interiors]
    return Polygon(exter, inters)

gdf["geometry"] = gdf["geometry"].apply(snap_polygon)

# Remove any invalid or empty geometries that might result
gdf = gdf[gdf.is_valid & ~gdf.is_empty].copy()
print(f"After snapping & validity check: {len(gdf)} rows")

# ── 5. Dissolve / union by zoning_class ─────────────────────────────────
# Since the task says "one merged geometry per canonical zoning class",
# we need to union all parcels of the same class into a single geometry.

dissolved = gdf.dissolve(by="zoning_class", aggfunc=None)
dissolved = dissolved.reset_index()

# Compute area in the layer's CRS (EPSG:26331 — metres, so area is in m²)
dissolved["area_m2"] = dissolved.geometry.area

# Keep only the two requested columns
result = dissolved[["zoning_class", "area_m2", "geometry"]].copy()
result = result.sort_values("zoning_class").reset_index(drop=True)

print("\nFinal result:")
print(result[["zoning_class", "area_m2"]].to_string())

# ── 6. Write output ─────────────────────────────────────────────────────
result.to_file("/work/zoning_aggregated.gpkg", layer="zoning_aggregated", driver="GPKG")
print("\n✅ Written to /work/zoning_aggregated.gpkg")

Wrote /work/solve.py (3006 chars).

0:43

Let me check how many unique canonical classes we'd get from the mapping, and also look more carefully at the data to make sure my mapping covers everything.

0:43

Bash cd /work && python3 -c "

cd /work && python3 -c "
import geopandas as gpd
gdf = gpd.read_file('lagos_zoning_legacy.gpkg')
# Check ALL unique zoning_class values
vals = gdf['zoning_class'].unique()
for v in sorted(vals, key=str):
    print(repr(v))
"

''
'\t'
'   '
'AGRI.'
'AGRICULTURAL'
'Agri.'
'Agricultural'
'COMM.'
'COMMERCIAL'
'Comm.'
'Commercial'
'INDUS.'
'INDUSTRIAL'
'Indus.'
'Industrial'
'RESI.'
'RESIDENTIAL'
'Resi.'
'Residential'
'agri.'
'agricultural'
'comm.'
'commercial'
'indus.'
'industrial'
nan
'resi.'
'residential'

0:55

Let me review the mapping more carefully. The canonicalise function should handle all these variations. Let me check:

- '' (empty), '   ' (whitespace), '\t' (tab), nan → should be dropped (done in step 2)
- 'RESIDENTIAL', 'residential', 'Residential', 'RESI.', 'resi.', 'Resi.' → "Residential"
- 'COMMERCIAL', 'commercial', 'Commercial', 'COMM.', 'comm.', 'Comm.' → "Commercial"
- 'INDUSTRIAL', 'industrial', 'Industrial', 'INDUS.', 'indus.', 'Indus.' → "Industrial"
- 'AGRICULTURAL', 'agricultural', 'Agricultural', 'AGRI.', 'agri.', 'Agri.' → "Agricultural"

My current canonicalise function does:
1. v = str(val).strip().lower()
2. Then checks exact matches for lower-cased values.

Let me trace through:
- 'residential'.lower() = 'residential' → matches 'residential' → Returns "Residential" ✓
- 'RESIDENTIAL'.lower() = 'residential' → matches 'residential' → Returns "Residential" ✓
- 'Residential'.lower() = 'residential' → matches 'residential' → Returns "Residential" ✓
- 'RESI.'.lower() = 'resi.' → matches 'resi.' → Returns "Residential" ✓
- 'resi.'.lower() = 'resi.' → matches 'resi.' → Returns "Residential" ✓
- 'Resi.'.lower() = 'resi.' → matches 'resi.' → Returns "Residential" ✓

Good. But what about the '.title()' fallback at the end? That would turn 'resi.' into 'Resi.' which isn't our canonical. But since 'resi.' is matched earlier, it won't reach there.

Wait, what about 'RESI.'? lower() = 'resi.' → matches 'resi.' → "Residential" ✓

Good. But I want to make sure I handle all 4 classes. Let me also think about the fallback: `return val.strip().title()` - this would turn 'RESIDENTIAL' (if not caught) to 'Residential'. But since it IS caught by the lower check, it's fine. But just to be safe, maybe I should have a more robust approach.

Actually, let me also think: what if a value like '  Residential  ' came in? str(val).strip().lower() = 'residential' → matches. Good.

Let me also double check: the empty strings '', '   ', '\t' - after stripping, those become '' which is falsy, and since they were already dropped in step 2, they won't reach canonicalise. And NaN: pd.notna check drops them. So step 2 handles that.

Let me run the script and see.

0:55

Bash cd /work && python solve.py

Read 10080 rows
After dropping blank/NaN: 10030 rows
Canonicalised zoning_class distribution:
zoning_class
Residential     2508
Commercial      2508
Industrial      2507
Agricultural    2507
Name: count, dtype: int64
After snapping & validity check: 10000 rows

Traceback (most recent call last):
  File "/work/solve.py", line 64, in <module>
    dissolved = gdf.dissolve(by="zoning_class", aggfunc=None)
  File "/usr/local/lib/python3.14/site-packages/geopandas/geodataframe.py", line 2274, in dissolve
    aggregated_data = data.groupby(**groupby_kwargs).agg(aggfunc, **kwargs)
  File "/usr/local/lib/python3.14/site-packages/pandas/core/groupby/generic.py", line 2281, in aggregate
    relabeling, func, columns, order = reconstruct_func(func, **kwargs)
                                       ~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^
  File "/usr/local/lib/python3.14/site-packages/pandas/core/apply.py", line 1770, in reconstruct_func
    raise TypeError("Must provide 'func' or tuples of '(column, aggfunc).")
TypeError: Must provide 'func' or tuples of '(column, aggfunc).