spa-l2-lagos-hotspot-overlaps

Rerun

import geopandas as gpd

# Check the CRS and columns of the input files
landuse = gpd.read_file('lagos_landuse.geojson')
hex_grid = gpd.read_file('lagos_hex_grid.geojson')

print("--- Landuse ---")
print(landuse.head())
print(landuse.crs)
print(landuse.columns)

print("\n--- Hex Grid ---")
print(hex_grid.head())
print(hex_grid.crs)
print(hex_grid.columns)

import geopandas as gpd
import pandas as pd
import numpy as np

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

    # 2. CRS handling
    # We need to work in a metric CRS for area calculations.
    # Nigeria's national grid (UTM Zone 31N or similar) is often EPSG:26331 or EPSG:32631.
    # Let's use EPSG:32631 (WGS 84 / UTM zone 31N) as it's standard for Lagos.
    target_crs = 'EPSG:32631'
    
    landuse_proj = landuse.to_crs(target_crs)
    hex_grid_proj = hex_grid.to_crs(target_crs)

    # 3. Filter slivers (< 100 m2)
    # We must count how many were filtered per hex later.
    # It's easier to calculate landuse area first.
    landuse_proj['landuse_area'] = landuse_proj.geometry.area
    
    # Pre-calculate sliver mask
    is_sliver = landuse_proj['landuse_area'] < 100
    n_slivers_total = is_sliver.sum()
    
    # Let's actually track slivers per hex. 
    # To do this properly, we need to know which hex a sliver belongs to.
    # But a sliver might overlap multiple hexes. 
    # The prompt says "n_slivers_filtered (integer)" in the ranking file.
    # Since it's "per cell", we should count how many slivers overlap each cell.
    
    landuse_clean = landuse_proj[~is_sliver].copy()

    # 4. Spatial Join
    # Join hexes with cleaned landuse
    # We need to calculate area-weighted mean: sum(density_i * area_intersect_i) / sum(area_intersect_i)
    # intersection = overlay(hex, landuse)
    
    # To count slivers per hex:
    # We can do an intersection between hex_grid and the slivers.
    slivers = landuse_proj[is_sliver].copy()
    hex_slivers = gpd.sjoin(slivers, hex_grid_proj, how='inner', predicate='intersects')
    sliver_counts = hex_slivers.groupby('hex_id').size().reset_index(name='n_slivers_filtered')

    # Now the main calculation:
    # Intersection of hex_grid and landuse_clean
    intersections = gpd.overlay(hex_grid_proj, landuse_clean, how='intersection')
    
    # Area of intersection in metric CRS
    intersections['inter_area'] = intersections.geometry.area
    
    # Weighted density: density * area
    intersections['weighted_pop'] = intersections['pop_density'] * intersections['inter_area']
    
    # Group by hex_id to get the weighted mean
    grouped = intersections.groupby('hex_id').agg(
        total_weighted_pop=('weighted_pop', 'sum'),
        total_inter_area=('inter_area', 'sum'),
        n_overlap_polygons=('id', 'count') # Using 'id' from landuse as proxy for unique polygons
    ).reset_index()
    
    # Area-weighted mean = sum(pop * area) / sum(area)
    # If a hex has no intersection area (though overlay handles this), avoid division by zero.
    grouped['area_weighted_density'] = grouped['total_weighted_pop'] / grouped['total_inter_area']

    # 5. Select Top 10%
    # The ranking is based on area_weighted_density.
    # "Top 10% of hex cells ranked by area-weighted mean population density"
    # Usually this means the 10% of cells with the highest density.
    
    num_to_keep = max(1, int(np.ceil(len(grouped) * 0.10)))
    
    # Sort by density descending
    grouped = grouped.sort_values('area_weighted_density', ascending=False)
    
    # We take the top 10%
    top_hexes = grouped.head(num_to_keep).copy()
    
    # 6. Ranking
    # "rank must be a unique integer starting at 1 where ascending rank corresponds to descending area_weighted_density"
    top_hexes['rank'] = range(1, len(top_hexes) + 1)
    
    # 7. Finalize DataFrames
    # Merge sliver counts back
    top_hexes = top_hexes.merge(sliver_counts, on='hex_id', how='left')
    top_hexes['n_slivers_filtered'] = top_hexes['n_slivers_filtered'].fillna(0).astype(int)
    
    # The files must share EXACTLY the same set of hex_id values.
    # The instructions say: "I need two output files that share exactly the same set of hex_id values, one row per cell."
    # This refers to the top 10% cells.
    
    # Prepare hotspot_ranking.parquet
    # columns: hex_id, rank, area_weighted_density, n_overlap_polygons, n_slivers_filtered
    ranking_df = top_hexes[['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered']]
    
    # Prepare hotspots.geoparquet
    # columns: hex_id, rank, area_weighted_density, geometry
    # Needs to be in Nigeria's national grid (EPSG:32631)
    # We need the geometry from the original hex_grid.
    
    hotspots_gdf = hex_grid_proj.merge(top_hexes[['hex_id', 'rank', 'area_weighted_density']], on='hex_id')
    hotspots_gdf = hotspots_gdf[['hex_id', 'rank', 'area_weighted_density', 'geometry']]
    
    # 8. Save files
    ranking_df.to_parquet('hotspot_ranking.parquet', index=False)
    hotspots_gdf.to_geoparquet('hotspots.geoparquet')

    print("Success.")
    print(f"Top 10% size: {len(top_hexes)}")
    print(f"Rankings: {ranking_df.head()}")

