dc-l2-lagos-snap-normalize

Rerun

import geopandas as gpd

# Load the data to inspect columns and values
gdf = gpd.read_file('lagos_zoning_legacy.gpkg')

print("Columns:", gdf.columns.tolist())
print("\nFirst few rows:")
print(gdf.head())
print("\nUnique zoning_class values before cleaning:")
print(gdf['zoning_class'].unique())
print("\nCRS:", gdf.crs)

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

def clean_zoning_class(val):
    if pd.isna(val):
        return None
    s = str(val).strip().lower()
    if not s:
        return None
    
    # Mapping of variations to canonical TitleCase names
    # Based on the inspection:
    # 'RESIDENTIAL', 'residential', 'Resi.', 'resi.', 'RESI.', 'Residential' -> 'Residential'
    # 'Comm.', 'comm.', 'COMM.', 'Commercial', 'COMMERCIAL', 'commercial' -> 'Commercial'
    # 'Indus.', 'indus.', 'INDUS.', 'Industrial', 'INDUSTRIAL', 'industrial' -> 'Industrial'
    # 'AGRI.', 'Agricultural', 'AGRICULTURAL', 'agricultural', 'Agri.', 'agri.' -> 'Agricultural'
    
    mapping = {
        'resi.': 'Residential',
        'residential': 'Residential',
        'res': 'Residential',
        'comm.': 'Commercial',
        'commercial': 'Commercial',
        'comm': 'Commercial',
        'indus.': 'Industrial',
        'industrial': 'Industrial',
        'indus': 'Industrial',
        'agri.': 'Agricultural',
        'agricultural': 'Agricultural',
        'agri': 'Agricultural'
    }
    
    # Try exact match in mapping, otherwise fallback to Title Case if it's not a known abbr
    # But since we want to be safe, let's look for substrings or common abbreviations
    
    for key, canonical in mapping.items():
        if s == key or s.startswith(key):
            return canonical
            
    # If no mapping matches, try title case as a fallback
    return s.title()

def solve():
    # 1. Load data
    gdf = gpd.read_file('lagos_zoning_legacy.gpkg')
    print(f"Original row count: {len(gdf)}")

    # 2. Clean zoning_class
    # Drop rows with null or whitespace-only classes
    gdf = gdf[gdf['zoning_class'].astype(str).str.strip() != '']
    gdf = gdf.dropna(subset=['zoning_class'])
    
    # Apply cleaning function
    gdf['zoning_class'] = gdf['zoning_class'].apply(clean_zoning_class)
    
    # Drop any that became None after cleaning
    gdf = gdf.dropna(subset=['zoning_class'])
    print(f"Row count after cleaning class: {len(gdf)}")
    print("New unique classes:", gdf['zoning_class'].unique())

    # 3. Fix vertex precision issues
    # A common way to fix tiny gaps/overlaps is to round the coordinates 
    # and then use buffer(0) or union to resolve geometry issues.
    # We'll round to a reasonable precision (e.g., 2 decimal places if metric)
    
    def fix_geometry(geom):
        if geom is None:
            return None
        # Rounding coordinates to fix precision gaps
        # This effectively snaps vertices to a grid
        from shapely.ops import snap
        # But simpler: round coordinates manually or use a small buffer/union
        # Let's try rounding first
        return make_valid(geom)

    # Actually, a better approach for "tiny gaps/overlaps" in GIS:
    # 1. Round coordinates to a specific precision to snap them.
    # 2. Dissolve (union) by zoning class.
    
    def round_coords(geom, precision=2):
        from shapely.geometry import mapping
        import json
        # Using a more robust way to round coords:
        if geom.is_empty:
            return geom
        
        from shapely.wkt import dumps, loads
        # This is a bit hacky but effective for rounding
        # Alternatively, use a custom function to walk the tree
        # Let's use the buffer(0) trick after rounding
        return geom

    # We will use a more standard approach:
    # 1. Ensure valid geometries
    gdf['geometry'] = gdf['geometry'].apply(lambda x: make_valid(x))
    
    # 2. To handle gaps/overlaps, we will dissolve. 
    # To snap vertices, we'll round the coordinates of the geometries.
    def snap_to_grid(geom, precision=0.01):
        from shapely.geometry import MultiPolygon, Polygon, LineString, Point
        
        def _round_coords(coords):
            return tuple(round(c, 2) for c in coords) # 2 decimal places is ~1cm in most UTMs

        if geom.is_empty:
            return geom
        
        # We can use the simple approach: buffer(0) often fixes small issues
        # but doesn't snap. Snapping is better.
        # Let's just round the coords manually via a trick
        import numpy as np
        # This is complex to do for all geom types manually, 
        # so we'll use a simplified approach: buffer(0) + dissolution
        return geom.buffer(0)

    gdf['geometry'] = gdf['geometry'].apply(lambda x: make_valid(x).buffer(0))

    # 4. Merge (Dissolve) by zoning_class
    # First, dissolve to get one geometry per class
    # We also need to calculate the total area
    
    # Dissolve and aggregate area
    # Note: Dissolving can be slow on large datasets.
    # We'll group by zoning_class and union the geometries.
    
    dissolved = gdf.dissolve(by='zoning_class', aggfunc={'area_m2': 'sum'})
    
    # After dissolving, the geometries might still have tiny gaps if they weren't perfectly snapped.
    # However, the prompt asks to "fix the vertex precision problems" AND "produce one merged geometry".
    # If we dissolve, the gaps between polygons of the same class will disappear.
    # The only thing left are gaps between DIFFERENT classes (which aren't part of the same dissolved polygon)
    # OR gaps/overlaps within the same class that weren't closed.
    
    # To truly fix vertex precision for the whole layer, we should snap all vertices to a grid first.
    # Let's refine the snap process.
    
    print("Proceeding with dissolution...")
    
    # To address the "vertex precision" specifically before dissolving:
    # A common trick is to round coordinates.
    from shapely import affinity
    
    def precision_fix(geom, decimals=2):
        # This is a bit heavy, but let's try to round the coordinates
        # A faster way: use a small buffer
        return geom.buffer(0.001).buffer(-0.001) # This is a "cleaning" trick but can shrink polygons
    
    # Actually, let's just use the dissolve. Dissolve(by='zoning_class') 
    # is the standard way to "produce one merged geometry per canonical zoning class".
    
    # Final check on columns
    # Result should have: zoning_class (index after dissolve), area_m2
    # We need to reset index to make zoning_class a column
    result = dissolved.reset_index()
    
    # Reorder/Select columns
    result = result[['zoning_class', 'area_m2', 'geometry']]
    
    # Calculate area_m2 again if the original area_m2 was just for parcels
    # The user wants the area of the merged geometry.
    # Let's recalculate area from the merged geometry to be precise.
    result['area_m2'] = result['geometry'].area
    
    # Final output selection
    result = result[['zoning_class', 'area_m2', 'geometry']]

    # Write to file
    result.to_file('zoning_aggregated.gpkg', driver="GPKG")
    print("Successfully wrote zoning_aggregated.gpkg")
    print(result)

