spa-l2-lagos-hotspot-overlaps

Rerun

import geopandas as gpd

# Inspect landuse
landuse = gpd.read_file('lagos_landuse.geojson')
print("Landuse columns:", landuse.columns)
print("Landuse head:\n", landuse.head())
print("Landuse geometry type:", landuse.geometry.type.unique())

# Inspect hex grid
hex_grid = gpd.read_file('lagos_hex_grid.geojson')
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())

import geopandas as gpd
import pandas as pd

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

    # 1. Filter tiny polygons < 100 m^2
    # Assuming input is in WGS84 (degrees) based on the print output. 
    # We need to project to a metric system to calculate area accurately.
    # Nigeria's national grid is typically EPSG:26331 (UTM Zone 31N) or similar.
    # Let's check if we can use a generic metric CRS or if we should find a specific one.
    # For area calculation, projecting to a local UTM is best.
    
    # Let's try to find a suitable UTM zone for Lagos (around 6.5 N, 3.4 E) -> UTM 31N (EPSG:32631)
    projected_crs = 'EPSG:32631'
    
    landuse_proj = landuse.to_crs(projected_crs)
    hex_grid_proj = hex_grid.to_crs(projected_crs)

    # Identify slivers (polygons < 100 m2)
    # We need to keep track of which ones were filtered.
    # However, the requirement says "Tiny polygons ... should be excluded".
    # Does this mean they are excluded from the calculation? Yes.
    
    # We need to count how many slivers are filtered per hex. 
    # This implies we need to know which landuse polygons intersect which hexes.
    
    # To count n_slivers_filtered, we need to see which polygons in the original landuse
    # would have fallen into the hex if they weren't slivers.
    
    # Let's first define the mask for valid landuse polygons.
    landuse_proj['area'] = landuse_proj.geometry.area
    is_sliver = landuse_proj['area'] < 100
    n_slivers_total = is_sliver.sum()
    
    # The task says "n_slivers_filtered" is a column in hotspot_ranking.parquet.
    # This likely means the number of sliver polygons that intersected that specific hex.
    
    # Filter landuse
    landuse_valid = landuse_proj[~is_sliver].copy()
    landuse_slivers = landuse_proj[is_sliver].copy()
    
    # 2. Spatial Join: Intersect hexes with landuse
    # We need area-weighted mean population density.
    # Area-weighted mean = sum(pop_density_i * area_i_in_hex) / total_area_of_hex_covered_by_landuse
    # Wait, the requirement is "area-weighted mean population density across overlapping land-use polygons".
    # Usually, this means for each hex, we look at all landuse polygons that overlap it, 
    # find the intersection area, and weight the density by that intersection area.
    
    # We use the valid landuse for the density calculation.
    # But we need to count slivers per hex.
    
    # Step A: Calculate intersection areas for valid landuse
    # We'll perform a spatial join to find intersections.
    # Using overlay is better for getting the actual intersection geometries.
    
    # Note: hex_grid_proj might be large, and landuse_valid might be large.
    # overlay(how='intersection') can be expensive.
    
    # Let's do it more efficiently.
    # For each hex, we want to compute:
    # density = sum(density_i * area(intersection_i)) / sum(area(intersection_i))
    # wait, "area-weighted mean population density across overlapping land-use polygons".
    # This often means: (sum_{i in hex} density_i * area(intersection_i)) / (sum_{i in hex} area(intersection_i))
    # Actually, the sum of intersection areas might not equal the hex area if there are gaps.
    # But the denominator should be the sum of the intersection areas to get a weighted mean of the densities present.
    
    # Let's perform overlay.
    intersections = gpd.overlay(hex_grid_proj[['hex_id', 'geometry']], landuse_valid[['pop_density', 'geometry']], how='intersection')
    
    # Calculate area of each intersection piece
    intersections['inter_area'] = intersections.geometry.area
    
    # Calculate weighted density components
    intersections['weighted_pop'] = intersections['pop_density'] * intersections['inter_area']
    
    # Group by hex_id
    hex_stats = intersections.groupby('hex_id').agg(
        sum_weighted_pop=('weighted_pop', 'sum'),
        sum_inter_area=('inter_area', 'sum'),
        n_overlap_polygons=('pop_density', 'count') # This might be wrong if one polygon creates multiple intersection parts.
        # The requirement "n_overlap_polygons" usually refers to the number of unique landuse polygons.
    ).reset_index()
    
    # Wait, n_overlap_polygons needs to be the count of unique landuse polygons.
    # Let's adjust the intersection to keep track of landuse IDs.
    
    # Re-doing intersection with landuse ID
    landuse_valid_with_id = landuse_proj.copy()
    landuse_valid_with_id['landuse_id'] = landuse_valid_with_id['id']
    
    intersections = gpd.overlay(hex_grid_proj[['hex_id', 'geometry']], landuse_valid_with_id[['landuse_id', 'pop_density', 'geometry']], how='intersection')
    intersections['inter_area'] = intersections.geometry.area
    intersections['weighted_pop'] = intersections['pop_density'] * intersections['inter_area']
    
    hex_stats = intersections.groupby('hex_id').agg(
        sum_weighted_pop=('weighted_pop', 'sum'),
        sum_inter_area=('inter_area', 'sum'),
        n_overlap_polygons=('landuse_id', 'nunique')
    ).reset_index()
    
    hex_stats['area_weighted_density'] = hex_stats['sum_weighted_pop'] / hex_stats['sum_inter_area']
    
    # Step B: Count slivers per hex
    # We need to know how many slivers intersect each hex.
    sliver_intersections = gpd.overlay(hex_grid_proj[['hex_id', 'geometry']], landuse_slivers[['geometry']], how='intersection')
    # For slivers, we just need to count how many original sliver polygons were involved.
    # Actually, overlay(how='intersection') with landuse_slivers will give us the intersection pieces.
    # We want to count how many sliver polygons intersect each hex.
    # A spatial join is better for this.
    
    sliver_counts = gpd.sjoin(hex_grid_proj[['hex_id', 'geometry']], landuse_slivers[['geometry']], how='left', predicate='intersects')
    # The sjoin will duplicate hex_id for every intersection.
    # We need to count unique landuse indices (or just count how many indices are not null).
    # But landuse_slivers doesn't have a unique ID in the snippet above, let's add one.
    landuse_slivers['sli_id'] = range(len(landuse_slivers))
    sliver_counts = gpd.sjoin(hex_grid_proj[['hex_id', 'geometry']], landuse_slivers[['sli_id', 'geometry']], how='left', predicate='intersects')
    
    # Count non-null sli_id
    hex_sliver_counts = sliver_counts.groupby('hex_id')['sli_id'].nunique().reset_index()
    hex_sliver_counts.columns = ['hex_id', 'n_slivers_filtered']
    
    # Step C: Combine everything
    # Start with hex_grid to ensure we have all hexes if needed? 
    # No, the prompt says "top 10 % of hex cells ranked...". 
    # This implies we only care about hexes that have some landuse? 
    # Or do we rank all hexes? "top 10 % of hex cells ranked by area-weighted mean population density".
    # Usually, if a hex has no landuse, its density is 0.
    
    # Let's merge stats back to the hex_grid.
    final_df = hex_grid[['hex_id', 'geometry']].merge(hex_stats, on='hex_id', how='left')
    final_df = final_df.merge(hex_sliver_counts, on='hex_id', how='left')
    
    # Fill NaNs
    final_df['area_weighted_density'] = final_df['area_weighted_density'].fillna(0)
    final_df['n_overlap_polygons'] = final_df['n_overlap_polygons'].fillna(0).astype(int)
    final_df['n_slivers_filtered'] = final_df['n_slivers_filtered'].fillna(0).astype(int)
    
    # Now find the top 10%
    # Rank descending area_weighted_density.
    # "rank must be a unique integer starting at 1 where ascending rank corresponds to descending area_weighted_density"
    # This means rank 1 is the highest density.
    final_df = final_df.sort_values('area_weighted_density', ascending=False)
    final_df['rank'] = range(1, len(final_df) + 1)
    
    # Top 10%
    n_top = max(1, int(len(final_df) * 0.10))
    hotspots = final_df.head(n_top).copy()
    
    # Ensure the rank is from the full sorted list? 
    # "top 10 % of hex cells ranked ... rank must be a unique integer starting at 1"
    # This is ambiguous. Does it mean the rank is relative to all cells, or 1 to N_top?
    # "rank must be a unique integer starting at 1 where ascending rank corresponds to descending area_weighted_density"
    # This usually implies the rank is the position in the sorted list.
    # And since the output files "share exactly the same set of hex_id values", 
    # the rank should be the same in both.
    
    # Let's check the requirement: "rank must be a unique integer starting at 1 where ascending rank corresponds to descending area_weighted_density"
    # This applies to the "top 10 %". 
    # If I have 100 cells and take top 10, do their ranks become 1-10 or their original ranks?
    # "rank must be a unique integer starting at 1" suggests 1, 2, 3...
    # Let's re-rank the hotspots.
    hotspots = hotspots.sort_values('area_weighted_density', ascending=False)
    hotspots['rank'] = range(1, len(hotspots) + 1)

    # Output 1: hotspots.geoparquet
    # columns: hex_id, rank, area_weighted_density, geometry
    hotspots_geo = hotspots[['hex_id', 'rank', 'area_weighted_density', 'geometry']]
    # Re-project to national grid if needed. The prompt says "in Nigeria's national grid".
    # Nigeria's national grid is typically EPSG:26331 (UTM 31N).
    # But wait, the input is in WGS84 (lat/lon). 
    # The prompt says "in Nigeria's national grid". I'll use EPSG:26331.
    # However, let's check if the input hex_grid is already in some grid. 
    # It's in WGS84. 
    # I will project the output to EPSG:26331.
    hotspots_geo = hotspots_geo.to_crs('EPSG:26331')
    hotspots_geo.to_parquet('hotspots.geoparquet')
    
    # Output 2: hotspot_ranking.parquet
    # columns: hex_id, rank, area_weighted_density, n_overlap_polygons, n_slivers_filtered
    hotspots_rank = hotspots[['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered']]
    hotspots_rank.to_parquet('hotspot_ranking.parquet')
    
    print(f"Successfully created files. Hotspots count: {len(hotspots)}")