if __name__ == "__main__":
    solve()

import geopandas as gpd
import pandas as pd

# Load files
hotspots_gdf = gpd.read_file('hotspots.geoparquet')
ranking_df = pd.read_parquet('hotspot_ranking.parquet')

print("--- Hotspots GDF ---")
print(hotspots_gdf.head())
print(hotspots_gdf.columns)
print(hotspots_gdf.crs)

print("\n--- Ranking DF ---")
print(ranking_df.head())
print(ranking_df.columns)

# Verify shared hex_id
shared_ids = set(hotspots_gdf['hex_id']) == set(ranking_df['hex_id'])
print(f"\nShared hex_id set: {shared_ids}")

# Verify rank is unique integer starting at 1
rank_check = (ranking_df['rank'].dtype == 'int64') and (ranking_df['rank'].min() == 1) and (ranking_df['rank'].is_unique)
print(f"Rank check: {rank_check}")

# Verify n_overlap_polygons and n_slivers_filtered are integer
print(f"n_overlap_polygons type: {ranking_df['n_overlap_polygons'].dtype}")
print(f"n_slivers_filtered type: {ranking_df['n_slivers_filtered'].dtype}")

import pandas as pd
import geopandas as gpd

# Load files
# For GeoParquet, use pandas.read_parquet or geopandas.read_file (if supported)
# Since pyogrio had issues, let's try pandas for the geometry/metadata if needed or just use geopandas again properly.
# Actually, geopandas.read_file should work if it's a proper GeoParquet.
# Let's check what the file actually is.

import os
print(f"File size: {os.path.getsize('hotspots.geoparquet')}")

# Try reading with pandas first to see if it's a valid parquet
try:
    df = pd.read_parquet('hotspots.geoparquet')
    print("Read with pandas success")
    print(df.head())
except Exception as e:
    print(f"Read with pandas failed: {e}")

# Try reading with geopandas properly
try:
    gdf = gpd.read_file('hotspots.geoparquet')
    print("Read with geopandas success")
    print(gdf.head())
except Exception as e:
    print(f"Read with geopandas failed: {e}")

# Also check ranking_df
try:
    df_ranking = pd.read_parquet('hotspot_ranking.parquet')
    print("Read ranking success")
    print(df_ranking.head())
except Exception as e:
    print(f"Read ranking failed: {e}")

import pandas as pd
import geopandas as gpd
from shapely import wkb

# Read with pandas
df = pd.read_parquet('hotspots.geoparquet')

# The geometry column seems to be WKB bytes. Let's convert it to actual geometries.
# Note: GeoParquet usually stores geometry in a way that geopandas understands via pyogrio or fiona.
# If it was saved with to_parquet, it might just be a regular parquet with a WKB column.
# But 'hotspots.geoparquet' implies it should be a GeoParquet.

# Let's re-run the solve.py but use the proper geopandas to_parquet for GeoParquet.
# Actually, for modern GeoPandas, to_parquet() IS the way to save GeoParquet.
# The error in verify_outputs_v1 was probably due to pyogrio's driver detection.

# Let's try to manually reconstruct it to be sure for the user.
# I'll rewrite solve.py to use a more explicit approach if needed.
# But let's first see if we can fix the output.
print("Column types:")
print(df.dtypes)
print(df.head())

