spa-l2-lagos-hotspot-overlaps

Rerun

import geopandas as gpd

landuse = gpd.read_file('lagos_landuse.geojson')
hex_grid = gpd.read_file('lagos_hex_grid.geojson')

print("Landuse columns:", landuse.columns)
print("Landuse head:\n", landuse.head())
print("Landuse geometry type:", landuse.geometry.type.unique())
print("Landuse CRS:", landuse.crs)

print("\nHex grid columns:", hex_grid.columns)
print("Hex grid head:\n", hex_grid.head())
print("Hex grid geometry type:", hex_grid.geometry.type.unique())
print("Hex grid CRS:", hex_grid.crs)

print("\nLanduse area stats:")
print(landuse.geometry.area.describe())

import geopandas as gpd
import pandas as pd

def solve():
    # 1. Load data
    landuse = gpd.read_file('lagos_landuse.geojson')
    hex_grid = gpd.read_file('lagos_hex_grid.geojson')

    # 2. Setup CRS. Need a projected CRS for area calculations.
    # Nigeria's national grid is typically UTM Zone 31N (EPSG:32631) or similar.
    # Let's use EPSG:32631 for Lagos.
    target_crs = 'EPSG:32631'
    
    landuse = landuse.to_crs(target_crs)
    hex_grid = hex_grid.to_crs(target_crs)

    # 3. Pre-processing: filter tiny polygons (< 100 m2)
    # We need to track how many were filtered.
    landuse['orig_area'] = landuse.geometry.area
    landuse['is_sliver'] = landuse['orig_area'] < 100
    
    # Number of slivers to be filtered (per hex) is tricky. 
    # The requirement says "n_slivers_filtered" in the parquet.
    # This implies we count how many slivers were present in each hex before filtering them out.

    # 4. Spatial Join: hex_grid <-> landuse
    # We want to find the intersection of each hex and each landuse polygon.
    # overlay is better for calculating intersection areas.
    intersections = gpd.overlay(hex_grid, landuse, how='intersection')
    
    # The intersection geometry is the part of the landuse polygon inside the hex.
    intersections['inter_area'] = intersections.geometry.area

    # 5. Calculate weighted density
    # area_weighted_density = sum(pop_density_i * area_i) / sum(area_i)
    # where area_i is the intersection area.
    
    # We need to track:
    # - n_overlap_polygons: count of landuse polygons overlapping the hex (excluding slivers?) 
    #   Actually, the instruction says "n_overlap_polygons (integer), and n_slivers_filtered (integer)".
    #   It doesn't explicitly say if n_overlap_polygons should include slivers. 
    #   Let's assume n_overlap_polygons is the count of ALL landuse polygons that intersect the hex.
    #   And n_slivers_filtered is the count of polygons < 100m2 that intersect the hex.
    
    # To get these counts correctly, let's first do a simple spatial join to get counts.
    # But wait, overlay might be slow. Let's check the counts.
    
    # First, get sliver status for the landuse
    # We already have landuse['is_sliver']
    
    # For each hex, we want to calculate the density based on non-sliver intersections.
    # But we also need to count the slivers that *would* have been there.
    
    # Let's use the intersections result.
    # Note: intersections contains 'is_sliver' from the landuse.
    
    # Filter out slivers for the density calculation
    valid_intersections = intersections[~intersections['is_sliver']].copy()
    
    # Group by hex_id to calculate weighted density
    # Grouping by hex_id to get the sum of (pop_density * inter_area) and sum of inter_area
    grouped = valid_intersections.groupby('hex_id').agg(
        sum_weighted_density=('pop_density', lambda x: (x * valid_intersections.loc[x.index, 'inter_area']).sum()),
        total_inter_area=('inter_area', 'sum')
    )
    
    # Wait, the lambda above is slightly wrong in usage with groupby. Let's fix it.
    valid_intersections['weighted_val'] = valid_intersections['pop_density'] * valid_intersections['inter_area']
    grouped = valid_intersections.groupby('hex_id').agg(
        sum_weighted_val=('weighted_val', 'sum'),
        total_inter_area=('inter_area', 'sum')
    )
    
    grouped['area_weighted_density'] = grouped['sum_weighted_val'] / grouped['total_inter_area']

    # 6. Get counts
    # n_overlap_polygons: total number of landuse polygons intersecting the hex
    # n_slivers_filtered: total number of sliver polygons intersecting the hex
    
    # We can use the intersections dataframe (before filtering slivers) for counts.
    counts = intersections.groupby('hex_id').agg(
        n_overlap_polygons=('id', 'count'),
        n_slivers_filtered=('is_sliver', 'sum')
    )
    
    # Merge density and counts
    results = grouped.join(counts, how='left').fillna(0)
    results['n_overlap_polygons'] = results['n_overlap_polygons'].astype(int)
    results['n_slivers_filtered'] = results['n_slivers_filtered'].astype(int)
    
    # 7. Select top 10%
    # The requirement: "top 10 % of hex cells ranked by area-weighted mean population density"
    # This implies we take the top 10% of the cells that actually have a calculated density.
    # If a hex has no valid landuse, its density is effectively 0 or NaN.
    
    # Sort by density descending
    results = results.sort_values('area_weighted_density', ascending=False)
    
    # Calculate top 10% threshold
    # "top 10% of hex cells" - usually refers to the total number of cells in the hex grid,
    # but ranking is based on density. If we only rank cells with density, we might miss some.
    # However, "ranked by area-weighted mean population density" suggests we look at the ones with data.
    # Let's take the top 10% of the hex_grid.
    
    num_top_cells = max(1, int(len(hex_grid) * 0.1))
    top_cells_ids = results.index[:num_top_cells]
    
    top_results = results.loc[top_cells_ids].copy()
    
    # Assign Rank
    # "rank must be a unique integer starting at 1 where ascending rank corresponds to descending area_weighted_density"
    top_results['rank'] = range(1, len(top_results) + 1)
    
    # 8. Finalize Outputs
    # hotspot_ranking.parquet: hex_id, rank, area_weighted_density, n_overlap_polygons, n_slivers_filtered
    ranking_df = top_results.reset_index()[['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered']]
    ranking_df.to_parquet('hotspot_ranking.parquet')
    
    # hotspots.geoparquet: hex_id, rank, area_weighted_density, geometry
    # We need to join the geometry from hex_grid
    hotspots_gdf = hex_grid.merge(top_results[['rank', 'area_weighted_density']], left_on='hex_id', right_index=True)
    # Keep only required columns
    hotspots_gdf = hotspots_gdf[['hex_id', 'rank', 'area_weighted_density', 'geometry']]
    
    # Convert back to the original CRS (or national grid)? 
    # Instruction says "Nigeria's national grid". 
    # If I used 32631, that is a valid projection for Nigeria. 
    # Let's stick with the projection used for calculation to be safe, 
    # unless I should re-project to something else. 
    # Usually, "national grid" implies a specific CRS like EPSG:26331 or similar. 
    # But EPSG:32631 is standard. Let's assume 32631 is fine.
    # Actually, the input was 4326. Let's check if there's a specific Nigerian CRS.
    # For the purpose of this task, if no specific EPSG is given, 32631 is a good choice for Lagos.
    
    hotspots_gdf.to_geoparquet('hotspots.geoparquet')
    
    print(f"Successfully processed {len(top_results)} hotspots.")
    print(f"Total hexes: {len(hex_grid)}")

