dc-l2-lagos-snap-normalize
Map (reference ← swipe → agent)
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:04
Let me start by exploring the input file to understand its structure and data.
0:04
Bash
ls -la /work/lagos_zoning_legacy.gpkg
-rw-r--r-- 1 agent agent 2527232 Jun 18 23:05 /work/lagos_zoning_legacy.gpkg
0:04
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import geopandas as gpd
gdf = gpd.read_file('lagos_zoning_legacy.gpkg')
print('Columns:', gdf.columns.tolist())
print('CRS:', gdf.crs)
print('Shape:', gdf.shape)
print('dtypes:')
print(gdf.dtypes)
print()
print('First 5 rows:')
print(gdf.head())
"
Columns: ['parcel_id', 'lga_source', 'zoning_class', 'area_m2', 'geometry'] CRS: EPSG:26331 Shape: (10080, 5) dtypes: parcel_id int64 lga_source str zoning_class str area_m2 float64 geometry geometry dtype: object First 5 rows: parcel_id ... geometry 0 1 ... POLYGON ((540000 720000, 540010 720000, 540010... 1 2 ... POLYGON ((540010 720000, 540020 720000, 540020... 2 3 ... POLYGON ((540020 720000, 540030 720000, 540030... 3 4 ... POLYGON ((540030 720000, 540040 720000, 540040... 4 5 ... POLYGON ((540040 720000, 540050 720000, 540050... [5 rows x 5 columns]
0:10
Let me explore the data more - look at unique zoning_class values, blank entries, and vertex precision issues.
0:10
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import geopandas as gpd
gdf = gpd.read_file('lagos_zoning_legacy.gpkg')
print('Unique zoning_class values:')
vals = gdf['zoning_class'].unique()
for v in sorted(vals, key=str):
print(repr(v))
print(f'\nTotal unique: {len(vals)}')
print(f'\nBlank entries: {(gdf[\"zoning_class\"].isna() | gdf[\"zoning_class\"].str.strip().eq(\"\")).sum()}')
print(f'\nNaN count: {gdf[\"zoning_class\"].isna().sum()}')
"
Unique zoning_class values: '' '\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' Total unique: 28 Blank entries: 50 NaN count: 12
0:16
Let me understand the data: 1. **zoning_class values** - There are many variants: - Blank entries: '' (empty), '\t' (tab), ' ' (spaces), and NaN (12) - Residential variants: 'RESI.', 'RESIDENTIAL', 'Resi.', 'Residential', 'resi.', 'residential' - Commercial variants: 'COMM.', 'COMMERCIAL', 'Comm.', 'Commercial', 'comm.', 'commercial' - Industrial variants: 'INDUS.', 'INDUSTRIAL', 'Indus.', 'Industrial', 'indus.', 'industrial' - Agricultural variants: 'AGRI.', 'AGRICULTURAL', 'Agri.', 'Agricultural', 'agri.', 'agricultural' 2. **Canonical classes**: residential, commercial, industrial, agricultural (probably Title Case: Residential, Commercial, Industrial, Agricultural) 3. **Vertex precision issues** - causing tiny gaps/overlaps between adjacent parcels. Need to snap/round coordinates to a certain precision. 4. **Output**: One merged geometry per canonical zoning class, with two columns: zoning_class and area_m2. Let me look at a few more geometry details. The CRS is EPSG:26331 (Arc 1950 / UTM zone 31S). Coordinates are around 540000, 720000 which suggests UTM with 1m precision. The vertex precision issues likely involve coordinates with many decimal places. Let me explore the geometry precision.
0:16
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import geopandas as gpd
gdf = gpd.read_file('lagos_zoning_legacy.gpkg')
# Check some geometry coordinates for precision
import numpy as np
# Get all coordinates and see their decimal precision
coords = []
for geom in gdf.geometry.head(20):
if geom is not None:
for x, y in geom.exterior.coords:
coords.append((x, y))
print('Sample coords:')
for c in coords[:15]:
print(c)
# Check what decimal precision exists
xs = np.array([c[0] for c in coords])
ys = np.array([c[1] for c in coords])
print(f'\nX min decimals: {xs}')
print(f'Y min decimals: {ys}')
"
Sample coords: (540000.0, 720000.0) (540010.00002, 720000.000005) (540010.000005, 720010.000015) (540000.000025, 720010.00002) (540000.0, 720000.0) (540010.00001, 720000.000025) (540020.0, 720000.000005) (540020.00002, 720010.00001) (540010.000005, 720010.00002) (540010.00001, 720000.000025) (540020.000025, 720000.00002) (540030.000015, 720000.0) (540030.0, 720010.00001) (540020.00002, 720010.000015) (540020.000025, 720000.00002) X min decimals: [540000. 540010.00002 540010.000005 540000.000025 540000. 540010.00001 540020. 540020.00002 540010.000005 540010.00001 540020.000025 540030.000015 540030. 540020.00002 540020.000025 540030.000005 540040.000025 540040.000015 540030. 540030.000005 540040.00002 540050.000005 540050.00003 540040.000015 540040.00002 540050. 540060.00002 540060.00001 540050.00003 540050. 540060.000015 540070. 540070.00002 540060.00001 540060.000015 540070.00003 540080.000015 540080. 540070.000025 540070.00003 540080.00001 540090.00003 540090.000015 540080.000005 540080.00001 540090.000025 540100.00001 540100.00003 540090.000015 540090.000025 540100.000005 540110.000025 540110.00001 540100.00003 540100.000005 540110.000015 540120.000005 540120.000025 540110.00001 540110.000015 540120.00003 540130.00002 540130.000005 540120.000025 540120.00003 540130.00001 540140. 540140.00002 540130.000005 540130.00001 540140.000025 540150.00001 540150. 540140.00002 540140.000025 540150.000005 540160.000025 540160.000015 540150. 540150.000005 540160.00002 540170.000005 540170.000025 540160.000015 540160.00002 540170. 540180.00002 540180.000005 540170.00003 540170. 540180.000015 540190. 540190.00002 540180.00001 540180.000015 540190.00003 540200.000015 540200. 540190.00002 540190.00003 ] Y min decimals: [720000. 720000.000005 720010.000015 720010.00002 720000. 720000.000025 720000.000005 720010.00001 720010.00002 720000.000025 720000.00002 720000. 720010.00001 720010.000015 720000.00002 720000.00002 720000.000025 720010.000005 720010.000015 720000.00002 720000.000015 720000.000025 720010. 720010.00001 720000.000015 720000.000015 720000.00002 720010. 720010.000005 720000.000015 720000.00001 720000.00002 720010.000025 720010.000005 720000.00001 720000.000005 720000.000015 720010.000025 720010. 720000.000005 720000.000005 720000.00001 720010.00002 720010. 720000.000005 720000. 720000.00001 720010.000015 720010.000025 720000. 720000. 720000.000005 720010.000015 720010.00002 720000. 720000.000025 720000.000005 720010.00001 720010.00002 720000.000025 720000.00002 720000. 720010.00001 720010.000015 720000.00002 720000.00002 720000. 720010.000005 720010.000015 720000.00002 720000.000015 720000.000025 720010.000005 720010.00001 720000.000015 720000.000015 720000.00002 720010. 720010.00001 720000.000015 720000.00001 720000.00002 720010.000025 720010.000005 720000.00001 720000.00001 720000.000015 720010.000025 720010. 720000.00001 720000.000005 720000.000015 720010.00002 720010. 720000.000005 720000. 720000.00001 720010.00002 720010.000025 720000. ]
0:32
Let me understand the data: 1. **zoning_class cleanup**: There are 28 unique values including blank, whitespace-only, and NaN. The canonical classes are: - Agricultural → "Agricultural" - Commercial → "Commercial" - Industrial → "Industrial" - Residential → "Residential" Variations: 'AGRI.', 'AGRICULTURAL', 'Agri.', 'Agricultural', 'agri.', 'agricultural' → "Agricultural" Similar for other classes. Blank/whitespace-only/NaN entries should be dropped. 2. **Vertex precision**: Coordinates have up to 5-6 decimal places (like 0.000005, 0.000025) causing tiny gaps. I should round coordinates to a reasonable precision - since CRS is EPSG:26331 (A metric projection for Nigeria), rounding to e.g., 3 decimal places (millimeter precision) or even 2 (cm) or 0 (1m) would work. Let me check the scale... These look like they're in meters. Rounding to 0 decimal places (1m) would clean up most issues. But maybe 1 decimal place would be sufficient (0.1m = 10cm). Actually, let me think about this more carefully. The coordinates like 540000.0, 540010.00002 - the variations are tiny (0.00002 m = 0.02mm). So rounding to 1 or 2 decimal places would fix the precision issues while staying accurate enough. Let me use 1 decimal point (0.1m precision) or even 0 (1m precision). For the state portal, rounding to integer (1m) precision is probably fine for a zoning layer. The approach: - For each geometry, round all coordinates using shapely's `set_precision` or manually round coordinates. - Then use `dissolve` by zoning_class with `make_valid=True` (or snap/tiny gap fix). - Actually better: use `snap_to_grid` or `set_precision` with a grid size, then unary_union within each class. Let me use `shapely.set_precision` with a grid size of 0.01 or 0.1 to snap coordinates, then dissolve. Actually, the proper way to fix vertex precision issues causing tiny gaps and overlaps is to snap coordinates to a grid and then snap adjacent geometries. Let me use `snap_to_grid` or `set_precision` with a small grid size, then use a small snap tolerance when dissolving. Let me write a comprehensive script. 3. **Aggregation**: After cleaning classes and fixing geometry, dissolve by zoning_class into one merged geometry, compute area_m2. Let me write the script.
0:32
Write
/work/solve.py
content · 82 lines · py
#!/usr/bin/env python3
"""
Clean up lagos_zoning_legacy.gpkg:
1. Normalize zoning_class to canonical TitleCase values
2. Drop rows with blank/whitespace/NaN zoning_class
3. Fix vertex precision (snap to 0.01m grid)
4. Dissolve by zoning_class, compute area_m2
5. Write to zoning_aggregated.gpkg
"""
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely import set_precision, make_valid
# --- Read data ---
gdf = gpd.read_file('/work/lagos_zoning_legacy.gpkg')
print(f"Input rows: {len(gdf)}")
print(f"Input CRS: {gdf.crs}")
# --- Step 1: Clean zoning_class ---
# Build a normalisation map
norm_map = {}
for raw in gdf['zoning_class'].dropna().unique():
s = str(raw).strip().lower()
if s in ('',):
continue
if s.startswith('agri'):
norm_map[raw] = 'Agricultural'
elif s.startswith('comm'):
norm_map[raw] = 'Commercial'
elif s.startswith('indu'):
norm_map[raw] = 'Industrial'
elif s.startswith('resi'):
norm_map[raw] = 'Residential'
print(f"Normalisation map: {norm_map}")
# Apply normalisation
gdf['zoning_class'] = gdf['zoning_class'].map(lambda x: norm_map.get(x, None) if pd.notna(x) else None)
# Drop rows where zoning_class is None, empty, or whitespace-only
# (map would have returned None for unrecognised, and NaN stays NaN)
mask = gdf['zoning_class'].notna() & (gdf['zoning_class'].astype(str).str.strip() != '')
gdf = gdf[mask].copy()
print(f"Rows after dropping blanks: {len(gdf)}")
# Verify
print(f"Unique classes: {gdf['zoning_class'].unique()}")
# --- Step 2: Fix vertex precision ---
# Snap all coordinates to a 0.01 m grid to eliminate tiny gaps/overlaps
grid_size = 0.01
print(f"Snapping geometry to {grid_size} m grid ...")
gdf['geometry'] = gdf['geometry'].apply(lambda geom: set_precision(geom, grid_size, mode='valid_output'))
# Also ensure validity
gdf['geometry'] = gdf['geometry'].apply(lambda geom: make_valid(geom) if geom is not None else geom)
gdf = gdf[gdf['geometry'].notna() & ~gdf['geometry'].is_empty].copy()
print(f"Rows after geometry fix: {len(gdf)}")
# --- Step 3: Dissolve by zoning_class ---
print("Dissolving by zoning_class ...")
dissolved = gdf.dissolve(by='zoning_class', aggfunc='first')
# The dissolve result has zoning_class as index; reset to columns
dissolved = dissolved.reset_index()
# Make geometries valid and simplified (snapped)
dissolved['geometry'] = dissolved['geometry'].apply(lambda geom: make_valid(geom) if geom is not None else geom)
# Compute area in square metres
dissolved['area_m2'] = dissolved.geometry.area
# Keep only the two requested columns
result = dissolved[['zoning_class', 'area_m2', 'geometry']].copy()
print(f"Result rows: {len(result)}")
print(result[['zoning_class', 'area_m2']])
# --- Write output ---
result.to_file('/work/zoning_aggregated.gpkg', layer='zoning_aggregated', driver='GPKG', engine='pyogrio')
print("Written to /work/zoning_aggregated.gpkg")Wrote /work/solve.py (2938 chars).
0:37
Let me run the script to see if it works.
0:37
Bash
cd /work && python solve.py 2>&1
Input rows: 10080
Input CRS: EPSG:26331
Normalisation map: {'RESIDENTIAL': 'Residential', 'residential': 'Residential', 'Resi.': 'Residential', 'resi.': 'Residential', 'RESI.': 'Residential', 'Residential': 'Residential', 'Comm.': 'Commercial', 'comm.': 'Commercial', 'COMM.': 'Commercial', 'Commercial': 'Commercial', 'COMMERCIAL': 'Commercial', 'commercial': 'Commercial', 'Indus.': 'Industrial', 'indus.': 'Industrial', 'INDUS.': 'Industrial', 'Industrial': 'Industrial', 'INDUSTRIAL': 'Industrial', 'industrial': 'Industrial', 'AGRI.': 'Agricultural', 'Agricultural': 'Agricultural', 'AGRICULTURAL': 'Agricultural', 'agricultural': 'Agricultural', 'Agri.': 'Agricultural', 'agri.': 'Agricultural'}
Rows after dropping blanks: 10030
Unique classes: <ArrowStringArray>
['Residential', 'Commercial', 'Industrial', 'Agricultural']
Length: 4, dtype: str
Snapping geometry to 0.01 m grid ...
/work/solve.py:59: UserWarning: GeoSeries.notna() previously returned False for both missing (None) and empty geometries. Now, it only returns False for missing values. Since the calling GeoSeries contains empty geometries, the result has changed compared to previous versions of GeoPandas.
Given a GeoSeries 's', you can use '~s.is_empty & s.notna()' to get back the old behaviour.
To further ignore this warning, you can do:
import warnings; warnings.filterwarnings('ignore', 'GeoSeries.notna', UserWarning)
gdf = gdf[gdf['geometry'].notna() & ~gdf['geometry'].is_empty].copy()
Rows after geometry fix: 10000
Dissolving by zoning_class ...
Result rows: 4
zoning_class area_m2
0 Agricultural 250000.0
1 Commercial 250000.0
2 Industrial 250000.0
3 Residential 250000.0
Written to /work/zoning_aggregated.gpkg
0:44
The script ran successfully. Let me verify the output by reading it back.
0:44
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import geopandas as gpd
gdf = gpd.read_file('zoning_aggregated.gpkg')
print(gdf)
print()
print('CRS:', gdf.crs)
"
zoning_class area_m2 geometry 0 Agricultural 250000.0 POLYGON ((540500 720580, 540500 720590, 540500... 1 Commercial 250000.0 POLYGON ((540500 720080, 540500 720090, 540500... 2 Industrial 250000.0 POLYGON ((540000 720580, 540000 720590, 540000... 3 Residential 250000.0 POLYGON ((540000 720080, 540000 720090, 540000... CRS: EPSG:26331
1:04
The output looks good. Let me verify the data is clean by checking: 1. All 4 canonical classes are present (Agricultural, Commercial, Industrial, Residential) 2. The geometry is clean (snapped to 0.01m grid) 3. The area_m2 columns look reasonable (250,000 m² each) 4. CRS is preserved (EPSG:26331) Let me also do a quick sanity check on the geometry precision to confirm the snapping worked.
1:04
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import geopandas as gpd
gdf = gpd.read_file('zoning_aggregated.gpkg')
# Verify coordinates are snapped to 0.01
from shapely import get_coordinates
import numpy as np
coords = get_coordinates(gdf.geometry.values)
print('All coords divisible by 0.01?', np.allclose(coords % 0.01, 0, atol=1e-12))
print('Total coords:', len(coords))
print('Sample:', coords[:10])
"
All coords divisible by 0.01? False Total coords: 804 Sample: [[540500. 720580.] [540500. 720590.] [540500. 720600.] [540500. 720610.] [540500. 720620.] [540500. 720630.] [540500. 720640.] [540500. 720650.] [540500. 720660.] [540500. 720670.]]
1:08
The coordinates look like integer values (540500, 720580, etc.). Those are divisible by 0.01. Let me check more carefully - the check `% 0.01` with floating point arithmetic is tricky. Let me verify properly. Actually 540500 % 0.01 should be 0 since 540500 is an integer, which is divisible by 0.01. But floating point representation of 0.01 might cause issues. Let me check differently.
1:08
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import numpy as np
# Check if coordinates are multiples of 0.01
vals = np.array([540500., 720580., 540500., 720590.])
# Convert to centi-units
check = vals / 0.01
print('In 0.01 units:', check)
print('Close to int?', np.allclose(check, np.round(check), atol=1e-9))
"
In 0.01 units: [54050000. 72058000. 54050000. 72059000.] Close to int? True
1:13
The coordinates are indeed divisible by 0.01 — it's just a floating point comparison nuance. The values like 540500.0 are evenly divisible by 0.01. Let me verify more carefully. Actually, 540500.0 / 0.01 = 54050000.0 which is an integer. So the snapping worked. The `np.allclose` failure was because of floating point modulo arithmetic on already clean values — they are exact multiples. Let me double-check on the original data to be sure: Actually the original data had coordinates like 540010.00002, 720000.000005 etc. After snapping to 0.01, we got 540500.0, 720580.0 which are clean. The issue was just that my `% 0.01` check was imprecise due to floating point. The output file `zoning_aggregated.gpkg` is ready with 4 rows (Agricultural, Commercial, Industrial, Residential), each with clean geometry snapped to 0.01m, and area_m2 computed properly. Let me verify the output is complete and correct.
1:13
1:13
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import geopandas as gpd
gdf = gpd.read_file('zoning_aggregated.gpkg')
print('Columns:', gdf.columns.tolist())
print('Rows:', len(gdf))
print()
for _, row in gdf.iterrows():
print(f\"{row['zoning_class']:15s} {row['area_m2']:>10.1f} m² geom type: {row['geometry'].geom_type}\")
"
Columns: ['zoning_class', 'area_m2', 'geometry'] Rows: 4 Agricultural 250000.0 m² geom type: Polygon Commercial 250000.0 m² geom type: Polygon Industrial 250000.0 m² geom type: Polygon Residential 250000.0 m² geom type: Polygon
1:18
Everything looks good. Let me summarize what was done.
1:18