dc-l2-lagos-snap-normalize

Rerun

import geopandas as gpd

# Load the data
gdf = gpd.read_file('lagos_zoning_legacy.gpkg')

print("Columns:", gdf.columns.tolist())
print("First few rows:")
print(gdf.head())
print("\nValue counts for zoning_class:")
print(gdf['zoning_class'].value_counts(dropna=False))
print("\nCRS:", gdf.crs)

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

def clean_zoning_class(val):
    if pd.isna(val):
        return None
    s = str(val).strip().lower()
    if not s:
        return None
    
    # Mapping dictionary for variations
    mapping = {
        'residential': 'Residential',
        'resi.': 'Residential',
        'resi': 'Residential',
        'commercial': 'Commercial',
        'comm.': 'Commercial',
        'comm': 'Commercial',
        'industrial': 'Industrial',
        'indus.': 'Industrial',
        'indus': 'Industrial',
        'agricultural': 'Agricultural',
        'agri.': 'Agricultural',
        'agri': 'Agricultural'
    }
    
    # Try direct mapping first
    if s in mapping:
        return mapping[s]
    
    # Fallback to Title Case if not in mapping, but let's see if we missed something
    return s.title()

def solve():
    # 1. Load data
    print("Loading data...")
    gdf = gpd.read_file('lagos_zoning_legacy.gpkg')

    # 2. Clean zoning_class
    print("Cleaning zoning_class...")
    gdf['zoning_class'] = gdf['zoning_class'].apply(clean_zoning_class)
    
    # Drop rows where zoning_class is None or whitespace-only (already handled by None in clean_zoning_class)
    gdf = gdf.dropna(subset=['zoning_class'])
    gdf = gdf[gdf['zoning_class'].str.strip() != '']

    # 3. Fix vertex precision issues (tiny gaps/overlaps)
    # A common way is to round coordinates to a certain precision or use a snap-to-grid approach.
    # Since it's a metric CRS (EPSG:26331), 1mm is 0.001. Let's round to 1cm (0.01) or 1mm (0.001).
    # Let's try rounding to 1mm (0.001) to avoid losing too much detail but snap vertices.
    print("Fixing vertex precision...")
    def snap_geometry(geom, precision=0.001):
        if geom is None:
            return None
        # Rounding coordinates is a simple way to snap to a grid
        from shapely.ops import transform
        def round_coords(x, y, z=None):
            return (round(x / precision) * precision, round(y / precision) * precision)
        return transform(lambda x, y, z=None: round_coords(x, y), geom)

    gdf['geometry'] = gdf['geometry'].apply(lambda g: snap_geometry(g))

    # 4. Merge geometries by zoning_class
    print("Merging geometries by zoning_class...")
    # We need to dissolve the geometries based on zoning_class.
    # To ensure we handle overlaps/gaps properly after snapping, 
    # we should ideally use a buffer(0) or union. 
    # Dissolve in geopandas is a union operation.
    
    # Group by zoning_class and dissolve
    dissolved = gdf.dissolve(by='zoning_class', aggfunc={'area_m2': 'sum'})
    
    # The dissolve operation results in 'zoning_class' being the index.
    # We need to move it back to a column and ensure 'area_m2' is preserved.
    dissolved = dissolved.reset_index()
    
    # Note: The area_m2 in the prompt might mean the area of the merged polygon.
    # "each row carrying two columns: zoning_class and a numeric area_m2"
    # If they want the total area of the resulting polygon, we should recalculate it.
    print("Recalculating merged area...")
    dissolved['area_m2'] = dissolved['geometry'].area

    # 5. Save to file
    print(f"Saving to zoning_aggregated.gpkg...")
    # Ensure columns are in correct order/subset
    output_gdf = dissolved[['zoning_class', 'area_m2', 'geometry']]
    output_gdf.to_file('zoning_aggregated.gpkg', driver='GPKG')
    
    print("Done!")
    print(output_gdf.head())

if __name__ == "__main__":
    solve()