dc-l2-lagos-snap-normalize

Rerun

import geopandas as gpd

gdf = gpd.read_file('lagos_zoning_legacy.gpkg')
print(gdf.head())
print(gdf.info())
print(gdf['zoning_class'].unique())

import geopandas as gpd
import numpy as np

gdf = gpd.read_file('lagos_zoning_legacy.gpkg')
print("Columns:", gdf.columns)
print("First few rows:")
print(gdf.head())
print("\nUnique zoning classes:")
print(gdf['zoning_class'].unique())

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

def clean_zoning_class(val):
    if pd.isna(val) or str(val).strip() == "":
        return None
    
    val = str(val).strip().lower()
    
    # Mapping variations to canonical TitleCase values
    mapping = {
        'residential': 'Residential',
        'resi.': 'Residential',
        'res': 'Residential',
        'commercial': 'Commercial',
        'comm.': 'Commercial',
        'comm': 'Commercial',
        'industrial': 'Industrial',
        'indus.': 'Industrial',
        'indus': 'Industrial',
        'agricultural': 'Agricultural',
        'agri.': 'Agricultural',
        'agri': 'Agricultural'
    }
    
    # Try to find a match in the mapping, else TitleCase it
    # We use common prefixes to catch variations like 'Resi.'
    for key, canonical in mapping.items():
        if val.startswith(key):
            return canonical
            
    return val.capitalize()

def solve():
    print("Loading data...")
    gdf = gpd.read_file('lagos_zoning_legacy.gpkg')
    
    print("Cleaning zoning classes...")
    gdf['zoning_class'] = gdf['zoning_class'].apply(clean_zoning_class)
    
    # Drop rows with None/NaN/Empty zoning_class
    gdf = gdf.dropna(subset=['zoning_class'])
    gdf = gdf[gdf['zoning_class'].str.strip() != ""]
    
    print(f"Unique classes after cleaning: {gdf['zoning_class'].unique()}")
    
    # Fix vertex precision and overlaps/gaps
    # A common technique is to snap geometries to a grid or to each other.
    # Given the requirement "fix the vertex precision problems... produce one merged geometry per canonical zoning class"
    # We will group by zoning_class and use unary_union, which dissolves boundaries.
    # However, if there are small gaps, unary_union won't fill them.
    # We can snap geometries to a grid to align vertices.
    
    print("Snapping geometries to a grid to fix precision issues...")
    # Using a precision of 0.1 meters (assuming the CRS is in meters)
    # This helps align vertices that are slightly off.
    gdf['geometry'] = gdf['geometry'].buffer(0) # Fix invalid geometries
    
    # To address gaps/overlaps before union, we can use a small buffer or snap.
    # Let's try snapping all geometries to a 0.01 grid.
    # But for simplicity and robustness in a 'cleaning' task, 
    # we can group by class and use unary_union. 
    # To handle gaps, we can apply a tiny buffer then dissolve, but that might overlap.
    # Better approach: snap vertices to a grid.
    
    from shapely import wkb
    # Snap to 1cm grid
    precision = 0.01
    gdf['geometry'] = gdf['geometry'].apply(lambda geom: snap(geom, geom, precision)) # This is identity, not helpful.
    
    # Let's use a standard trick: buffer(0) and then unary_union by group.
    # To handle the "tiny gaps/overlaps" explicitly, we can use a very small buffer.
    # Or, we can snap the whole layer to a grid.
    
    # Actually, the requirement says "fix the vertex precision problems, and produce one merged geometry per canonical zoning class".
    # A common way to fix precision is to use the 'shapely.set_precision' if available, or coordinate rounding.
    
    def snap_to_grid(geom, grid_size=0.01):
        # Rounding coordinates to the grid size
        if geom.is_empty:
            return geom
        return lambda x, y: (round(x / grid_size) * grid_size, round(y / grid_size) * grid_size)
    
    # More robust: use shapely.set_precision (available in newer shapely)
    try:
        from shapely import set_precision
        gdf['geometry'] = gdf['geometry'].apply(lambda g: set_precision(g, grid_size=0.01))
    except ImportError:
        print("shapely.set_precision not available, using coordinate rounding.")
        def round_coords(geom):
            from shapely.ops import transform
            return transform(lambda x, y: (round(x, 2), round(y, 2)), geom)
        gdf['geometry'] = gdf['geometry'].apply(round_coords)

    print("Dissolving geometries by zoning class...")
    # Dissolve by zoning_class
    # We want to merge geometries first, then calculate area.
    # If there are gaps, unary_union on the snapped geometries should close many of them.
    
    aggregated_data = []
    
    for z_class, group in gdf.groupby('zoning_class'):
        merged_geom = group['geometry'].unary_union
        # If gaps still exist, unary_union won't fill them. 
        # But typically in these tasks, snapping to grid is the intended solution for "precision issues".
        
        area = merged_geom.area
        aggregated_data.append({
            'zoning_class': z_class,
            'area_m2': area,
            'geometry': merged_geom
        })
        
    print("Creating output GeoDataFrame...")
    out_gdf = gpd.GeoDataFrame(aggregated_data, crs=gdf.crs)
    
    # Ensure columns are in correct order and only contains requested
    # "each row carrying two columns: zoning_class and a numeric area_m2"
    # Wait, it says "produce one merged geometry per canonical zoning class"
    # and then "each row carrying two columns: zoning_class and a numeric area_m2".
    # Usually, a GeoDataFrame must have a geometry column. I'll include it.
    # If the user strictly wants ONLY two columns, they might want a non-spatial file, 
    # but it's a .gpkg, so it must have geometry. I will include the geometry.
    
    # Re-reading: "produce one merged geometry per canonical zoning class? ... write the result to zoning_aggregated.gpkg, with each row carrying two columns: zoning_class and a numeric area_m2."
    # This is slightly contradictory if it's a .gpkg. I'll include zoning_class, area_m2, and geometry.
    
    out_gdf = out_gdf[['zoning_class', 'area_m2', 'geometry']]
    
    print(f"Saving to zoning_aggregated.gpkg...")
    out_gdf.to_file('zoning_aggregated.gpkg', driver='GPKG')
    print("Done!")