if __name__ == "__main__":
    solve()

import geopandas as gpd
import pandas as pd

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

    # Project to UTM Zone 31N (EPSG:32631) for area calculations
    # This is a standard for Nigeria/Lagos.
    projected_crs = 'EPSG:32631'
    landuse_proj = landuse.to_crs(projected_crs)
    hex_grid_proj = hex_grid.to_crs(projected_crs)

    # 1. Identify slivers (< 100 m2)
    landuse_proj['area'] = landuse_proj.geometry.area
    is_sliver = landuse_proj['area'] < 100
    landuse_valid = landuse_proj[~is_sliver].copy()
    landuse_slivers = landuse_proj[is_sliver].copy()
    
    # 2. Calculate area-weighted density for each hex
    # We need to find intersection geometries between hexes and valid landuse
    # Overlay is necessary to get the area of the intersection
    
    # To be efficient, we first spatial join to filter candidates
    # then overlay only on those
    joined = gpd.sjoin(hex_grid_proj[['hex_id', 'geometry']], landuse_valid[['pop_density', 'geometry']], how='inner', predicate='intersects')
    
    # Actually, overlay is safer for area-weighted density.
    # Let's perform overlay.
    intersections = gpd.overlay(hex_grid_proj[['hex_id', 'geometry']], 
                                landuse_valid[['pop_density', 'id', 'geometry']], 
                                how='intersection')
    
    intersections['inter_area'] = intersections.geometry.area
    intersections['weighted_pop'] = intersections['pop_density'] * intersections['inter_area']
    
    # Group by hex_id
    hex_stats = intersections.groupby('hex_id').agg(
        sum_weighted_pop=('weighted_pop', 'sum'),
        sum_inter_area=('inter_area', 'sum'),
        n_overlap_polygons=('id', 'nunique')
    ).reset_index()
    
    hex_stats['area_weighted_density'] = hex_stats['sum_weighted_pop'] / hex_stats['sum_inter_area']
    
    # 3. Count slivers per hex
    # Use sjoin for efficiency
    landuse_slivers['sli_id'] = range(len(landuse_slivers))
    sliver_sjoin = gpd.sjoin(hex_grid_proj[['hex_id', 'geometry']], landuse_slivers[['sli_id', 'geometry']], how='left', predicate='intersects')
    hex_sliver_counts = sliver_sjoin.groupby('hex_id')['sli_id'].nunique().reset_index()
    hex_sliver_counts.columns = ['hex_id', 'n_slivers_filtered']
    
    # 4. Combine results
    # Start with the full hex grid to ensure we handle all cells correctly for ranking
    final_df = hex_grid_proj[['hex_id', 'geometry']].merge(hex_stats, on='hex_id', how='left')
    final_df = final_df.merge(hex_sliver_counts, on='hex_id', how='left')
    
    # Fill NaNs
    final_df['area_weighted_density'] = final_df['area_weighted_density'].fillna(0.0)
    final_df['n_overlap_polygons'] = final_df['n_overlap_polygons'].fillna(0).astype(int)
    final_df['n_slivers_filtered'] = final_df['n_slivers_filtered'].fillna(0).astype(int)
    
    # 5. Ranking and Hotspot selection
    # Sort by density descending
    final_df = final_df.sort_values('area_weighted_density', ascending=False).reset_index(drop=True)
    # Assign rank 1..N
    final_df['rank'] = range(1, len(final_df) + 1)
    
    # Identify top 10%
    n_total = len(final_df)
    n_top = max(1, int(n_total * 0.10))
    # Note: some people define top 10% as floor, some as ceil. 
    # Standard "top 10%" usually means taking the top 10% of the count.
    
    hotspots = final_df.head(n_top).copy()
    
    # 6. Output preparation
    # Output 1: hotspots.geoparquet
    # columns: hex_id, rank, area_weighted_density, geometry
    # Must be in Nigeria's national grid (EPSG:26331)
    hotspots_geo = hotspots[['hex_id', 'rank', 'area_weighted_density', 'geometry']].to_crs('EPSG:26331')
    hotspots_geo.to_parquet('hotspots.geoparquet')
    
    # Output 2: hotspot_ranking.parquet
    # columns: hex_id, rank, area_weighted_density, n_overlap_polygons, n_slivers_filtered
    hotspots_rank = hotspots[['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered']]
    hotspots_rank.to_parquet('hotspot_ranking.parquet')
    
    print(f"Total hexes: {n_total}")
    print(f"Hotspots: {len(hotspots)}")
    print(f"Top density: {hotspots['area_weighted_density'].max()}")
    print(f"Min density in hotspots: {hotspots['area_weighted_density'].min()}")

if __name__ == "__main__":
    solve()