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("\nFirst 5 rows:")
    print(gdf.head())
    print("\nValue counts for 'zoning_class':")
    if 'zoning_class' in gdf.columns:
        print(gdf['zoning_class'].value_counts(dropna=False))
    else:
        print("Column 'zoning_class' not found!")
    print("\nGeometry type:", gdf.geom_type.unique())
    print("\nCRS:", gdf.crs)

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

import geopandas as gpd

def inspect_data(file_path):
    gdf = gpd.read_file(file_path)
    print("Columns:", gdf.columns.tolist())
    print("\nFirst 5 rows:")
    print(gdf.head())
    print("\nValue counts for 'zoning_class':")
    if 'zoning_class' in gdf.columns:
        print(gdf['zoning_class'].value_counts(dropna=False))
    else:
        print("Column 'zoning_class' not found!")
    print("\nGeometry type:", gdf.geom_type.unique())
    print("\nCRS:", gdf.crs)

if __name__ == "__main__":
    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) or str(val).strip() == "":
        return None
    
    s = str(val).strip().lower()
    
    # Mapping variations to canonical TitleCase names
    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 to find a match in mapping, otherwise use TitleCase version of the cleaned string
    return mapping.get(s, s.capitalize())

def solve():
    input_file = '/work/lagos_zoning_legacy.gpkg'
    output_file = '/work/zoning_aggregated.gpkg'
    
    print(f"Reading {input_file}...")
    gdf = gpd.read_file(input_file)
    
    # 1. Clean zoning_class
    print("Cleaning zoning classes...")
    gdf['zoning_class'] = gdf['zoning_class'].apply(clean_zoning_class)
    
    # 2. Drop rows with blank/None classes
    initial_count = len(gdf)
    gdf = gdf.dropna(subset=['zoning_class'])
    gdf = gdf[gdf['zoning_class'].str.strip() != ""]
    print(f"Dropped {initial_count - len(gdf)} rows with invalid zoning classes.")

    # 3. Fix vertex precision problems (tiny gaps/overlaps)
    # We use buffer(0) and snap to a grid to fix topology issues.
    # A small buffer can help close gaps, but buffer(0) is often enough for cleaning.
    # To be more robust, let's use a small positive buffer then a negative one, 
    # or simply snap geometries to a small precision grid.
    print("Fixing geometry precision...")
    # Snapping to a grid of 0.01 (assuming units are meters, this is 1cm)
    gdf['geometry'] = gdf['geometry'].buffer(0.001).buffer(-0.001) # Slight expansion/contraction to snap
    # Or more simply:
    gdf['geometry'] = gdf['geometry'].apply(lambda geom: geom.simplify(0.01, preserve_topology=True))
    
    # 4. Produce one merged geometry per canonical zoning class
    print("Aggregating geometries by zoning class...")
    
    results = []
    for z_class, group in gdf.groupby('zoning_class'):
        # unary_union dissolves all boundaries between polygons in the group
        merged_geom = group.geometry.unary_union
        area = merged_geom.area
        results.append({
            'zoning_class': z_class,
            'area_m2': area,
            'geometry': merged_geom
        })
        
    agg_gdf = gpd.GeoDataFrame(results, crs=gdf.crs)
    
    # Ensure the geometry type is consistent (MultiPolygon if there are separate islands)
    # unary_union often returns Polygon or MultiPolygon.
    
    print(f"Writing results to {output_file}...")
    agg_gdf.to_file(output_file, 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
    
    s = str(val).strip().lower()
    
    # Mapping variations to canonical TitleCase names
    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 to find a match in mapping, otherwise use TitleCase version of the cleaned string
    return mapping.get(s, s.capitalize())

def solve():
    input_file = '/work/lagos_zoning_legacy.gpkg'
    output_file = '/work/zoning_aggregated.gpkg'
    
    print(f"Reading {input_file}...")
    gdf = gpd.read_file(input_file)
    
    # 1. Clean zoning_class
    print("Cleaning zoning classes...")
    gdf['zoning_class'] = gdf['zoning_class'].apply(clean_zoning_class)
    
    # 2. Drop rows with blank/None classes
    initial_count = len(gdf)
    gdf = gdf.dropna(subset=['zoning_class'])
    gdf = gdf[gdf['zoning_class'].str.strip() != ""]
    print(f"Dropped {initial_count - len(gdf)} rows with invalid zoning classes.")

    # 3. Fix vertex precision problems (tiny gaps/overlaps)
    # A robust way to fix small gaps/overlaps in a tiling dataset is to use a small buffer(0)
    # or to snap to a grid. 
    # To specifically address "tiny gaps and overlaps", we can use a tiny buffer(dist) 
    # followed by buffer(-dist) to align edges.
    print("Fixing geometry precision...")
    # We'll try a very small buffer to snap vertices.
    # Since the CRS is EPSG:26331 (meters), 0.01 is 1cm.
    gdf['geometry'] = gdf['geometry'].buffer(0.01).buffer(-0.01)
    
    # 4. Produce one merged geometry per canonical zoning class
    print("Aggregating geometries by zoning class...")
    
    results = []
    for z_class, group in gdf.groupby('zoning_class'):
        # Dissolve boundaries between polygons in the group
        merged_geom = group.geometry.unary_union
        area = merged_geom.area
        results.append({
            'zoning_class': z_class,
            'area_m2': float(area),
            'geometry': merged_geom
        })
        
    agg_gdf = gpd.GeoDataFrame(results, crs=gdf.crs)
    
    # Ensure we only keep the requested columns
    agg_gdf = agg_gdf[['zoning_class', 'area_m2', 'geometry']]
    
    print(f"Writing results to {output_file}...")
    agg_gdf.to_file(output_file, driver="GPKG")
    print("Done!")
    print("Resulting classes:")
    print(agg_gdf[['zoning_class', 'area_m2']])

if __name__ == "__main__":
    solve()

import geopandas as gpd
import pandas as pd

def inspect_output(file_path):
    gdf = gpd.read_file(file_path)
    print("Columns:", gdf.columns.tolist())
    print("\nContent:")
    print(gdf)
    print("\nGeometry type counts:")
    print(gdf.geom_type.value_counts())

if __name__ == "__main__":
    inspect_output('/work/zoning_aggregated.gpkg')