spa-l2-lagos-hotspot-overlaps

Rerun

import geopandas as gpd

# Check files and basic info
landuse = gpd.read_file('lagos_landuse.geojson')
hex_grid = gpd.read_file('lagos_hex_grid.geojson')

print("Landuse info:")
print(landuse.info())
print(landuse.head())
print("\nHexgrid info:")
print(hex_grid.info())
print(hex_grid.head())

print("\nCRS Landuse:", landuse.crs)
print("CRS Hexgrid:", hex_grid.crs)

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')

    # Check CRS
    # They are in EPSG:4326 (degrees). Need a projected CRS for area calculations.
    # Nigeria uses EPSG:26331 (UTM zone 31N) or similar. 
    # Let's use a generic UTM or an appropriate one for Lagos.
    # Lagos is around 6.5 N, 3.4 E. EPSG:32631 is UTM 31N.
    target_crs = 'EPSG:32631'

    # 2. Pre-process Landuse
    # Filter slivers (< 100 m2)
    # Must project to metric CRS first
    landuse_proj = landuse.to_crs(target_crs)
    landuse_proj['area'] = landuse_proj.geometry.area
    
    n_slivers_filtered = (landuse_proj['area'] < 100).sum()
    landuse_clean = landuse_proj[landuse_proj['area'] >= 100].copy()
    
    print(f"Original landuse: {len(landuse)}")
    print(f"Slivers filtered: {n_slivers_filtered}")
    print(f"Clean landuse: {len(landuse_clean)}")

    # 3. Project Hexgrid
    hex_grid_proj = hex_grid.to_crs(target_crs)
    hex_grid_proj['hex_area'] = hex_grid_proj.geometry.area

    # 4. Spatial Join
    # We need to calculate area-weighted mean population density.
    # Area-weighted mean = Sum(pop_density_i * intersection_area_i) / Total_hex_area_covered_by_polygons?
    # Or Sum(pop_density_i * intersection_area_i) / Total_hex_area_of_hex? 
    # Usually, "area-weighted mean population density across overlapping land-use polygons" 
    # implies we only consider the parts of the hex covered by the polygons.
    # However, the prompt says "across overlapping land-use polygons".
    # A standard way is: Mean Density = sum(density_i * area_i) / sum(area_i)
    
    # Perform intersection to get the area of each landuse piece within each hex
    # Using overlay is more precise for area calculations
    intersections = gpd.overlay(hex_grid_proj[['hex_id', 'geometry']], 
                                landuse_clean[['pop_density', 'geometry']], 
                                how='intersection')
    
    # Calculate the area of the intersection pieces
    intersections['inter_area'] = intersections.geometry.area
    
    # Calculate products
    intersections['weighted_pop'] = intersections['pop_density'] * intersections['inter_area']
    
    # Group by hex_id to get the weighted average
    # n_overlap_polygons is the count of landuse polygons overlapping the hex
    # We must be careful: the overlay might split one landuse polygon into multiple pieces in one hex.
    # The prompt says "n_overlap_polygons" (integer). This usually refers to the original landuse polygons.
    
    # Let's get the count of original polygons per hex first.
    # We can use the spatial join for that to be safe, or use the IDs from the overlay if they exist.
    # The overlay 'intersections' will have the 'id' from landuse if it was in landuse_clean.
    # Wait, landuse_clean is a subset of landuse. 
    # Let's re-check columns of landuse. It has 'id'.
    
    # Note: overlay preserves columns from both.
    # Let's group and aggregate.
    agg_results = intersections.groupby('hex_id').agg(
        sum_weighted_pop=('weighted_pop', 'sum'),
        sum_inter_area=('inter_area', 'sum'),
        n_overlap_polygons=('id', 'nunique') # Number of unique landuse polygons
    ).reset_index()
    
    # The prompt asks for 'n_slivers_filtered'. We already calculated it.
    # We also need 'n_overlap_polygons' in the final parquet.
    
    # The prompt also mentions 'n_slivers_filtered' (integer).
    # Let's check how it's used. It's in the parquet.
    
    # Calculate area-weighted density
    # Avoid division by zero
    agg_results['area_weighted_density'] = agg_results['sum_weighted_pop'] / agg_results['sum_inter_area']
    
    # 5. Identify top 10%
    # "top 10 % of hex cells ranked by area-weighted mean population density"
    # This implies we sort by density descending and take the top 10% of the cells.
    
    # Get the hex cells that have at least some overlap? 
    # Or all hex cells? Usually "top 10% of hex cells" refers to the set of cells.
    # However, if a hex has no overlap, its density is effectively 0 or undefined.
    # Let's merge the aggregated results back to the full hex_grid.
    
    final_hex_stats = hex_grid_proj.merge(agg_results, on='hex_id', how='left')
    final_hex_stats['area_weighted_density'] = final_hex_stats['area_weighted_density'].fillna(0)
    final_hex_stats['n_overlap_polygons'] = final_hex_stats['n_overlap_polygons'].fillna(0).astype(int)
    
    # Rank them
    # "rank must be a unique integer starting at 1 where ascending rank corresponds to descending area_weighted_density"
    # Use method='first' to ensure unique ranks for ties.
    final_hex_stats = final_hex_stats.sort_values('area_weighted_density', ascending=False)
    final_hex_stats['rank'] = range(1, len(final_hex_stats) + 1)
    
    # Select top 10%
    num_top = max(1, int(len(final_hex_stats) * 0.10))
    top_10_percent = final_hex_stats.head(num_top).copy()
    
    # 6. Prepare outputs
    # hotspots.geoparquet: hex_id, rank, area_weighted_density, geometry (in Nigeria's national grid)
    # The prompt says "Nigeria's national grid". I'll use EPSG:26331 as it's common, 
    # but I should probably check if there's a specific one. 
    # Actually, the input was EPSG:4326. If it doesn't specify, and "national grid" is mentioned,
    # I will use the projected CRS I used for calculation (32631) or transform back to 4326?
    # Wait, "Nigeria's national grid" often refers to a specific CRS. 
    # However, without a specific code, 32631 (UTM 31N) is the most standard for Lagos.
    # Let's assume the output geometry should be in the same CRS as the calculation (metric) 
    # OR if "national grid" implies something else, I might need to be careful.
    # Most GIS professionals use UTM for "grids" in such briefs.
    
    # Re-reading: "hotspots.geoparquet ... in Nigeria's national grid"
    # I will use EPSG:26331.
    
    hotspots = top_10_percent[['hex_id', 'rank', 'area_weighted_density', 'geometry']].copy()
    hotspots.to_geoparquet('hotspots.geoparquet')
    
    # hotspot_ranking.parquet: hex_id, rank, area_weighted_density, n_overlap_polygons, n_slivers_filtered
    hotspot_ranking = top_10_percent[['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons']].copy()
    hotspot_ranking['n_slivers_filtered'] = n_slivers_filtered
    hotspot_ranking.to_parquet('hotspot_ranking.parquet')

    print("Success!")
    print(f"hotspots.geoparquet shape: {hotspots.shape}")
    print(f"hotspot_ranking.parquet shape: {hotspot_ranking.shape}")

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')

    # Check CRS
    # They are in EPSG:4326 (degrees). Need a projected CRS for area calculations.
    # Nigeria uses EPSG:26331 (UTM zone 31N) or similar. 
    # Let's use a generic UTM or an appropriate one for Lagos.
    # Lagos is around 6.5 N, 3.4 E. EPSG:32631 is UTM 31N.
    target_crs = 'EPSG:32631'

    # 2. Pre-process Landuse
    # Filter slivers (< 100 m2)
    # Must project to metric CRS first
    landuse_proj = landuse.to_crs(target_crs)
    landuse_proj['area'] = landuse_proj.geometry.area
    
    n_slivers_filtered = int((landuse_proj['area'] < 100).sum())
    landuse_clean = landuse_proj[landuse_proj['area'] >= 100].copy()
    
    print(f"Original landuse: {len(landuse)}")
    print(f"Slivers filtered: {n_slivers_filtered}")
    print(f"Clean landuse: {len(landuse_clean)}")

    # 3. Project Hexgrid
    hex_grid_proj = hex_grid.to_crs(target_crs)
    hex_grid_proj['hex_area'] = hex_grid_proj.geometry.area

    # 4. Spatial Join
    # We need to calculate area-weighted mean population density.
    # Area-weighted mean = Sum(pop_density_i * intersection_area_i) / Total_hex_area_covered_by_polygons?
    # Or Sum(pop_density_i * intersection_area_i) / Total_hex_area_of_hex? 
    # Usually, "area-weighted mean population density across overlapping land-use polygons" 
    # implies we only consider the parts of the hex covered by the polygons.
    # However, the prompt says "across overlapping land-use polygons".
    # A standard way is: Mean Density = sum(density_i * area_i) / sum(area_i)
    
    # Perform intersection to get the area of each landuse piece within each hex
    # Using overlay is more precise for area calculations
    intersections = gpd.overlay(hex_grid_proj[['hex_id', 'geometry']], 
                                landuse_clean[['id', 'pop_density', 'geometry']], 
                                how='intersection')
    
    # Calculate the area of the intersection pieces
    intersections['inter_area'] = intersections.geometry.area
    
    # Calculate products
    intersections['weighted_pop'] = intersections['pop_density'] * intersections['inter_area']
    
    # Group by hex_id to get the weighted average
    # n_overlap_polygons is the count of original polygons overlapping the hex
    # We can use the IDs from the overlay if they exist.
    
    # Note: overlay preserves columns from both.
    # Let's group and aggregate.
    agg_results = intersections.groupby('hex_id').agg(
        sum_weighted_pop=('weighted_pop', 'sum'),
        sum_inter_area=('inter_area', 'sum'),
        n_overlap_polygons=('id', 'nunique')
    ).reset_index()
    
    # Calculate area-weighted density
    # Avoid division by zero
    agg_results['area_weighted_density'] = agg_results['sum_weighted_pop'] / agg_results['sum_inter_area']
    
    # 5. Identify top 10%
    # Get the hex cells that have at least some overlap? 
    # Or all hex cells? Usually "top 10% of hex cells" refers to the set of cells.
    
    # Merge the aggregated results back to the full hex_grid.
    final_hex_stats = hex_grid_proj.merge(agg_results, on='hex_id', how='left')
    final_hex_stats['area_weighted_density'] = final_hex_stats['area_weighted_density'].fillna(0.0)
    final_hex_stats['n_overlap_polygons'] = final_hex_stats['n_overlap_polygons'].fillna(0).astype(int)
    
    # Rank them
    # "rank must be a unique integer starting at 1 where ascending rank corresponds to descending area_weighted_density"
    # Use method='first' to ensure unique ranks for ties.
    final_hex_stats = final_hex_stats.sort_values('area_weighted_density', ascending=False)
    final_hex_stats['rank'] = range(1, len(final_hex_stats) + 1)
    
    # Select top 10%
    num_top = max(1, int(len(final_hex_stats) * 0.10))
    top_10_percent = final_hex_stats.head(num_top).copy()
    
    # 6. Prepare outputs
    # hotspots.geoparquet: hex_id, rank, area_weighted_density, geometry (in Nigeria's national grid)
    # Using the projected CRS used for calculation.
    hotspots = gpd.GeoDataFrame(top_10_percent[['hex_id', 'rank', 'area_weighted_density', 'geometry']], 
                                crs=target_crs)
    hotspots.to_parquet('hotspots.geoparquet')
    
    # hotspot_ranking.parquet: hex_id, rank, area_weighted_density, n_overlap_polygons, n_slivers_filtered
    hotspot_ranking = top_10_percent[['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons']].copy()
    hotspot_ranking['n_slivers_filtered'] = n_slivers_filtered
    hotspot_ranking.to_parquet('hotspot_ranking.parquet')

    print("Success!")
    print(f"hotspots.geoparquet shape: {hotspots.shape}")
    print(f"hotspot_ranking.parquet shape: {hotspot_ranking.shape}")

if __name__ == "__main__":
    solve()