dc-l2-lagos-snap-normalize

Rerun

import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Polygon, MultiPolygon
from shapely.ops import unary_union, snap
from shapely.validation import make_valid

# ------------------------------------------------------------
# 1. Read data
# ------------------------------------------------------------
gdf = gpd.read_file('/work/lagos_zoning_legacy.gpkg')
print(f"Input rows: {len(gdf)}")

# ------------------------------------------------------------
# 2. Canonical zoning_class mapping
# ------------------------------------------------------------
class_map = {
    'residential': 'Residential',
    'resi.':       'Residential',
    'RESIDENTIAL': 'Residential',
    'RESI.':       'Residential',
    'Resi.':       'Residential',
    'commercial':  'Commercial',
    'comm.':       'Commercial',
    'COMMERCIAL':  'Commercial',
    'COMM.':       'Commercial',
    'Comm.':       'Commercial',
    'industrial':  'Industrial',
    'indus.':      'Industrial',
    'INDUSTRIAL':  'Industrial',
    'INDUS.':      'Industrial',
    'Indus.':      'Industrial',
    'agricultural':'Agricultural',
    'agri.':       'Agricultural',
    'AGRICULTURAL':'Agricultural',
    'AGRI.':       'Agricultural',
    'Agricultural':'Agricultural',
    'Agri.':       'Agricultural',
}

# Apply mapping (case-insensitive fallback for any unlisted)
def normalize_class(val):
    if val is None or (isinstance(val, str) and val.strip() == ''):
        return None
    v = str(val).strip()
    # Direct map lookup
    if v in class_map:
        return class_map[v]
    # Case-insensitive fallback: try lowercased
    v_lower = v.lower()
    for k, canon in class_map.items():
        if k.lower() == v_lower:
            return canon
    # Capitalize first letter as fallback
    return v.capitalize()

gdf['zoning_class'] = gdf['zoning_class'].apply(normalize_class)

# Drop rows with None (blank / whitespace / NaN)
before = len(gdf)
gdf = gdf.dropna(subset=['zoning_class'])
after = len(gdf)
print(f"Dropped {before - after} rows with blank/whitespace classes")
print(f"Remaining: {len(gdf)}")

# Check distribution
print("\nCanonical class counts:")
print(gdf['zoning_class'].value_counts())

# ------------------------------------------------------------
# 3. Fix vertex precision – snap to a 0.001 grid (mm)
# ------------------------------------------------------------
# Use Shapely's snap + coordinate rounding to eliminate micro-gaps
# First round coordinates to 3 decimal places
def round_geom(geom, ndigits=3):
    """Round all coordinates in a geometry to fixed decimal places."""
    if geom is None:
        return None
    return shapely.set_precision(geom, grid_size=10**-ndigits)

import shapely

gdf['geometry'] = gdf['geometry'].apply(lambda g: round_geom(g, ndigits=3))

# ------------------------------------------------------------
# 4. Snap adjacent geometries within a tiny tolerance
# ------------------------------------------------------------
# Dissolve by class first (temporary) to get a merged per-class shape
# But we need to fix overlaps/gaps between classes too.
# Strategy: use unary_union per class after snapping, which will
# cleanly merge overlapping/snapped boundaries.

# Ensure valid geometries
gdf['geometry'] = gdf['geometry'].apply(lambda g: make_valid(g) if g is not None else g)

# ------------------------------------------------------------
# 5. Dissolve by zoning_class – merge geometries
# ------------------------------------------------------------
dissolved = gdf.dissolve(by='zoning_class', aggfunc='mean')

# ------------------------------------------------------------
# 6. Compute area_m2 for each merged geometry
# ------------------------------------------------------------
# CRS is EPSG:26331 (projected, meters) – area from geometry is in m²
dissolved['area_m2'] = dissolved.geometry.area

# Round area to 2 decimal places for cleanliness
dissolved['area_m2'] = dissolved['area_m2'].round(2)

# Keep only the two requested columns
result = dissolved[['area_m2']].copy()
result.index.name = 'zoning_class'
result = result.reset_index()

