dc-l2-lagos-snap-normalize

Rerun

import geopandas as gpd

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

# Inspect the data
print("Initial columns:", gdf.columns.tolist())
print("Initial head:")
print(gdf.head())
print("\nValue counts for zoning_class:")
print(gdf['zoning_class'].value_counts(dropna=False))
print("\nGeometry types:", gdf.geometry.type.unique())
print("\nCRS:", gdf.crs)

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

def clean_zoning_class(val):
    if pd.isna(val):
        return None
    s = str(val).strip().lower()
    if not s:
        return None
    
    # Mapping dictionary for spelling variations
    # Based on inspection: 
    # Residential: 'residential', 'resi.'
    # Commercial: 'commercial', 'comm.'
    # Industrial: 'industrial', 'indus.'
    # Agricultural: 'agricultural', 'agri.'
    
    mapping = {
        'residential': 'Residential',
        'resi.': 'Residential',
        'commercial': 'Commercial',
        'comm.': 'Commercial',
        'industrial': 'Industrial',
        'indus.': 'Industrial',
        'agricultural': 'Agricultural',
        'agri.': 'Agricultural'
    }
    
    return mapping.get(s, s.title())

def solve():
    # 1. Load
    gdf = gpd.read_file('lagos_zoning_legacy.gpkg')
    
    # 2. Clean zoning_class
    gdf['zoning_class'] = gdf['zoning_class'].apply(clean_zoning_class)
    
    # 3. Drop rows with blank/None classes
    gdf = gdf.dropna(subset=['zoning_class'])
    gdf = gdf[gdf['zoning_class'].str.strip() != '']
    
    print("After cleaning classes, value counts:")
    print(gdf['zoning_class'].value_counts())

    # 4. Fix vertex precision problems (gaps and overlaps)
    # A common way to fix topology issues in GIS is to use a snap-to-grid or buffer(0) approach.
    # However, "tiny gaps and overlaps" are often solved by a small buffer(0) or by snapping to a grid.
    # Since we want to merge them anyway, let's snap the coordinates to a certain precision first.
    # EPSG:26331 is in meters. Snapping to 1cm (0.01m) should help.
    
    def snap_geometry(geom, precision=0.01):
        if geom.is_empty:
            return geom
        # Rounding coordinates is a simple way to snap to a grid
        def round_coords(coords):
            return tuple(round(c, 2) for c in coords)
        
        if isinstance(geom, Polygon):
            return Polygon(round_coords(geom.exterior.coords))
        elif isinstance(geom, MultiPolygon):
            return MultiPolygon([Polygon(round_coords(p.exterior.coords)) for p in geom.geoms])
        return geom

    # Actually, a better approach for gaps/overlaps is buffering by a tiny amount or using a grid.
    # Let's try snapping all coordinates to the nearest 0.01m.
    gdf.geometry = gdf.geometry.apply(lambda x: x.buffer(0.001).buffer(-0.001)) 
    # Wait, buffer(0) is more standard for fixing invalid geometries.
    # Let's use a coordinate rounding approach.
    
    from shapely import wkt
    
    def snap_to_grid(geom, grid_size=0.01):
        if geom.is_empty:
            return geom
        # This is a bit complex to do manually for all types, 
        # but for Polygons we can use the fact that we'll dissolve later.
        # Let's use shapely's rounding if available or just buffer(0).
        return geom.buffer(0)

    gdf.geometry = gdf.geometry.buffer(0)
    
    # Actually, the most robust way to handle "gaps and overlaps" for dissolution
    # is to snap geometries to a grid.
    
    def snap_to_grid_geometry(geom, precision=0.01):
        # This approach snaps all vertices to a grid
        from shapely.ops import transform
        import numpy as np
        
        def snap_func(x, y):
            return np.round(x / precision) * precision, np.round(y / precision) * precision
        
        # We'll use a more direct approach:
        # Just use the built-in precision if possible.
        # Since we are using shapely, we can use the `wkt` or coordinate manipulation.
        pass

    # Let's try the simplest way to fix precision: round the coordinates
    from shapely.affinity import scale
    
    def round_coords_to_precision(geom, precision=0.01):
        # Shapely doesn't have a direct "round" function for all geometries easily
        # but we can use the 'apply' on coordinates.
        # A trick is to use buffer(0) which often cleans up topology.
        return geom.buffer(0)

    # For the specific issue of gaps/overlaps, let's try snapping to 1cm grid
    # A simple way to snap:
    from shapely.geometry import mapping
    from shapely.ops import unary_union

    # Instead of complex snapping, let's use the fact that we want to aggregate.
    # We will group by zoning_class, then unary_union the geometries.
    # unary_union is excellent at merging adjacent polygons and resolving overlaps.
    
    results = []
    for z_class, group in gdf.groupby('zoning_class'):
        # Merge all geometries in this group
        merged_geom = group.geometry.unary_union
        # Calculate total area
        total_area = merged_geom.area
        results.append({
            'zoning_class': z_class,
            'area_m2': total_area,
            'geometry': merged_geom
        })
    
    output_gdf = gpd.GeoDataFrame(results, crs=gdf.crs)
    
    # 5. Write to file
    output_gdf.to_file('zoning_aggregated.gpkg', driver="GPKG")
    print("Successfully wrote zoning_aggregated.gpkg")

if __name__ == "__main__":
    solve()

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

def clean_zoning_class(val):
    if pd.isna(val):
        return None
    s = str(val).strip().lower()
    if not s:
        return None
    
    # Mapping dictionary for spelling variations
    mapping = {
        'residential': 'Residential',
        'resi.': 'Residential',
        'commercial': 'Commercial',
        'comm.': 'Commercial',
        'industrial': 'Industrial',
        'indus.': 'Industrial',
        'agricultural': 'Agricultural',
        'agri.': 'Agricultural'
    }
    
    return mapping.get(s, s.title())

def snap_to_grid(geom, precision=0.01):
    def snap_func(x, y):
        return np.round(x / precision) * precision, np.round(y / precision) * precision
    return transform(snap_func, geom)

def solve():
    # 1. Load
    gdf = gpd.read_file('lagos_zoning_legacy.gpkg')
    
    # 2. Clean zoning_class
    gdf['zoning_class'] = gdf['zoning_class'].apply(clean_zoning_class)
    
    # 3. Drop rows with blank/None classes
    gdf = gdf.dropna(subset=['zoning_class'])
    gdf = gdf[gdf['zoning_class'].str.strip() != '']
    
    print("After cleaning classes, value counts:")
    print(gdf['zoning_class'].value_counts())

    # 4. Fix vertex precision problems (gaps and overlaps)
    print("Snapping geometries to 0.01m grid...")
    gdf.geometry = gdf.geometry.apply(lambda x: snap_to_grid(x, 0.01))
    
    # 5. Dissolve by zoning_class
    print("Dissolving geometries by zoning_class...")
    # Dissolve handles the merging of polygons into one per class
    dissolved = gdf.dissolve(by='zoning_class')
    
    # 6. Calculate area of the merged geometry
    dissolved['area_m2'] = dissolved.geometry.area
    
    # Prepare final GeoDataFrame
    output_gdf = dissolved.reset_index()[['zoning_class', 'area_m2', 'geometry']]
    
    # 7. Write to file
    output_gdf.to_file('zoning_aggregated.gpkg', driver="GPKG")
    print("Successfully wrote zoning_aggregated.gpkg")
    print("\nFinal aggregated data:")
    print(output_gdf)

if __name__ == "__main__":
    solve()