spa-l2-lagos-hotspot-overlaps

Rerun

     1	import geopandas as gpd
     2	import pandas as pd
     3	import numpy as np
     4	
     5	# National grid CRS for Nigeria (Minna / Nigeria West) — metric, covers Lagos
     6	NIGERIA_GRID = "EPSG:26391"
     7	
     8	# ── 1. Load data ──────────────────────────────────────────────────────────────
     9	print("Loading land-use polygons …")
    10	lu = gpd.read_file("/work/lagos_landuse.geojson")
    11	print(f"  {len(lu)} polygons, CRS={lu.crs}")
    12	
    13	print("Loading hex grid …")
    14	hex_grid = gpd.read_file("/work/lagos_hex_grid.geojson")
    15	print(f"  {len(hex_grid)} cells, CRS={hex_grid.crs}")
    16	
    17	# ── 2. Reproject to metric CRS ────────────────────────────────────────────────
    18	print(f"\nReprojecting to {NIGERIA_GRID} …")
    19	lu = lu.to_crs(NIGERIA_GRID)
    20	hex_grid = hex_grid.to_crs(NIGERIA_GRID)
    21	
    22	# ── 3. Filter slivers (area < 100 m²) ─────────────────────────────────────────
    23	# Compute area while still in metric CRS
    24	lu["area_m2"] = lu.geometry.area
    25	n_before = len(lu)
    26	lu_big = lu[lu["area_m2"] >= 100.0].copy()
    27	n_slivers = n_before - len(lu_big)
    28	print(f"\nSlivers (< 100 m²): {n_slivers} polygons filtered out")
    29	print(f"Remaining land-use polygons: {len(lu_big)}")
    30	
    31	# ── 4. Spatial join: hex → overlapping land-use polygons ──────────────────────
    32	# Use sjoin to get all hex-polygon overlaps
    33	print("\nSpatial join (hex × land-use) …")
    34	joined = gpd.sjoin(
    35	    hex_grid[["hex_id", "geometry"]],
    36	    lu_big[["pop_density", "area_m2", "geometry"]],
    37	    how="inner",
    38	    predicate="intersects",
    39	)
    40	print(f"  {len(joined)} overlap relationships")
    41	
    42	# ── 5. Clip each polygon to its hex cell, compute clipped area ────────────────
    43	# For area-weighted mean we need the area of the *intersection*,
    44	# not the full polygon area.
    45	print("Computing intersection areas …")
    46	
    47	# We'll join again but get the actual intersection geometry
    48	# Using overlay for efficiency
    49	overlay = gpd.overlay(
    50	    hex_grid[["hex_id", "geometry"]],
    51	    lu_big[["pop_density", "geometry"]],
    52	    how="intersection",
    53	    keep_geom_type=False,
    54	)
    55	print(f"  {len(overlay)} intersection fragments")
    56	
    57	# Area of each fragment
    58	overlay["intersect_area_m2"] = overlay.geometry.area
    59	
    60	# ── 6. Compute area-weighted mean density per hex ─────────────────────────────
    61	print("Computing area-weighted mean density per hex cell …")
    62	
    63	# For each hex: sum(pop_density * intersect_area) / sum(intersect_area)
    64	grp = overlay.groupby("hex_id")
    65	agg = grp.agg(
    66	    sum_weighted_density=("intersect_area_m2", lambda s: np.nansum(s.values)),
    67	    # Actually we need weighted sum: density * area
    68	)
    69	# Let's compute properly
    70	overlay["w"] = overlay["pop_density"] * overlay["intersect_area_m2"]
    71	hex_agg = overlay.groupby("hex_id").agg(
    72	    weighted_sum=("w", "sum"),
    73	    total_intersect_area=("intersect_area_m2", "sum"),
    74	    n_overlap_polygons=("pop_density", "count"),
    75	)
    76	hex_agg["area_weighted_density"] = (
    77	    hex_agg["weighted_sum"] / hex_agg["total_intersect_area"]
    78	)
    79	
    80	# ── 7. Count slivers filtered per hex ─────────────────────────────────────────
    81	# Slivers: land-use polygons < 100 m² that overlapped the hex
    82	# We can find these by looking at the original lu polygons with area < 100
    83	lu_slivers = lu[lu["area_m2"] < 100.0].copy()
    84	if len(lu_slivers) > 0:
    85	    sliver_join = gpd.sjoin(
    86	        hex_grid[["hex_id", "geometry"]],
    87	        lu_slivers[["geometry"]],
    88	        how="inner",
    89	        predicate="intersects",
    90	    )
    91	    n_slivers_per_hex = sliver_join.groupby("hex_id").size().rename("n_slivers_filtered")
    92	else:
    93	    n_slivers_per_hex = pd.Series(dtype=int)
    94	
    95	hex_agg = hex_agg.join(n_slivers_per_hex)
    96	hex_agg["n_slivers_filtered"] = hex_agg["n_slivers_filtered"].fillna(0).astype(int)
    97	
    98	# ── 8. Sort by density descending and take top 10 % ───────────────────────────
    99	hex_agg = hex_agg.sort_values("area_weighted_density", ascending=False)
   100	n_total = len(hex_agg)
   101	n_top = max(1, int(np.round(n_total * 0.10)))
   102	print(f"\nTotal hex cells: {n_total}")
   103	print(f"Top 10 %: {n_top} cells")
   104	
   105	top_hex_ids = hex_agg.index[:n_top]
   106	top = hex_agg.loc[top_hex_ids].copy()
   107	
   108	# ── 9. Create rank ────────────────────────────────────────────────────────────
   109	top["rank"] = np.arange(1, len(top) + 1, dtype=int)
   110	
   111	# ── 10. Build output dataframes ───────────────────────────────────────────────
   112	# hotspots.geoparquet — with geometry
   113	result_geo = hex_grid[hex_grid["hex_id"].isin(top_hex_ids)].copy()
   114	result_geo = result_geo.set_index("hex_id")
   115	result_geo = result_geo.join(top[["rank", "area_weighted_density"]])
   116	result_geo = result_geo.reset_index()
   117	result_geo = result_geo[["hex_id", "rank", "area_weighted_density", "geometry"]]
   118	result_geo = result_geo.sort_values("rank")
   119	print(f"\nhotspots.geoparquet: {len(result_geo)} rows")
   120	
   121	# hotspot_ranking.parquet — no geometry
   122	result_plain = top[["rank", "area_weighted_density", "n_overlap_polygons", "n_slivers_filtered"]].copy()
   123	result_plain = result_plain.reset_index()
   124	result_plain = result_plain[["hex_id", "rank", "area_weighted_density", "n_overlap_polygons", "n_slivers_filtered"]]
   125	result_plain = result_plain.sort_values("rank")
   126	print(f"hotspot_ranking.parquet: {len(result_plain)} rows")
   127	
   128	# ── 11. Write outputs ─────────────────────────────────────────────────────────
   129	result_geo.to_parquet("/work/hotspots.geoparquet", index=False)
   130	print("  ✓ hotspots.geoparquet written")
   131	
   132	result_plain.to_parquet("/work/hotspot_ranking.parquet", index=False)
   133	print("  ✓ hotspot_ranking.parquet written")
   134	
   135	# ── 12. Sanity checks ─────────────────────────────────────────────────────────
   136	print("\n─── Sanity checks ───")
   137	# Check hex_ids match
   138	ids_geo = set(result_geo["hex_id"])
   139	ids_plain = set(result_plain["hex_id"])
   140	assert ids_geo == ids_plain, "hex_id sets differ!"
   141	print(f"✓ Both files share the same {len(ids_geo)} hex_id values")
   142	
   143	# Check rank is unique 1..n
   144	assert result_geo["rank"].is_monotonic_increasing, "ranks not sorted ascending"
   145	assert list(result_geo["rank"]) == list(range(1, len(result_geo) + 1)), "ranks not 1..n"
   146	print("✓ Ranks are unique 1..n")
   147	
   148	# Check ordering: descending density
   149	assert all(
   150	    result_geo["area_weighted_density"].iloc[i] >= result_geo["area_weighted_density"].iloc[i+1]
   151	    for i in range(len(result_geo) - 1)
   152	), "density not descending by rank"
   153	print("✓ Density is descending with rank")
   154	
   155	# Verify CRS of geo output
   156	print(f"✓ Output geometry CRS: {result_geo.crs}")
   157	
   158	# Print top 5 rows
   159	print("\nTop 5 ranked cells:")
   160	print(result_geo.head(5).to_string(index=False))
   161	print()
   162	print(result_plain.head(5).to_string(index=False))

     1	import geopandas as gpd
     2	import pandas as pd
     3	import numpy as np
     4	
     5	# National grid CRS for Nigeria (Minna / Nigeria West) — metric, covers Lagos
     6	NIGERIA_GRID = "EPSG:26391"
     7	
     8	# ── 1. Load data ──────────────────────────────────────────────────────────────
     9	print("Loading land-use polygons …")
    10	lu = gpd.read_file("/work/lagos_landuse.geojson")
    11	print(f"  {len(lu)} polygons, CRS={lu.crs}")
    12	
    13	print("Loading hex grid …")
    14	hex_grid = gpd.read_file("/work/lagos_hex_grid.geojson")
    15	print(f"  {len(hex_grid)} cells, CRS={hex_grid.crs}")
    16	
    17	# ── 2. Reproject to metric CRS ────────────────────────────────────────────────
    18	print(f"\nReprojecting to {NIGERIA_GRID} …")
    19	lu = lu.to_crs(NIGERIA_GRID)
    20	hex_grid = hex_grid.to_crs(NIGERIA_GRID)
    21	
    22	# ── 3. Filter slivers (area < 100 m²) ─────────────────────────────────────────
    23	# Compute area while still in metric CRS
    24	lu["area_m2"] = lu.geometry.area
    25	n_before = len(lu)
    26	lu_big = lu[lu["area_m2"] >= 100.0].copy()
    27	n_slivers = n_before - len(lu_big)
    28	print(f"\nSlivers (< 100 m²): {n_slivers} polygons filtered out")
    29	print(f"Remaining land-use polygons: {len(lu_big)}")
    30	
    31	# ── 4. Spatial join: hex → overlapping land-use polygons ──────────────────────
    32	# Use sjoin to get all hex-polygon overlaps
    33	print("\nSpatial join (hex × land-use) …")
    34	joined = gpd.sjoin(
    35	    hex_grid[["hex_id", "geometry"]],
    36	    lu_big[["pop_density", "area_m2", "geometry"]],
    37	    how="inner",
    38	    predicate="intersects",
    39	)
    40	print(f"  {len(joined)} overlap relationships")
    41	
    42	# ── 5. Clip each polygon to its hex cell, compute clipped area ────────────────
    43	# For area-weighted mean we need the area of the *intersection*,
    44	# not the full polygon area.
    45	print("Computing intersection areas …")
    46	
    47	# We'll join again but get the actual intersection geometry
    48	# Using overlay for efficiency
    49	overlay = gpd.overlay(
    50	    hex_grid[["hex_id", "geometry"]],
    51	    lu_big[["pop_density", "geometry"]],
    52	    how="intersection",
    53	    keep_geom_type=False,
    54	)
    55	print(f"  {len(overlay)} intersection fragments")
    56	
    57	# Area of each fragment
    58	overlay["intersect_area_m2"] = overlay.geometry.area
    59	
    60	# ── 6. Compute area-weighted mean density per hex ─────────────────────────────
    61	print("Computing area-weighted mean density per hex cell …")
    62	
    63	# For each hex: sum(pop_density * intersect_area) / sum(intersect_area)
    64	grp = overlay.groupby("hex_id")
    65	agg = grp.agg(
    66	    sum_weighted_density=("intersect_area_m2", lambda s: np.nansum(s.values)),
    67	    # Actually we need weighted sum: density * area
    68	)
    69	# Let's compute properly
    70	overlay["w"] = overlay["pop_density"] * overlay["intersect_area_m2"]
    71	hex_agg = overlay.groupby("hex_id").agg(
    72	    weighted_sum=("w", "sum"),
    73	    total_intersect_area=("intersect_area_m2", "sum"),
    74	    n_overlap_polygons=("pop_density", "count"),
    75	)
    76	hex_agg["area_weighted_density"] = (
    77	    hex_agg["weighted_sum"] / hex_agg["total_intersect_area"]
    78	)
    79	
    80	# ── 7. Count slivers filtered per hex ─────────────────────────────────────────
    81	# Slivers: land-use polygons < 100 m² that overlapped the hex
    82	# We can find these by looking at the original lu polygons with area < 100
    83	lu_slivers = lu[lu["area_m2"] < 100.0].copy()
    84	if len(lu_slivers) > 0:
    85	    sliver_join = gpd.sjoin(
    86	        hex_grid[["hex_id", "geometry"]],
    87	        lu_slivers[["geometry"]],
    88	        how="inner",
    89	        predicate="intersects",
    90	    )
    91	    n_slivers_per_hex = sliver_join.groupby("hex_id").size().rename("n_slivers_filtered")
    92	else:
    93	    n_slivers_per_hex = pd.Series(dtype=int)
    94	
    95	hex_agg = hex_agg.join(n_slivers_per_hex)
    96	hex_agg["n_slivers_filtered"] = hex_agg["n_slivers_filtered"].fillna(0).astype(int)
    97	
    98	# ── 8. Merge back to all hex cells, rank all, take top 10 % ──────────────────
    99	# Cells with no overlap get NaN density — they rank at the bottom
   100	all_hex = hex_grid[["hex_id"]].copy()
   101	all_hex = all_hex.set_index("hex_id")
   102	all_hex = all_hex.join(hex_agg[["area_weighted_density", "n_overlap_polygons", "n_slivers_filtered"]])
   103	all_hex["n_overlap_polygons"] = all_hex["n_overlap_polygons"].fillna(0).astype(int)
   104	all_hex["n_slivers_filtered"] = all_hex["n_slivers_filtered"].fillna(0).astype(int)
   105	
   106	# Sort: NaN densities go to the bottom
   107	all_hex = all_hex.sort_values("area_weighted_density", ascending=False, na_position="last")
   108	
   109	n_total = len(all_hex)
   110	n_top = max(1, int(np.round(n_total * 0.10)))
   111	print(f"\nTotal hex cells: {n_total}")
   112	print(f"Top 10 %: {n_top} cells")
   113	
   114	top_hex_ids = all_hex.index[:n_top]
   115	top = all_hex.loc[top_hex_ids].copy()
   116	
   117	top_hex_ids = hex_agg.index[:n_top]
   118	top = hex_agg.loc[top_hex_ids].copy()
   119	
   120	# ── 9. Create rank ────────────────────────────────────────────────────────────
   121	top["rank"] = np.arange(1, len(top) + 1, dtype=int)
   122	
   123	# ── 10. Build output dataframes ───────────────────────────────────────────────
   124	# hotspots.geoparquet — with geometry
   125	result_geo = hex_grid[hex_grid["hex_id"].isin(top_hex_ids)].copy()
   126	result_geo = result_geo.set_index("hex_id")
   127	result_geo = result_geo.join(top[["rank", "area_weighted_density"]])
   128	result_geo = result_geo.reset_index()
   129	result_geo = result_geo[["hex_id", "rank", "area_weighted_density", "geometry"]]
   130	result_geo = result_geo.sort_values("rank")
   131	print(f"\nhotspots.geoparquet: {len(result_geo)} rows")
   132	
   133	# hotspot_ranking.parquet — no geometry
   134	result_plain = top[["rank", "area_weighted_density", "n_overlap_polygons", "n_slivers_filtered"]].copy()
   135	result_plain = result_plain.reset_index()
   136	result_plain = result_plain[["hex_id", "rank", "area_weighted_density", "n_overlap_polygons", "n_slivers_filtered"]]
   137	result_plain = result_plain.sort_values("rank")
   138	print(f"hotspot_ranking.parquet: {len(result_plain)} rows")
   139	
   140	# ── 11. Write outputs ─────────────────────────────────────────────────────────
   141	result_geo.to_parquet("/work/hotspots.geoparquet", index=False)
   142	print("  ✓ hotspots.geoparquet written")
   143	
   144	result_plain.to_parquet("/work/hotspot_ranking.parquet", index=False)
   145	print("  ✓ hotspot_ranking.parquet written")
   146	
   147	# ── 12. Sanity checks ─────────────────────────────────────────────────────────
   148	print("\n─── Sanity checks ───")
   149	# Check hex_ids match
   150	ids_geo = set(result_geo["hex_id"])
   151	ids_plain = set(result_plain["hex_id"])
   152	assert ids_geo == ids_plain, "hex_id sets differ!"
   153	print(f"✓ Both files share the same {len(ids_geo)} hex_id values")
   154	
   155	# Check rank is unique 1..n
   156	assert result_geo["rank"].is_monotonic_increasing, "ranks not sorted ascending"
   157	assert list(result_geo["rank"]) == list(range(1, len(result_geo) + 1)), "ranks not 1..n"
   158	print("✓ Ranks are unique 1..n")
   159	
   160	# Check ordering: descending density
   161	assert all(
   162	    result_geo["area_weighted_density"].iloc[i] >= result_geo["area_weighted_density"].iloc[i+1]
   163	    for i in range(len(result_geo) - 1)
   164	), "density not descending by rank"
   165	print("✓ Density is descending with rank")
   166	
   167	# Verify CRS of geo output
   168	print(f"✓ Output geometry CRS: {result_geo.crs}")
   169	
   170	# Print top 5 rows
   171	print("\nTop 5 ranked cells:")
   172	print(result_geo.head(5).to_string(index=False))
   173	print()
   174	print(result_plain.head(5).to_string(index=False))