print("\nResult:")
print(result)
print(f"\nTotal area: {result['area_m2'].sum():.2f} m²")

# ------------------------------------------------------------
# 7. Write output
# ------------------------------------------------------------
result.to_file('/work/zoning_aggregated.gpkg', layer='zoning_aggregated', driver='GPKG')
print("\nWritten to /work/zoning_aggregated.gpkg")

import geopandas as gpd
import pandas as pd
import numpy as np
import shapely
from shapely.validation import make_valid

# ------------------------------------------------------------
# 1. Read data
# ------------------------------------------------------------
gdf = gpd.read_file('/work/lagos_zoning_legacy.gpkg')
print(f"Input rows: {len(gdf)}")

# ------------------------------------------------------------
# 2. Canonical zoning_class mapping
# ------------------------------------------------------------
class_map = {
    'residential': 'Residential',
    'resi.':       'Residential',
    'RESIDENTIAL': 'Residential',
    'RESI.':       'Residential',
    'Resi.':       'Residential',
    'commercial':  'Commercial',
    'comm.':       'Commercial',
    'COMMERCIAL':  'Commercial',
    'COMM.':       'Commercial',
    'Comm.':       'Commercial',
    'industrial':  'Industrial',
    'indus.':      'Industrial',
    'INDUSTRIAL':  'Industrial',
    'INDUS.':      'Industrial',
    'Indus.':      'Industrial',
    'agricultural':'Agricultural',
    'agri.':       'Agricultural',
    'AGRICULTURAL':'Agricultural',
    'AGRI.':       'Agricultural',
    'Agricultural':'Agricultural',
    'Agri.':       'Agricultural',
}

def normalize_class(val):
    if val is None:
        return None
    if isinstance(val, str) and val.strip() == '':
        return None
    # Also catch 'nan' strings or blank-only strings
    v = str(val).strip()
    if v.lower() == 'nan' or v == '':
        return None
    # Direct map lookup
    if v in class_map:
        return class_map[v]
    # Case-insensitive fallback
    v_lower = v.lower()
    for k, canon in class_map.items():
        if k.lower() == v_lower:
            return canon
    # Fallback: TitleCase
    return v.capitalize()

gdf['zoning_class'] = gdf['zoning_class'].apply(normalize_class)

# Drop rows with None (blank / whitespace / NaN)
before = len(gdf)
gdf = gdf.dropna(subset=['zoning_class'])
after = len(gdf)
print(f"Dropped {before - after} rows with blank/whitespace classes")
print(f"Remaining: {len(gdf)}")

# Check distribution
print("\nCanonical class counts:")
print(gdf['zoning_class'].value_counts())

# ------------------------------------------------------------
# 3. Fix vertex precision – snap to a 0.001 grid
# ------------------------------------------------------------
def snap_to_grid(geom, grid_size=0.001):
    """Round all coordinates in a geometry to a fixed grid."""
    if geom is None:
        return None
    return shapely.set_precision(geom, grid_size=grid_size, mode='valid_output')

gdf['geometry'] = gdf['geometry'].apply(lambda g: snap_to_grid(g, 0.001))

# Ensure valid geometries
gdf['geometry'] = gdf['geometry'].apply(lambda g: make_valid(g) if g is not None else g)

# ------------------------------------------------------------
# 4. Dissolve by zoning_class – merge geometries
# ------------------------------------------------------------
dissolved = gdf.dissolve(by='zoning_class', aggfunc='mean')

# ------------------------------------------------------------
# 5. Compute area_m2 for each merged geometry
# ------------------------------------------------------------
# CRS is EPSG:26331 (projected, metres) – area from geometry is in m²
dissolved['area_m2'] = dissolved.geometry.area.round(2)

# Keep only the two requested columns
result = dissolved[['area_m2']].copy()
result.index.name = 'zoning_class'
result = result.reset_index()

print("\nResult:")
print(result)
print(f"\nTotal area: {result['area_m2'].sum():.2f} m²")