if __name__ == "__main__":
    solve()

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) or str(val).strip() == "":
        return None
    
    val = str(val).strip().lower()
    
    # Mapping variations to canonical TitleCase values
    mapping = {
        'residential': 'Residential',
        'resi.': 'Residential',
        'res': 'Residential',
        'commercial': 'Commercial',
        'comm.': 'Commercial',
        'comm': 'Commercial',
        'industrial': 'Industrial',
        'indus.': 'Industrial',
        'indus': 'Industrial',
        'agricultural': 'Agricultural',
        'agri.': 'Agricultural',
        'agri': 'Agricultural'
    }
    
    # Check for exact matches or prefix matches
    for key, canonical in mapping.items():
        if val == key or val.startswith(key + '.') or val.startswith(key + ' '):
            return canonical
            
    # Fallback: try to match by start of word to handle things like 'Resi.'
    if val.startswith('resi'): return 'Residential'
    if val.startswith('comm'): return 'Commercial'
    if val.startswith('indu'): return 'Industrial'
    if val.startswith('agri'): return 'Agricultural'
    
    return val.capitalize()

def solve():
    print("Loading data...")
    gdf = gpd.read_file('lagos_zoning_legacy.gpkg')
    
    print("Cleaning zoning classes...")
    gdf['zoning_class'] = gdf['zoning_class'].apply(clean_zoning_class)
    
    # Drop rows with None/NaN/Empty zoning_class
    gdf = gdf.dropna(subset=['zoning_class'])
    gdf = gdf[gdf['zoning_class'].str.strip() != ""]
    
    print(f"Unique classes after cleaning: {gdf['zoning_class'].unique()}")
    
    # Fix vertex precision by snapping to a grid. 
    # Using a slightly larger grid to ensure tiny gaps/overlaps are closed.
    try:
        from shapely import set_precision
        print("Using shapely.set_precision...")
        gdf['geometry'] = gdf['geometry'].apply(lambda g: set_precision(g, grid_size=0.05))
    except ImportError:
        print("shapely.set_precision not available, using coordinate rounding.")
        def round_coords(geom):
            from shapely.ops import transform
            return transform(lambda x, y: (round(x, 2), round(y, 2)), geom)
        gdf['geometry'] = gdf['geometry'].apply(round_coords)

    print("Dissolving geometries by zoning class...")
    # Ensure geometry is valid before dissolving
    gdf['geometry'] = gdf['geometry'].buffer(0)
    
    aggregated_data = []
    for z_class, group in gdf.groupby('zoning_class'):
        # unary_union dissolves boundaries
        merged_geom = group['geometry'].unary_union
        area = merged_geom.area
        aggregated_data.append({
            'zoning_class': z_class,
            'area_m2': area,
            'geometry': merged_geom
        })
        
    print("Creating output GeoDataFrame...")
    out_gdf = gpd.GeoDataFrame(aggregated_data, crs=gdf.crs)
    
    # The instruction asks for two columns (zoning_class and area_m2) 
    # but since it's a .gpkg, the geometry column must be present.
    # We'll include zoning_class, area_m2 and the geometry.
    out_gdf = out_gdf[['zoning_class', 'area_m2', 'geometry']]
    
    print(f"Saving to zoning_aggregated.gpkg...")
    out_gdf.to_file('zoning_aggregated.gpkg', driver='GPKG')
    print("Done!")

if __name__ == "__main__":
    solve()