if __name__ == "__main__":
    solve()

import geopandas as gpd
import pandas as pd
from shapely.validation import make_valid
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
    
    # Mapping of variations to canonical TitleCase names
    mapping = {
        'resi.': 'Residential',
        'residential': 'Residential',
        'res': 'Residential',
        'comm.': 'Commercial',
        'commercial': 'Commercial',
        'comm': 'Commercial',
        'indus.': 'Industrial',
        'industrial': 'Industrial',
        'indus': 'Industrial',
        'agri.': 'Agricultural',
        'agricultural': 'Agricultural',
        'agri': 'Agricultural'
    }
    
    # Check for exact match in mapping
    if s in mapping:
        return mapping[s]
        
    # Fallback for things like "RESIDENTIAL" (if not caught by mapping)
    # or other variations that might just need Title Case.
    # We check if the string is basically a known prefix
    if s.startswith('resi'): return 'Residential'
    if s.startswith('comm'): return 'Commercial'
    if s.startswith('indus'): return 'Industrial'
    if s.startswith('agri'): return 'Agricultural'

    return s.title()

def solve():
    # 1. Load data
    gdf = gpd.read_file('lagos_zoning_legacy.gpkg')
    print(f"Original row count: {len(gdf)}")

    # 2. Clean zoning_class
    # Drop rows with null or whitespace-only classes
    gdf = gdf[gdf['zoning_class'].astype(str).str.strip() != '']
    gdf = gdf.dropna(subset=['zoning_class'])
    
    # Apply cleaning function
    gdf['zoning_class'] = gdf['zoning_class'].apply(clean_zoning_class)
    
    # Drop any that became None after cleaning
    gdf = gdf.dropna(subset=['zoning_class'])
    print(f"Row count after cleaning class: {len(gdf)}")
    print("New unique classes:", gdf['zoning_class'].unique())

    # 3. Fix vertex precision problems
    # To fix tiny gaps/overlaps, we round the coordinates to a specific precision.
    # This 'snaps' vertices to a grid.
    def round_coords(geom, decimals=2):
        if geom.is_empty:
            return geom
        # Use a transformation to round coordinates
        # A simple way to do this is to use a small buffer or manual coordinate rounding
        # But for many geometries, let's use the coordinate rounding approach
        def _round_point(x, y):
            return (round(x, decimals), round(y, decimals))

        if geom.geom_type == 'Polygon':
            return Polygon([_round_point(x, y) for x, y in geom.exterior.coords]) + \
                   MultiPolygon([Polygon([_round_point(x, y) for x, y in poly.exterior.coords]) 
                                 for poly in geom.interiors]) # This is wrong for interiors
        # Using a better way: just use a tiny buffer to close gaps before dissolving
        return geom

    # Let's try a simpler approach for "vertex precision":
    # 1. Make valid
    # 2. Buffer(0) to fix topology
    # 3. Use a small buffer to merge touching edges if they have gaps
    # Actually, the most standard way to handle this in GIS for 'cleaning' is to 
    # snap the layer to a grid. Let's use a simple rounding function.
    
    from shapely.ops import transform

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

    print("Snapping to grid...")
    gdf['geometry'] = gdf['geometry'].apply(lambda x: snap_to_grid(make_valid(x), 0.01))

    # 4. Merge (Dissolve) by zoning_class
    print("Dissolving by zoning_class...")
    # We dissolve to merge all geometries of the same class into one
    dissolved = gdf.dissolve(by='zoning_class')
    
    # After dissolve, ensure we have one geometry per class and calculate area
    # Dissolve by default takes the first 'area_m2' or can be aggregated.
    # The user wants the area of the merged geometry.
    
    # Reset index to make 'zoning_class' a column again
    result = dissolved.reset_index()
    
    # Calculate the true area of the merged geometry
    result['area_m2'] = result['geometry'].area
    
    # Final column selection: zoning_class, area_m2, and geometry
    # (geometry is included as it's a spatial file)
    result = result[['zoning_class', 'area_m2', 'geometry']]

    # Write to file
    result.to_file('zoning_aggregated.gpkg', driver="GPKG")
    print("Successfully wrote zoning_aggregated.gpkg")
    print(result)

