dc-l2-lagos-snap-normalize

Rerun

import geopandas as gpd

def inspect_data(file_path):
    gdf = gpd.read_file(file_path)
    print("Columns:", gdf.columns.tolist())
    print("\nValue counts for 'zoning_class':")
    print(gdf['zoning_class'].value_counts(dropna=False))
    print("\nGeometry type:", gdf.geometry.type.unique())
    print("\nShape:", gdf.shape)
    print("\nCRS:", gdf.crs)
    return gdf

if __name__ == "__main__":
    gdf = inspect_data('lagos_zoning_legacy.gpkg')

import geopandas as gpd
import pandas as pd
from shapely.geometry import Polygon, MultiPolygon
from shapely.ops import unary_union

def clean_zoning_class(val):
    if not isinstance(val, str):
        return None
    s = val.strip().lower()
    if not s:
        return None
    
    # Mapping variations to canonical TitleCase names
    # Based on inspection: Residential, Commercial, Industrial, Agricultural
    if s in ['residential', 'resi.', 'resi']:
        return 'Residential'
    elif s in ['commercial', 'comm.', 'comm']:
        return 'Commercial'
    elif s in ['industrial', 'indus.', 'indus']:
        return 'Industrial'
    elif s in ['agricultural', 'agri.', 'agri']:
        return 'Agricultural'
    else:
        # Fallback: just Title Case if it's unknown but not empty
        return s.capitalize()

def solve():
    # 1. Load data
    gdf = gpd.read_file('lagos_zoning_legacy.gpkg')
    
    # 2. Clean zoning_class
    # Convert to string first to handle potential mixed types, then apply cleaning
    gdf['zoning_class'] = gdf['zoning_class'].astype(str).apply(clean_zoning_class)
    
    # Drop rows with None (which covers NaN, empty, and whitespace-only after cleaning)
    gdf = gdf.dropna(subset=['zoning_class'])
    
    # 3. Fix vertex precision and resolve gaps/overlaps
    # A common trick for "fixing" gaps/overlaps in vector data is to round coordinates 
    # and then dissolve. To handle the "gaps and overlaps", we can use a small buffer 
    # or simply rely on the dissolve which handles overlaps by merging geometries.
    # However, 'tiny gaps' usually require a small buffer or snapping.
    # We'll use a small buffer(0) or simplify to snap vertices, but a more robust 
    # way for "merging per class" is to dissolve.
    
    # For precision issues, let's round the coordinates to a reasonable precision (e.g., 2 decimal places for meters)
    # Given CRS is EPSG:26331 (meters), 2 decimal places is ~1cm precision.
    def snap_to_grid(geom, precision=0.01):
        # Using a simple approach: round coordinates
        # shapely's simplify can also help, but rounding is more deterministic for "precision issues"
        # However, the easiest way to fix "tiny gaps/overlaps" is to dissolve.
        # If gaps exist, dissolve will create a MultiPolygon with holes.
        # If they must be "fixed", we could buffer by a tiny amount then dissolve.
        return geom.buffer(0.01).buffer(-0.01) # This can be risky.
    
    # Let's try a safer approach: dissolve first, then clean up. 
    # The prompt says "fix the vertex precision problems, and produce one merged geometry per canonical zoning class".
    # We will round the coordinates to fix precision.
    
    def round_coords(geom, precision=0.01):
        if geom.is_empty:
            return geom
        # Use a more standard way to handle vertex precision: precision rounding via simplify
        return geom.simplify(precision, preserve_topology=True)

    gdf['geometry'] = gdf['geometry'].apply(lambda x: round_coords(x, 0.01))

    # 4. Aggregate (Dissolve)
    # Dissolve by 'zoning_class'
    # We need to calculate area_m2 for the aggregated geometry
    
    # Group by zoning_class and dissolve
    dissolved = gdf.dissolve(by='zoning_class', aggfunc={'area_m2': 'sum'})
    
    # Note: the 'dissolve' method in geopandas returns a GeoDataFrame where 
    # the grouping column becomes the index. We need to reset it and ensure 
    # 'area_m2' is the calculated area of the *new* geometry.
    
    # Let's re-calculate area for the final merged geometries to be sure.
    dissolved = dissolved.reset_index()
    dissolved['area_m2'] = dissolved.geometry.area
    
    # Keep only the required columns
    result = dissolved[['zoning_class', 'area_m2', 'geometry']]
    
    # 5. Write result
    result.to_file('zoning_aggregated.gpkg', driver='GPKG')
    print("Successfully created zoning_aggregated.gpkg")
    print(result)

if __name__ == "__main__":
    solve()