# ------------------------------------------------------------
# 6. Write output
# ------------------------------------------------------------
result.to_file('/work/zoning_aggregated.gpkg', layer='zoning_aggregated', driver='GPKG')
print("\nWritten to /work/zoning_aggregated.gpkg")

import geopandas as gpd
import pandas as pd
import numpy as np
import shapely
from shapely.validation import make_valid

# ------------------------------------------------------------
# 1. Read data
# ------------------------------------------------------------
gdf = gpd.read_file('/work/lagos_zoning_legacy.gpkg')
print(f"Input rows: {len(gdf)}")

# ------------------------------------------------------------
# 2. Canonical zoning_class mapping
# ------------------------------------------------------------
class_map = {
    'residential': 'Residential',
    'resi.':       'Residential',
    'RESIDENTIAL': 'Residential',
    'RESI.':       'Residential',
    'Resi.':       'Residential',
    'commercial':  'Commercial',
    'comm.':       'Commercial',
    'COMMERCIAL':  'Commercial',
    'COMM.':       'Commercial',
    'Comm.':       'Commercial',
    'industrial':  'Industrial',
    'indus.':      'Industrial',
    'INDUSTRIAL':  'Industrial',
    'INDUS.':      'Industrial',
    'Indus.':      'Industrial',
    'agricultural':'Agricultural',
    'agri.':       'Agricultural',
    'AGRICULTURAL':'Agricultural',
    'AGRI.':       'Agricultural',
    'Agricultural':'Agricultural',
    'Agri.':       'Agricultural',
}

def normalize_class(val):
    if val is None:
        return None
    if isinstance(val, str) and val.strip() == '':
        return None
    v = str(val).strip()
    if v.lower() == 'nan' or v == '':
        return None
    # Direct map lookup
    if v in class_map:
        return class_map[v]
    # Case-insensitive fallback
    v_lower = v.lower()
    for k, canon in class_map.items():
        if k.lower() == v_lower:
            return canon
    # Fallback: TitleCase
    return v.capitalize()

gdf['zoning_class'] = gdf['zoning_class'].apply(normalize_class)

# Drop rows with None (blank / whitespace / NaN)
before = len(gdf)
gdf = gdf.dropna(subset=['zoning_class'])
after = len(gdf)
print(f"Dropped {before - after} rows with blank/whitespace classes")
print(f"Remaining: {len(gdf)}")

print("\nCanonical class counts:")
print(gdf['zoning_class'].value_counts())

# ------------------------------------------------------------
# 3. Fix vertex precision – snap to a 0.001 grid
# ------------------------------------------------------------
def snap_to_grid(geom, grid_size=0.001):
    """Round all coordinates in a geometry to a fixed grid."""
    if geom is None:
        return None
    return shapely.set_precision(geom, grid_size=grid_size, mode='valid_output')

gdf['geometry'] = gdf['geometry'].apply(lambda g: snap_to_grid(g, 0.001))

# Ensure valid geometries
gdf['geometry'] = gdf['geometry'].apply(lambda g: make_valid(g) if g is not None else g)

# ------------------------------------------------------------
# 4. Dissolve by zoning_class – merge geometries
# ------------------------------------------------------------
# Keep only essential columns for dissolve
gdf_simple = gdf[['zoning_class', 'geometry']].copy()
dissolved = gdf_simple.dissolve(by='zoning_class')

# ------------------------------------------------------------
# 5. Compute area_m2 for each merged geometry
# ------------------------------------------------------------
# CRS is EPSG:26331 (projected, metres) – area from geometry is in m²
dissolved['area_m2'] = dissolved.geometry.area.round(2)

# Keep only the two requested columns
result = dissolved[['area_m2']].copy()
result.index.name = 'zoning_class'
result = result.reset_index()

print("\nResult:")
print(result)
print(f"\nTotal area: {result['area_m2'].sum():.2f} m²")

# ------------------------------------------------------------
# 6. Write output
# ------------------------------------------------------------
result.to_file('/work/zoning_aggregated.gpkg', layer='zoning_aggregated', driver='GPKG')
print("\nWritten to /work/zoning_aggregated.gpkg")