import geopandas as gpd
import pandas as pd
import numpy as np

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

    # 2. CRS handling
    target_crs = 'EPSG:32631'
    landuse_proj = landuse.to_crs(target_crs)
    hex_grid_proj = hex_grid.to_crs(target_crs)

    # 3. Filter slivers (< 100 m2)
    landuse_proj['landuse_area'] = landuse_proj.geometry.area
    is_sliver = landuse_proj['landuse_area'] < 100
    
    # Count slivers per hex
    slivers = landuse_proj[is_sliver].copy()
    hex_slivers = gpd.sjoin(slivers, hex_grid_proj, how='inner', predicate='intersects')
    sliver_counts = hex_slivers.groupby('hex_id').size().reset_index(name='n_slivers_filtered')

    landuse_clean = landuse_proj[~is_sliver].copy()

    # 4. Spatial Join / Intersection
    intersections = gpd.overlay(hex_grid_proj, landuse_clean, how='intersection')
    intersections['inter_area'] = intersections.geometry.area
    intersections['weighted_pop'] = intersections['pop_density'] * intersections['inter_area']
    
    grouped = intersections.groupby('hex_id').agg(
        total_weighted_pop=('weighted_pop', 'sum'),
        total_inter_area=('inter_area', 'sum'),
        n_overlap_polygons=('id', 'count')
    ).reset_index()
    
    grouped['area_weighted_density'] = grouped['total_weighted_pop'] / grouped['total_inter_area']

    # 5. Select Top 10%
    num_to_keep = max(1, int(np.ceil(len(grouped) * 0.10)))
    grouped = grouped.sort_values('area_weighted_density', ascending=False)
    top_hexes = grouped.head(num_to_keep).copy()
    
    # 6. Ranking
    top_hexes['rank'] = range(1, len(top_hexes) + 1)
    
    # 7. Finalize DataFrames
    top_hexes = top_hexes.merge(sliver_counts, on='hex_id', how='left')
    top_hexes['n_slivers_filtered'] = top_hexes['n_slivers_filtered'].fillna(0).astype(int)
    
    # Prepare ranking_df
    ranking_df = top_hexes[['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered']]
    
    # Prepare hotspots_gdf
    # We must use the original hex_grid to get the correct geometry
    # and make sure it's a GeoDataFrame.
    hotspots_gdf = hex_grid_proj.merge(top_hexes[['hex_id', 'rank', 'area_weighted_density']], on='hex_id')
    hotspots_gdf = hotspots_gdf[['hex_id', 'rank', 'area_weighted_density', 'geometry']]
    
    # 8. Save files
    # For GeoParquet, it is better to use the GeoDataFrame.to_parquet() method
    # which is supported in newer geopandas/pyarrow.
    ranking_df.to_parquet('hotspot_ranking.parquet', index=False)
    hotspots_gdf.to_parquet('hotspots.geoparquet')

    print("Success.")

if __name__ == "__main__":
    solve()

import pandas as pd
import geopandas as gpd

# Load files
try:
    hotspots_gdf = gpd.read_file('hotspots.geoparquet')
    print("GDF loaded successfully")
except Exception as e:
    print(f"GDF load failed: {e}")

try:
    ranking_df = pd.read_parquet('hotspot_ranking.parquet')
    print("Ranking DF loaded successfully")
except Exception as e:
    print(f"Ranking DF load failed: {e}")

if 'hotspots_gdf' in locals() and 'ranking_df' in locals():
    # Verify shared hex_id
    shared_ids = set(hotspots_gdf['hex_id']) == set(ranking_df['hex_id'])
    print(f"Shared hex_id set: {shared_ids}")

    # Verify rank
    rank_check = (ranking_df['rank'].dtype == 'int64') and (ranking_df['rank'].min() == 1) and (ranking_df['rank'].is_unique)
    print(f"Rank check: {rank_check}")

    # Verify column names
    print(f"GDF columns: {list(hotspots_gdf.columns)}")
    print(f"Ranking columns: {list(ranking_df.columns)}")
    
    # Verify values
    print(f"Top rank density: {ranking_df.iloc[0]['area_weighted_density']}")
    print(f"Bottom rank density: {ranking_df.iloc[-1]['area_weighted_density']}")

import pandas as pd
import geopandas as gpd

# Since geopandas.read_file is failing for the parquet file, let's try reading it with pandas
# and then converting it back to a GeoDataFrame.
df = pd.read_parquet('hotspots.geoparquet')
print("Pandas read successful")

# If it was saved with GeoPandas to_parquet, it should have a geometry column
if 'geometry' in df.columns:
    # Convert WKB/WKT to actual geometries if necessary
    from shapely import wkb
    # Check if it's bytes
    if isinstance(df['geometry'].iloc[0], bytes):
        df['geometry'] = df['geometry'].apply(lambda x: wkb.loads(x))
    
    gdf = gpd.GeoDataFrame(df, geometry='geometry')
    # Set CRS (we used EPSG:32631)
    gdf.set_crs('EPSG:32631', inplace=True)
    print("Converted to GeoDataFrame successfully")
    print(gdf.head())
    print(f"CRS: {gdf.crs}")
else:
    print("No geometry column found")

# Now let's check the ranking_df
ranking_df = pd.read_parquet('hotspot_ranking.parquet')
print("\nRanking DF head:")
print(ranking_df.head())

# Check if they have the same hex_ids
shared_ids = set(gdf['hex_id']) == set(ranking_df['hex_id'])
print(f"\nShared hex_id set: {shared_ids}")