if __name__ == "__main__":
    solve()

import geopandas as gpd
import pandas as pd
from shapely.validation import make_valid
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 = {
        'resi.': 'Residential',
        'residential': 'Residential',
        'res': 'Residential',
        'comm.': 'Commercial',
        'commercial': 'Commercial',
        'comm': 'Commercial',
        'indus.': 'Industrial',
        'industrial': 'Industrial',
        'indus': 'Industrial',
        'agri.': 'Agricultural',
        'agricultural': 'Agricultural',
        'agri': 'Agricultural'
    }
    
    if s in mapping:
        return mapping[s]
    
    if s.startswith('resi'): return 'Residential'
    if s.startswith('comm'): return 'Commercial'
    if s.startswith('indus'): return 'Industrial'
    if s.startswith('agri'): return 'Agricultural'

    return s.title()

def solve():
    # 1. Load data
    gdf = gpd.read_file('lagos_zoning_legacy.gpkg')
    print(f"Original row count: {len(gdf)}")

    # 2. Clean zoning_class
    gdf = gdf[gdf['zoning_class'].astype(str).str.strip() != '']
    gdf = gdf.dropna(subset=['zoning_class'])
    gdf['zoning_class'] = gdf['zoning_class'].apply(clean_zoning_class)
    gdf = gdf.dropna(subset=['zoning_class'])
    print(f"Row count after cleaning class: {len(gdf)}")

    # 3. Fix vertex precision problems
    # Snapping to a 0.01 grid to resolve gaps/overlaps
    def snap_to_grid(geom, precision=0.01):
        def _snap(x, y):
            return (round(x / precision) * precision, round(y / precision) * precision)
        return transform(_snap, geom)

    print("Fixing geometry precision...")
    # Make valid, then snap to grid, then buffer(0) to clean up any topology issues introduced by snapping
    gdf['geometry'] = gdf['geometry'].apply(lambda x: make_valid(x))
    gdf['geometry'] = gdf['geometry'].apply(lambda x: snap_to_grid(x, 0.01))
    gdf['geometry'] = gdf['geometry'].apply(lambda x: x.buffer(0))

    # 4. Merge (Dissolve) by zoning_class
    print("Dissolving by zoning_class...")
    dissolved = gdf.dissolve(by='zoning_class')
    
    result = dissolved.reset_index()
    # Recalculate area from the dissolved geometry
    result['area_m2'] = result['geometry'].area
    
    # Final column selection: zoning_class, area_m2, and geometry
    result = result[['zoning_class', 'area_m2', 'geometry']]

    # Write to file
    result.to_file('zoning_aggregated.gpkg', driver="GPKG")
    print("Successfully wrote zoning_aggregated.gpkg")
    print(result)

if __name__ == "__main__":
    solve()