spa-l2-lagos-hotspot-overlaps

Rerun

import geopandas as gpd

# Check the first few rows and CRS of both files
try:
    landuse = gpd.read_file('lagos_landuse.geojson')
    print("Landuse columns:", landuse.columns)
    print("Landuse CRS:", landuse.crs)
    print("Landuse first 5 rows:\n", landuse.head())
    print("Landuse area check:", landuse.geometry.area.describe())
except Exception as e:
    print(f"Error reading landuse: {e}")

try:
    hex_grid = gpd.read_file('lagos_hex_grid.geojson')
    print("\nHex grid columns:", hex_grid.columns)
    print("Hex grid CRS:", hex_grid.crs)
    print("Hex grid first 5 rows:\n", hex_grid.head())
except Exception as e:
    print(f"Error reading hex_grid: {e}")

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

    # Nigeria National Grid is EPSG:26332 (UTM zone 32N) or similar.
    # Let's use EPSG:3857 or a local UTM to ensure accurate area calculations.
    # Since the prompt says "Nigeria's national grid", I'll check for a likely candidate 
    # or just use a standard UTM zone for Nigeria if not specified. 
    # Actually, many GIS tasks assume EPSG:3857 for general work, 
    # but let's use EPSG:26332 (UTM 32N) which covers Lagos.
    target_crs = 'EPSG:26332'
    
    # 2. Filter slivers in landuse
    # We must project to metric CRS to filter by area in m^2
    landuse_proj = landuse.to_crs(target_crs)
    
    # Count slivers before filtering
    landuse['is_sliver'] = landuse_proj.geometry.area < 100
    n_slivers_filtered = landuse['is_sliver'].sum()
    
    landuse_clean = landuse[~landuse['is_sliver']].copy()
    landuse_clean_proj = landuse_clean.to_crs(target_crs)

    # 3. Spatial Join to calculate area-weighted density
    # We need the intersection of each hex with each landuse polygon
    # To do this efficiently, we intersect the hexes with landuse
    hex_grid_proj = hex_grid.to_crs(target_crs)
    
    # Perform intersection
    # This returns geometries that are the intersection of the two
    intersections = gpd.overlay(hex_grid_proj, landuse_clean_proj, how='intersection')
    
    # The intersection geometries now have the attributes of both
    # We need to calculate the area of these intersection pieces
    intersections['inter_area'] = intersections.geometry.area
    
    # We also need the original hex area to weight it? 
    # The prompt says "area-weighted mean population density across overlapping land-use polygons"
    # Formula for weighted mean: sum(density_i * area_i) / sum(area_i) 
    # where area_i is the area of the landuse polygon within the hex.
    
    # Calculate numerator: sum(pop_density * inter_area)
    intersections['weighted_pop'] = intersections['pop_density'] * intersections['inter_area']
    
    # Group by hex_id
    grouped = intersections.groupby('hex_id').agg(
        sum_weighted_pop=('weighted_pop', 'sum'),
        sum_inter_area=('inter_area', 'sum'),
        n_overlap_polygons=('id', 'count') # 'id' is from landuse
    ).reset_index()
    
    # Calculate the area-weighted mean population density
    # We divide by the sum of the intersection areas (which is the area of the hex covered by landuse)
    # Note: if sum_inter_area is 0, density should be 0.
    grouped['area_weighted_density'] = grouped['sum_weighted_pop'] / grouped['sum_inter_area']
    
    # 4. Filter top 10%
    # "top 10% of hex cells ranked by area-weighted mean population density"
    # This implies we take the 10% of cells with the highest density.
    threshold = grouped['area_weighted_density'].quantile(0.9)
    top_hexes = grouped[grouped['area_weighted_density'] >= threshold].copy()
    
    # 5. Ranking
    # "rank must be a unique integer starting at 1 where ascending rank corresponds to descending area_weighted_density"
    top_hexes = top_hexes.sort_values(by='area_weighted_density', ascending=False)
    top_hexes['rank'] = range(1, len(top_hexes) + 1)
    
    # Merge back n_slivers_filtered info. Since n_slivers_filtered is a single number 
    # for the whole dataset, we'll just add it as a constant column or handle it per row?
    # "carrying hex_id, rank, area_weighted_density, n_overlap_polygons (integer), and n_slivers_filtered (integer)"
    # This implies n_slivers_filtered is a column in the parquet.
    top_hexes['n_slivers_filtered'] = n_slivers_filtered
    
    # 6. Prepare outputs
    # Output 1: hotspots.geoparquet
    # columns: hex_id, rank, area_weighted_density, geometry
    # Must be in Nigeria's national grid (EPSG:26332)
    
    # We need to join the geometry back from the original hex_grid
    hotspots_geo = hex_grid[hex_grid['hex_id'].isin(top_hexes['hex_id'])].copy()
    hotspots_geo = hotspots_geo.merge(
        top_hexes[['hex_id', 'rank', 'area_weighted_density']], 
        on='hex_id', 
        how='inner'
    )
    
    # Project to target CRS
    hotspots_geo = hotspots_geo.to_crs(target_crs)
    
    # Reorder columns
    hotspots_geo = hotspots_geo[['hex_id', 'rank', 'area_weighted_density', 'geometry']]
    hotspots_geo.to_parquet('hotspots.geoparquet')
    
    # Output 2: hotspot_ranking.parquet
    # columns: hex_id, rank, area_weighted_density, n_overlap_polygons, n_slivers_filtered
    # No geometry
    ranking_df = top_hexes[['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered']]
    ranking_df.to_parquet('hotspot_ranking.parquet')

    print("Files created successfully.")
    print(f"Total top hexes: {len(top_hexes)}")
    print(f"Slivers filtered: {n_slivers_filtered}")

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

    # Nigeria National Grid is EPSG:26332 (UTM zone 32N)
    target_crs = 'EPSG:26332'
    
    # 2. Filter slivers in landuse
    landuse_proj = landuse.to_crs(target_crs)
    
    # Count slivers before filtering
    landuse['is_sliver'] = landuse_proj.geometry.area < 100
    n_slivers_filtered = int(landuse['is_sliver'].sum())
    
    landuse_clean = landuse[~landuse['is_sliver']].copy()
    landuse_clean_proj = landuse_clean.to_crs(target_crs)

    # 3. Spatial Join / Intersection
    hex_grid_proj = hex_grid.to_crs(target_crs)
    
    # Intersection: this yields polygons that are parts of both hex and landuse
    intersections = gpd.overlay(hex_grid_proj, landuse_clean_proj, how='intersection')
    
    # Calculate intersection area
    intersections['inter_area'] = intersections.geometry.area
    
    # Calculate weighted population component
    intersections['weighted_pop'] = intersections['pop_density'] * intersections['inter_area']
    
    # Group by hex_id
    # n_overlap_polygons: count of original landuse polygons inside the hex
    # Note: 'id' is the landuse identifier
    grouped = intersections.groupby('hex_id').agg(
        sum_weighted_pop=('weighted_pop', 'sum'),
        sum_inter_area=('inter_area', 'sum'),
        n_overlap_polygons=('id', 'count')
    ).reset_index()
    
    # Area-weighted mean: sum(density_i * area_i) / sum(area_i)
    # Where sum(area_i) is the total area within the hex that is covered by landuse polygons
    grouped['area_weighted_density'] = grouped['sum_weighted_pop'] / grouped['sum_inter_area']
    
    # 4. Filter top 10%
    # We need the top 10% of hexes by density.
    # To get exactly 10%, we use quantile
    threshold = grouped['area_weighted_density'].quantile(0.9)
    top_hexes_df = grouped[grouped['area_weighted_density'] >= threshold].copy()
    
    # 5. Ranking
    # Sort descending by density
    top_hexes_df = top_hexes_df.sort_values(by='area_weighted_density', ascending=False)
    # Assign unique rank starting at 1
    top_hexes_df['rank'] = range(1, len(top_hexes_df) + 1)
    
    # Add n_slivers_filtered
    top_hexes_df['n_slivers_filtered'] = n_slivers_filtered

    # 6. Prepare outputs
    # Output 1: hotspots.geoparquet
    # columns: hex_id, rank, area_weighted_density, geometry
    # Must match exactly the hex_ids in top_hexes_df
    
    # Join geometry back from the original hex_grid (which is in EPSG:4326)
    # But the prompt asks for the file in Nigeria's national grid
    hotspots_geo = hex_grid.merge(
        top_hexes_df[['hex_id', 'rank', 'area_weighted_density']], 
        on='hex_id', 
        how='inner'
    )
    hotspots_geo = hotspots_geo.to_crs(target_crs)
    hotspots_geo = hotspots_geo[['hex_id', 'rank', 'area_weighted_density', 'geometry']]
    hotspots_geo.to_parquet('hotspots.geoparquet')
    
    # Output 2: hotspot_ranking.parquet
    # columns: hex_id, rank, area_weighted_density, n_overlap_polygons, n_slivers_filtered
    # No geometry. 
    # IMPORTANT: Must share EXACTLY the same hex_id values as hotspots.geoparquet
    # Since we merged on hex_id using inner join, it should be correct.
    ranking_df = top_hexes_df[['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered']]
    ranking_df.to_parquet('hotspot_ranking.parquet')

    print("Files created successfully.")
    print(f"Total top hexes: {len(top_hexes_df)}")
    print(f"Slivers filtered: {n_slivers_filtered}")

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

    # Nigeria National Grid is EPSG:26332 (UTM zone 32N)
    target_crs = 'EPSG:26332'
    
    # 2. Filter slivers in landuse
    landuse_proj = landuse.to_crs(target_crs)
    
    # Count slivers before filtering
    landuse['is_sliver'] = landuse_proj.geometry.area < 100
    n_slivers_filtered = int(landuse['is_sliver'].sum())
    
    landuse_clean = landuse[~landuse['is_sliver']].copy()
    landuse_clean_proj = landuse_clean.to_crs(target_crs)

    # 3. Spatial Join / Intersection
    hex_grid_proj = hex_grid.to_crs(target_crs)
    
    # Intersection: this yields polygons that are parts of both hex and landuse
    intersections = gpd.overlay(hex_grid_proj, landuse_clean_proj, how='intersection')
    
    # Calculate intersection area
    intersections['inter_area'] = intersections.geometry.area
    
    # Calculate weighted population component
    intersections['weighted_pop'] = intersections['pop_density'] * intersections['inter_area']
    
    # Group by hex_id
    # n_overlap_polygons: count of original landuse polygons inside the hex
    # Note: 'id' is the landuse identifier
    grouped = intersections.groupby('hex_id').agg(
        sum_weighted_pop=('weighted_pop', 'sum'),
        sum_inter_area=('inter_area', 'sum'),
        n_overlap_polygons=('id', 'count')
    ).reset_index()
    
    # Area-weighted mean: sum(density_i * area_i) / sum(area_i)
    # Where sum(area_i) is the total area within the hex that is covered by landuse polygons
    grouped['area_weighted_density'] = grouped['sum_weighted_pop'] / grouped['sum_inter_area']
    
    # 4. Filter top 10%
    # We need the top 10% of hexes by density.
    # To get exactly 10%, we use quantile
    threshold = grouped['area_weighted_density'].quantile(0.9)
    top_hexes_df = grouped[grouped['area_weighted_density'] >= threshold].copy()
    
    # 5. Ranking
    # Sort descending by density
    top_hexes_df = top_hexes_df.sort_values(by='area_weighted_density', ascending=False)
    # Assign unique rank starting at 1
    top_hexes_df['rank'] = range(1, len(top_hexes_df) + 1)
    
    # Add n_slivers_filtered
    top_hexes_df['n_slivers_filtered'] = n_slivers_filtered

    # 6. Prepare outputs
    # Output 1: hotspots.geoparquet
    # columns: hex_id, rank, area_weighted_density, geometry
    # Must match exactly the hex_ids in top_hexes_df
    
    # Join geometry back from the original hex_grid (which is in EPSG:4326)
    # We'll merge the top_hexes_df with hex_grid to get geometry.
    # hex_grid contains the original geometries in EPSG:4326.
    
    hotspots_geo = hex_grid.merge(
        top_hexes_df[['hex_id', 'rank', 'area_weighted_density']], 
        on='hex_id', 
        how='inner'
    )
    hotspots_geo = hotspots_geo.to_crs(target_crs)
    hotspots_geo = hotspots_geo[['hex_id', 'rank', 'area_weighted_density', 'geometry']]
    hotspots_geo.to_parquet('hotspots.geoparquet')
    
    # Output 2: hotspot_ranking.parquet
    # columns: hex_id, rank, area_weighted_density, n_overlap_polygons, n_slivers_filtered
    # No geometry. 
    # IMPORTANT: Must share EXACTLY the same hex_id values as hotspots.geoparquet
    
    # To ensure they match exactly, let's define the final set of hex_ids
    final_hex_ids = hotspots_geo['hex_id'].unique()
    
    # Filter top_hexes_df to include only these hex_ids and re-order to match hotspots_geo
    # (Though order doesn't strictly matter for "sharing the same set", let's be safe)
    ranking_df = top_hexes_df[top_hexes_df['hex_id'].isin(final_hex_ids)].copy()
    
    # Ensure column order matches prompt exactly
    ranking_df = ranking_df[['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered']]
    ranking_df.to_parquet('hotspot_ranking.parquet')

    print("Files created successfully.")
    print(f"Total top hexes: {len(top_hexes_df)}")
    print(f"Slivers filtered: {n_slivers_filtered}")

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

    # Nigeria National Grid is EPSG:26332 (UTM zone 32N)
    target_crs = 'EPSG:26332'
    
    # 2. Filter slivers in landuse
    landuse_proj = landuse.to_crs(target_crs)
    
    # Count slivers before filtering
    landuse['is_sliver'] = landuse_proj.geometry.area < 100
    n_slivers_filtered = int(landuse['is_sliver'].sum())
    
    landuse_clean = landuse[~landuse['is_sliver']].copy()
    landuse_clean_proj = landuse_clean.to_crs(target_crs)

    # 3. Spatial Join / Intersection
    hex_grid_proj = hex_grid.to_crs(target_crs)
    
    # Intersection: this yields polygons that are parts of both hex and landuse
    intersections = gpd.overlay(hex_grid_proj, landuse_clean_proj, how='intersection')
    
    # Calculate intersection area
    intersections['inter_area'] = intersections.geometry.area
    
    # Calculate weighted population component
    intersections['weighted_pop'] = intersections['pop_density'] * intersections['inter_area']
    
    # Group by hex_id
    # n_overlap_polygons: count of original landuse polygons inside the hex
    # Note: 'id' is the landuse identifier
    grouped = intersections.groupby('hex_id').agg(
        sum_weighted_pop=('weighted_pop', 'sum'),
        sum_inter_area=('inter_area', 'sum'),
        n_overlap_polygons=('id', 'count')
    ).reset_index()
    
    # Area-weighted mean: sum(density_i * area_i) / sum(area_i)
    grouped['area_weighted_density'] = grouped['sum_weighted_pop'] / grouped['sum_inter_area']
    
    # 4. Filter top 10%
    # We need the top 10% of hexes by density.
    # To get exactly 10%, we use quantile
    threshold = grouped['area_weighted_density'].quantile(0.9)
    top_hexes_df = grouped[grouped['area_weighted_density'] >= threshold].copy()
    
    # 5. Ranking
    # Sort descending by density
    top_hexes_df = top_hexes_df.sort_values(by='area_weighted_density', ascending=False)
    # Assign unique rank starting at 1
    top_hexes_df['rank'] = range(1, len(top_hexes_df) + 1)
    
    # Add n_slivers_filtered
    top_hexes_df['n_slivers_filtered'] = n_slivers_filtered

    # 6. Prepare outputs
    # Output 1: hotspots.geoparquet
    # columns: hex_id, rank, area_weighted_density, geometry
    
    # Join geometry back from the original hex_grid (which is in EPSG:4326)
    # We'll merge the top_hexes_df with hex_grid to get geometry.
    
    # Ensure we use the exact same set of hex_ids.
    # hotspots_geo will be based on top_hexes_df
    hotspots_geo = hex_grid.merge(
        top_hexes_df[['hex_id', 'rank', 'area_weighted_density']], 
        on='hex_id', 
        how='inner'
    )
    hotspots_geo = hotspots_geo.to_crs(target_crs)
    hotspots_geo = hotspots_geo[['hex_id', 'rank', 'area_weighted_density', 'geometry']]
    hotspots_geo.to_parquet('hotspots.geoparquet')
    
    # Output 2: hotspot_ranking.parquet
    # columns: hex_id, rank, area_weighted_density, n_overlap_polygons, n_slivers_filtered
    
    # Since we used inner join to create hotspots_geo, we can just use the same
    # top_hexes_df which was already filtered.
    # Let's ensure the order is the same if we want to be very careful, 
    # but the prompt just says "share exactly the same set of hex_id values".
    
    ranking_df = top_hexes_df[['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered']]
    
    # Re-sorting ranking_df to ensure it's in the same order as hotspots_geo if needed, 
    # but for Parquet it's not strictly required.
    # Let's just ensure the IDs are present.
    ranking_df = ranking_df.sort_values('hex_id')
    hotspots_geo_sorted = hotspots_geo.sort_values('hex_id')
    
    print("Checking ID match after merge...")
    print(f"hotspots hex_ids: {len(hotspots_geo['hex_id'].unique())}")
    print(f"ranking hex_ids: {len(ranking_df['hex_id'].unique())}")

    ranking_df.to_parquet('hotspot_ranking.parquet')

    print("Files created successfully.")
    print(f"Total top hexes: {len(top_hexes_df)}")
    print(f"Slivers filtered: {n_slivers_filtered}")

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

    # Nigeria National Grid is EPSG:26332 (UTM zone 32N)
    target_crs = 'EPSG:26332'
    
    # 2. Filter slivers in landuse
    landuse_proj = landuse.to_crs(target_crs)
    
    # Count slivers before filtering
    landuse['is_sliver'] = landuse_proj.geometry.area < 100
    n_slivers_filtered = int(landuse['is_sliver'].sum())
    
    landuse_clean = landuse[~landuse['is_sliver']].copy()
    landuse_clean_proj = landuse_clean.to_crs(target_crs)

    # 3. Spatial Join / Intersection
    hex_grid_proj = hex_grid.to_crs(target_crs)
    
    # Intersection: this yields polygons that are parts of both hex and landuse
    intersections = gpd.overlay(hex_grid_proj, landuse_clean_proj, how='intersection')
    
    # Calculate intersection area
    intersections['inter_area'] = intersections.geometry.area
    
    # Calculate weighted population component
    intersections['weighted_pop'] = intersections['pop_density'] * intersections['inter_area']
    
    # Group by hex_id
    # n_overlap_polygons: count of original landuse polygons inside the hex
    # Note: 'id' is the landuse identifier
    grouped = intersections.groupby('hex_id').agg(
        sum_weighted_pop=('weighted_pop', 'sum'),
        sum_inter_area=('inter_area', 'sum'),
        n_overlap_polygons=('id', 'count')
    ).reset_index()
    
    # Area-weighted mean: sum(density_i * area_i) / sum(area_i)
    grouped['area_weighted_density'] = grouped['sum_weighted_pop'] / grouped['sum_inter_area']
    
    # 4. Filter top 10%
    # We need the top 10% of hexes by density.
    # To get exactly 10%, we use quantile
    threshold = grouped['area_weighted_density'].quantile(0.9)
    top_hexes_df = grouped[grouped['area_weighted_density'] >= threshold].copy()
    
    # 5. Ranking
    # Sort descending by density
    top_hexes_df = top_hexes_df.sort_values(by='area_weighted_density', ascending=False)
    # Assign unique rank starting at 1
    top_hexes_df['rank'] = range(1, len(top_hexes_df) + 1)
    
    # Add n_slivers_filtered
    top_hexes_df['n_slivers_filtered'] = n_slivers_filtered

    # 6. Prepare outputs
    # Output 1: hotspots.geoparquet
    # columns: hex_id, rank, area_weighted_density, geometry
    
    # Join geometry back from the original hex_grid (which is in EPSG:4326)
    # We'll merge the top_hexes_df with hex_grid to get geometry.
    
    # To ensure they match exactly, we'll perform the join on hex_id.
    # We want to keep only the hex_ids in top_hexes_df.
    
    # We'll sort top_hexes_df by hex_id to ensure consistent ordering for comparison
    top_hexes_df_sorted = top_hexes_df.sort_values('hex_id')
    
    hotspots_geo = hex_grid.merge(
        top_hexes_df[['hex_id', 'rank', 'area_weighted_density']], 
        on='hex_id', 
        how='inner'
    )
    hotspots_geo = hotspots_geo.to_crs(target_crs)
    hotspots_geo = hotspots_geo[['hex_id', 'rank', 'area_weighted_density', 'geometry']]
    hotspots_geo.to_parquet('hotspots.geoparquet')
    
    # Output 2: hotspot_ranking.parquet
    # columns: hex_id, rank, area_weighted_density, n_overlap_polygons, n_slivers_filtered
    
    # Using top_hexes_df directly which contains all the required columns
    ranking_df = top_hexes_df[['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered']]
    
    # Sort both to make comparison easy
    ranking_df = ranking_df.sort_values('hex_id')
    hotspots_geo_ids = hotspots_geo['hex_id'].sort_values()
    ranking_ids = ranking_df['hex_id'].sort_values()
    
    # Check match
    match = hotspots_geo_ids.equals(ranking_ids)
    print(f"ID match check: {match}")

    ranking_df.to_parquet('hotspot_ranking.parquet')

    print("Files created successfully.")
    print(f"Total top hexes: {len(top_hexes_df)}")
    print(f"Slivers filtered: {n_slivers_filtered}")

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

    # Nigeria National Grid is EPSG:26332 (UTM zone 32N)
    target_crs = 'EPSG:26332'
    
    # 2. Filter slivers in landuse
    landuse_proj = landuse.to_crs(target_crs)
    
    # Count slivers before filtering
    landuse['is_sliver'] = landuse_proj.geometry.area < 100
    n_slivers_filtered = int(landuse['is_sliver'].sum())
    
    landuse_clean = landuse[~landuse['is_sliver']].copy()
    landuse_clean_proj = landuse_clean.to_crs(target_crs)

    # 3. Spatial Join / Intersection
    hex_grid_proj = hex_grid.to_crs(target_crs)
    
    # Intersection: this yields polygons that are parts of both hex and landuse
    intersections = gpd.overlay(hex_grid_proj, landuse_clean_proj, how='intersection')
    
    # Calculate intersection area
    intersections['inter_area'] = intersections.geometry.area
    
    # Calculate weighted population component
    intersections['weighted_pop'] = intersections['pop_density'] * intersections['inter_area']
    
    # Group by hex_id
    # n_overlap_polygons: count of original landuse polygons inside the hex
    # Note: 'id' is the landuse identifier
    grouped = intersections.groupby('hex_id').agg(
        sum_weighted_pop=('weighted_pop', 'sum'),
        sum_inter_area=('inter_area', 'sum'),
        n_overlap_polygons=('id', 'count')
    ).reset_index()
    
    # Area-weighted mean: sum(density_i * area_i) / sum(area_i)
    grouped['area_weighted_density'] = grouped['sum_weighted_pop'] / grouped['sum_inter_area']
    
    # 4. Filter top 10%
    # We need the top 10% of hexes by density.
    # To get exactly 10%, we use quantile
    threshold = grouped['area_weighted_density'].quantile(0.9)
    top_hexes_df = grouped[grouped['area_weighted_density'] >= threshold].copy()
    
    # 5. Ranking
    # Sort descending by density
    top_hexes_df = top_hexes_df.sort_values(by='area_weighted_density', ascending=False)
    # Assign unique rank starting at 1
    top_hexes_df['rank'] = range(1, len(top_hexes_df) + 1)
    
    # Add n_slivers_filtered
    top_hexes_df['n_slivers_filtered'] = n_slivers_filtered

    # 6. Prepare outputs
    # Output 1: hotspots.geoparquet
    # columns: hex_id, rank, area_weighted_density, geometry
    
    # Create a dataframe that will be the source for both outputs to ensure exact match
    # This is our master table for the top hexes
    master_df = top_hexes_df.copy()

    # Create the geo version
    # We'll merge master_df with the original hex_grid to get the geometry in EPSG:4326,
    # then transform to target_crs.
    hotspots_geo = hex_grid.merge(
        master_df[['hex_id', 'rank', 'area_weighted_density']], 
        on='hex_id', 
        how='inner'
    )
    hotspots_geo = hotspots_geo.to_crs(target_crs)
    hotspots_geo = hotspots_geo[['hex_id', 'rank', 'area_weighted_density', 'geometry']]
    hotspots_geo.to_parquet('hotspots.geoparquet')
    
    # Create the ranking version
    # We'll use the same master_df.
    # The prompt specifies: hex_id, rank, area_weighted_density, n_overlap_polygons, n_slivers_filtered
    ranking_df = master_df[['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered']]
    
    # Final validation before writing
    # Ensure the hex_ids in both files are the same.
    # Because we used master_df to drive both, and hex_grid merge is an inner join on hex_id,
    # they MUST be identical.
    
    # Convert to sets for comparison to be absolutely sure.
    set_geo = set(hotspots_geo['hex_id'])
    set_rank = set(ranking_df['hex_id'])
    
    if set_geo != set_rank:
        raise ValueError(f"Mismatch! Geo IDs: {set_geo}, Rank IDs: {set_rank}")
    
    ranking_df.to_parquet('hotspot_ranking.parquet')

    print("Files created successfully.")
    print(f"Total top hexes: {len(master_df)}")
    print(f"Slivers filtered: {n_slivers_filtered}")

if __name__ == "__main__":
    solve()