dc-l2-lagos-snap-normalize

Rerun

import geopandas as gpd

def inspect_data(path):
    gdf = gpd.read_file(path)
    print("Columns:", gdf.columns.tolist())
    print("Head:\n", gdf.head())
    print("\nUnique zoning_class values:\n", gdf['zoning_class'].unique())
    print("\nGeometry type:", gdf.geometry.type.unique())
    print("\nCRS:", gdf.crs)
    return gdf

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

import geopandas as gpd

def inspect_data(path):
    gdf = gpd.read_file(path)
    print("Columns:", gdf.columns.tolist())
    print("Head:\n", gdf.head())
    print("\nUnique zoning_class values:\n", gdf['zoning_class'].unique())
    print("\nGeometry type:", gdf.geometry.type.unique())
    print("\nCRS:", gdf.crs)
    return gdf

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

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

def clean_zoning_class(val):
    if pd.isna(val) or str(val).strip() == "":
        return None
    
    s = str(val).strip().lower()
    
    # Mapping variations to canonical names
    # Based on inspection:
    # RESI (Residential)
    # COMM (Commercial)
    # INDUS (Industrial)
    # AGRI (Agricultural)
    
    if s.startswith(('resi', 'res')):
        return "Residential"
    elif s.startswith(('comm', 'com')):
        return "Commercial"
    elif s.startswith(('indus', 'ind')):
        return "Industrial"
    elif s.startswith(('agri', 'agr')):
        return "Agricultural"
    else:
        # Fallback to Title Case if it doesn't match prefixes but is a valid string
        return s.title()

def solve():
    input_path = '/work/lagos_zoning_legacy.gpkg'
    output_path = '/work/zoning_aggregated.gpkg'
    
    print(f"Reading {input_path}...")
    gdf = gpd.read_file(input_path)
    
    # 1. Clean zoning_class
    print("Cleaning zoning_class...")
    gdf['zoning_class'] = gdf['zoning_class'].apply(clean_zoning_class)
    
    # 2. Drop rows with blank or whitespace-only classes (which are now None)
    gdf = gdf.dropna(subset=['zoning_class'])
    
    # 3. Fix vertex precision issues
    # A common way to fix tiny gaps/overlaps is to round coordinates and then 
    # use a small buffer(0) or snap to grid.
    # Given the coordinate magnitudes (~540000), rounding to 1 or 2 decimals 
    # (cm/mm level) should be safe.
    print("Fixing vertex precision...")
    def snap_to_grid(geom, precision=0.01):
        # snap to grid by rounding coordinates
        if geom.is_empty:
            return geom
        return gpd.GeoSeries([geom]).simplify(precision, preserve_topology=True).iloc[0]

    # Instead of simplify, let's try rounding coordinates via a more robust method if possible,
    # but simplify(preserve_topology=True) is quite effective for cleaning up small artifacts.
    # However, a common GIS trick for "cleaning" geometry is buffer(0).
    # Let's use a combination: snap to grid via rounding (if we can) or simplify.
    # Since we can't easily round all coordinates in Shapely directly without iterating, 
    # let's use simplify with a very small tolerance to snap vertices.
    gdf['geometry'] = gdf['geometry'].simplify(0.01, preserve_topology=True)
    
    # Ensure geometries are valid
    gdf['geometry'] = gdf['geometry'].buffer(0)
    
    # 4. Produce one merged geometry per canonical zoning class
    print("Aggregating by zoning class...")
    
    # We need to group by zoning_class, dissolve the geometries, and calculate area.
    # Note: Dissolving might create small gaps if the precision fix wasn't perfect.
    # But the task asks for "one merged geometry per canonical zoning class".
    
    # We will group by zoning_class
    grouped = gdf.dissolve(by='zoning_class', aggfunc={'area_m2': 'sum'})
    
    # The dissolve operation keeps the column used for grouping.
    # The resulting gdf will have 'zoning_class' as index.
    # We need to reset index and ensure columns are correct.
    
    result_gdf = grouped.reset_index()
    
    # Re-calculate area_m2 on the dissolved geometries to ensure accuracy 
    # (in case the original area_m2 column was just a field and not the actual area 
    # of the geometry, or if dissolving changed things).
    # The prompt says: "each row carrying two columns: zoning_class and a numeric area_m2"
    # It's safer to calculate it from the new geometries.
    
    result_gdf['area_m2'] = result_gdf['geometry'].area
    
    # Final selection of columns
    result_gdf = result_gdf[['zoning_class', 'area_m2', 'geometry']]
    
    print(f"Writing to {output_path}...")
    result_gdf.to_file(output_path, driver="GPKG")
    
    print("Done!")
    print(result_gdf.head())

if __name__ == "__main__":
    solve()

import geopandas as gpd
import pandas as pd

def clean_zoning_class(val):
    if pd.isna(val) or str(val).strip() == "":
        return None
    
    s = str(val).strip().lower()
    
    if s.startswith(('resi', 'res')):
        return "Residential"
    elif s.startswith(('comm', 'com')):
        return "Commercial"
    elif s.startswith(('indus', 'ind')):
        return "Industrial"
    elif s.startswith(('agri', 'agr')):
        return "Agricultural"
    else:
        return s.title()

def solve():
    input_path = '/work/lagos_zoning_legacy.gpkg'
    output_path = '/work/zoning_aggregated.gpkg'
    
    gdf = gpd.read_file(input_path)
    
    # 1. Clean zoning_class and drop blanks
    gdf['zoning_class'] = gdf['zoning_class'].apply(clean_zoning_class)
    gdf = gdf.dropna(subset=['zoning_class'])
    
    # 2. Fix vertex precision issues (gaps/overlaps)
    # First, snap to a grid to ensure adjacent vertices align
    # We'll use a small tolerance for simplification to snap vertices.
    # Then buffer(0) to fix topology.
    gdf['geometry'] = gdf['geometry'].simplify(0.01, preserve_topology=True)
    gdf['geometry'] = gdf['geometry'].buffer(0)
    
    # 3. Aggregate by zoning class
    # We dissolve the geometries to merge them into one per class
    # and sum the area.
    # We use dissolve on 'zoning_class'.
    grouped = gdf.dissolve(by='zoning_class')
    
    # 4. Final cleaning of output
    # The prompt asks for 'zoning_class' and 'area_m2'
    # The dissolved gdf index is 'zoning_class'.
    result_gdf = grouped.reset_index()
    result_gdf['area_m2'] = result_gdf['geometry'].area
    
    # Ensure the output has exactly the requested columns
    final_gdf = result_gdf[['zoning_class', 'area_m2', 'geometry']]
    
    # Write to GeoPackage
    final_gdf.to_file(output_path, driver="GPKG")
    print("Success")

if __name__ == "__main__":
    solve()