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 before cleaning:")
    print(gdf['zoning_class'].unique())
    print("\nGeometry type:", gdf.geometry.type.unique())
    print("CRS:", 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 Polygon, MultiPolygon
from shapely.ops import unary_union

def clean_zoning_class(val):
    if pd.isna(val):
        return None
    s = str(val).strip().lower()
    if not s:
        return None
    
    # Define mappings for spelling variations
    mappings = {
        'resi.': 'Residential',
        'residential': 'Residential',
        'comm.': 'Commercial',
        'commercial': 'Commercial',
        'indus.': 'Industrial',
        'industrial': 'Industrial',
        'agri.': 'Agricultural',
        'agricultural': 'Agricultural'
    }
    
    # Try to find mapping
    for key, canonical in mappings.items():
        if s == key:
            return canonical
    
    # Fallback to TitleCase if no specific mapping found
    return s.capitalize()

def solve():
    # 1. Load data
    gdf = gpd.read_file('/work/lagos_zoning_legacy.gpkg')
    print("Loaded input data.")

    # 2. Clean zoning_class
    # First, apply the cleaning function
    gdf['zoning_class'] = gdf['zoning_class'].apply(clean_zoning_class)
    
    # Drop rows where zoning_class is None (handles NaN, '', ' ', etc.)
    gdf = gdf.dropna(subset=['zoning_class'])
    gdf = gdf[gdf['zoning_class'].str.strip() != '']
    
    print("Unique zoning classes after cleaning:", gdf['zoning_class'].unique())

    # 3. Fix vertex precision and topological issues
    # We use a small buffer(0) trick or snap to grid to fix gaps/overlaps
    # But since we are merging by class, we can unary_union the geometries of each class.
    # unary_union is very robust at fixing overlaps and gaps if we first snap.
    
    # Snap geometries to a grid to reduce precision artifacts
    # EPSG:26331 is metric. 1cm precision (0.01m) should be plenty for zoning.
    precision = 0.01
    gdf['geometry'] = gdf['geometry'].apply(lambda geom: geom.simplify(0, preserve_topology=True)) # placeholder
    # A better way for precision:
    from shapely import wkt
    
    # Re-snapping to grid
    def snap_to_grid(geom, precision):
        # This is a simple way to snap coordinates to a grid
        if geom.is_empty:
            return geom
        return gpd.GeoSeries([geom]).apply(lambda g: g.buffer(0)).iloc[0] # buffer(0) fixes some topology

    # Actually, the most robust way to "fix" gaps/overlaps is to dissolve the geometries.
    # But to ensure the "gaps" are filled, we might need to merge everything and then subtract?
    # No, the prompt says: "fix the vertex precision problems, and produce one merged geometry per canonical zoning class"
    # This implies the classes are somewhat well-behaved but have slivers.
    
    # Let's snap to grid first.
    def snap_coords(geom, grid_size):
        from shapely.affinity import scale, translate
        # A more direct way to snap coordinates:
        def snap_point(x, y):
            return round(x / grid_size) * grid_size, round(y / grid_size) * grid_size
        
        # We can use shapely.ops.transform
        from shapely.ops import transform
        return transform(lambda x, y: snap_point(x, y), geom)

    print("Snapping geometries to grid...")
    gdf['geometry'] = gdf['geometry'].apply(lambda g: snap_coords(g, precision))

    # 4. Merge geometries per zoning class
    print("Dissolving geometries by zoning_class...")
    # Dissolve will merge all geometries in a class into one (or many MultiPolygons if they aren't contiguous)
    # The prompt says "produce one merged geometry per canonical zoning class"
    # This usually means a single MultiPolygon per class.
    
    dissolved = gdf.dissolve(by='zoning_class', aggfunc={'area_m2': 'sum'})
    
    # Note: dissolve might create separate polygons if they aren't touching. 
    # If the user wants exactly ONE row per class, we need to ensure they are MultiPolygons.
    # dissolve() with a single geometry output per group usually returns MultiPolygons if they aren't contiguous.
    
    # Check if we need to force single geometry per class (it should be MultiPolygon already)
    # Let's check types
    print("Dissolved geometry types:", dissolved.geometry.type.unique())

    # 5. Final check on area_m2
    # The prompt asks for area_m2. The input has area_m2.
    # If we dissolve, we should probably recalculate the area to be safe,
    # or use the sum of input areas if we assume no overlaps.
    # However, since we are fixing overlaps/gaps, calculating from the new geometry is more accurate.
    
    dissolved['area_m2'] = dissolved.geometry.area

    # 6. Finalizing columns
    # Requirement: "each row carrying two columns: zoning_class and a numeric area_m2"
    # dissolve() puts 'zoning_class' in the index.
    dissolved = dissolved.reset_index()
    result = dissolved[['zoning_class', 'area_m2', 'geometry']]
    
    # The prompt specifically says "each row carrying two columns: zoning_class and a numeric area_m2"
    # It doesn't explicitly say "exclude geometry", but you can't have a GeoPackage without geometry.
    # I will provide the geometry as it is a GeoPackage.
    
    # Final selection as per instruction
    # "produce one merged geometry per canonical zoning class... with each row carrying two columns: zoning_class and a numeric area_m2"
    # This is a bit ambiguous about whether the geometry column should be included in the final file.
    # In a GPKG, the geometry column is mandatory.
    
    print("Saving to zoning_aggregated.gpkg...")
    result.to_file('/work/zoning_aggregated.gpkg', driver='GPKG')
    print("Done.")

if __name__ == "__main__":
    solve()

import geopandas as gpd
import pandas as pd
from shapely.ops import transform
import numpy as np

def clean_zoning_class(val):
    if pd.isna(val):
        return None
    s = str(val).strip().lower()
    if not s:
        return None
    
    mappings = {
        'resi.': 'Residential',
        'residential': 'Residential',
        'comm.': 'Commercial',
        'commercial': 'Commercial',
        'indus.': 'Industrial',
        'industrial': 'Industrial',
        'agri.': 'Agricultural',
        'agricultural': 'Agricultural'
    }
    
    if s in mappings:
        return mappings[s]
    
    # Fallback: try to handle the case where it might be 'Resi' or something similar
    # but let's stick to the provided mapping first.
    # If not in mapping, try title case
    return s.title()

def snap_coords(geom, grid_size):
    def snap_point(x, y):
        return round(x / grid_size) * grid_size, round(y / grid_size) * grid_size
    return transform(lambda x, y: snap_point(x, y), geom)

def solve():
    # 1. Load data
    gdf = gpd.read_file('/work/lagos_zoning_legacy.gpkg')
    print(f"Loaded {len(gdf)} rows.")

    # 2. Clean zoning_class
    gdf['zoning_class'] = gdf['zoning_class'].apply(clean_zoning_class)
    gdf = gdf.dropna(subset=['zoning_class'])
    # Filter out any empty strings that might have survived
    gdf = gdf[gdf['zoning_class'].str.strip() != '']
    
    print("Unique zoning classes after cleaning:", gdf['zoning_class'].unique())

    # 3. Fix vertex precision and topological issues
    # Snapping to 1cm (0.01m) to fix tiny gaps and overlaps.
    precision = 0.01
    print("Snapping geometries to grid...")
    gdf['geometry'] = gdf['geometry'].apply(lambda g: snap_coords(g, precision))

    # 4. Merge geometries per zoning class
    # We use unary_union for each group to dissolve and clean up overlaps/gaps.
    # Dissolve with a groupby/agg is standard, but unary_union is better for cleaning geometry.
    
    print("Aggregating geometries...")
    
    results = []
    for cls, group in gdf.groupby('zoning_class'):
        # unary_union dissolves all geometries in the group into a single geometry
        merged_geom = group.geometry.unary_union
        # Recalculate area from the merged geometry to be precise
        merged_area = merged_geom.area
        results.append({
            'zoning_class': cls,
            'area_m2': merged_area,
            'geometry': merged_geom
        })
    
    # Create the final GeoDataFrame
    result_gdf = gpd.GeoDataFrame(results, crs=gdf.crs)
    
    # Ensure 'zoning_class' is a string (it might have become an object/category)
    result_gdf['zoning_class'] = result_gdf['zoning_class'].astype(str)
    result_gdf['area_m2'] = result_gdf['area_m2'].astype(float)

    print("Finalizing result...")
    print(result_gdf[['zoning_class', 'area_m2']])
    
    # 5. Save the result
    result_gdf.to_file('/work/zoning_aggregated.gpkg', driver='GPKG')
    print("Successfully saved to /work/zoning_aggregated.gpkg")

if __name__ == "__main__":
    solve()