if __name__ == "__main__":
    solve()

import geopandas as gpd
import pandas as pd

def solve():
    # 1. Load data
    landuse = gpd.read_file('lagos_landuse.geojson')
    hex_grid = gpd.read_file('lagos_hex_grid.geojson')

    # 2. Setup CRS. Need a projected CRS for area calculations.
    # Nigeria's national grid is typically UTM Zone 31N (EPSG:32631) or similar.
    target_crs = 'EPSG:32631'
    
    landuse = landuse.to_crs(target_crs)
    hex_grid = hex_grid.to_crs(target_crs)

    # 3. Pre-processing: filter tiny polygons (< 100 m2)
    landuse['orig_area'] = landuse.geometry.area
    landuse['is_sliver'] = landuse['orig_area'] < 100

    # 4. Spatial Join: hex_grid <-> landuse
    # overlay is better for calculating intersection areas.
    intersections = gpd.overlay(hex_grid, landuse, how='intersection')
    intersections['inter_area'] = intersections.geometry.area

    # 5. Calculate weighted density
    # area_weighted_density = sum(pop_density_i * area_i) / sum(area_i)
    # where area_i is the intersection area.
    
    # Filter out slivers for the density calculation
    valid_intersections = intersections[~intersections['is_sliver']].copy()
    valid_intersections['weighted_val'] = valid_intersections['pop_density'] * valid_intersections['inter_area']
    
    # Grouping to get density
    density_stats = valid_intersections.groupby('hex_id').agg(
        sum_weighted_val=('weighted_val', 'sum'),
        total_inter_area=('inter_area', 'sum')
    )
    density_stats['area_weighted_density'] = density_stats['sum_weighted_val'] / density_stats['total_inter_area']

    # 6. Get counts
    # n_overlap_polygons: total number of landuse polygons intersecting the hex
    # n_slivers_filtered: total number of sliver polygons intersecting the hex
    counts = intersections.groupby('hex_id').agg(
        n_overlap_polygons=('id', 'count'),
        n_slivers_filtered=('is_sliver', 'sum')
    )
    
    # Merge density and counts
    # We want to ensure all hexes are considered, even if density is 0
    full_stats = density_stats[['area_weighted_density']].join(counts, how='outer').fillna(0)
    full_stats['n_overlap_polygons'] = full_stats['n_overlap_polygons'].astype(int)
    full_stats['n_slivers_filtered'] = full_stats['n_slivers_filtered'].astype(int)
    
    # 7. Select top 10%
    # Sort by density descending
    full_stats = full_stats.sort_values('area_weighted_density', ascending=False)
    
    # Assign Rank
    full_stats['rank'] = range(1, len(full_stats) + 1)
    
    # The instruction: "top 10 % of hex cells ranked by area-weighted mean population density"
    num_top_cells = max(1, int(len(hex_grid) * 0.1))
    top_results = full_stats.head(num_top_cells).copy()
    
    # 8. Finalize Outputs
    # hotspot_ranking.parquet: hex_id, rank, area_weighted_density, n_overlap_polygons, n_slivers_filtered
    ranking_df = top_results.reset_index()[['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered']]
    ranking_df.to_parquet('hotspot_ranking.parquet')
    
    # hotspots.geoparquet: hex_id, rank, area_weighted_density, geometry
    hotspots_gdf = hex_grid.merge(top_results[['rank', 'area_weighted_density']], on='hex_id')
    hotspots_gdf = hotspots_gdf[['hex_id', 'rank', 'area_weighted_density', 'geometry']]
    hotspots_gdf.to_parquet('hotspots.geoparquet', index=False)
    
    print(f"Successfully processed {len(top_results)} hotspots.")
    print(f"Total hexes: {len(hex_grid)}")

if __name__ == "__main__":
    solve()