spa-l2-lagos-hotspot-overlaps
Map (reference ← swipe → agent)
0:00
Hot-spot map for the Lagos housing brief. `lagos_landuse` carries land-use polygons with a population density attribute; `lagos_hex_grid` is a hex grid. I need the top 10 % of hex cells ranked by area-weighted mean population density across overlapping land-use polygons. Tiny polygons under 100 m² are noise and should be excluded.
I need two output files that share exactly the same set of `hex_id` values, one row per cell. `hotspots.geoparquet` is a GeoParquet (Polygon or MultiPolygon geometry) in Nigeria's national grid with columns `hex_id`, `rank`, `area_weighted_density`, and `geometry`. `hotspot_ranking.parquet` is a plain Parquet with no geometry, carrying `hex_id`, `rank`, `area_weighted_density`, `n_overlap_polygons` (integer), and `n_slivers_filtered` (integer). In both files `rank` must be a unique integer starting at 1 where ascending rank corresponds to descending `area_weighted_density`.
0:00
System info
{
"adapter": "openrouter",
"model": "deepseek/deepseek-v4-flash",
"base_url": "https://openrouter.ai/api/v1"
}
0:04
Let me start by understanding the task: 1. I have two input files: `lagos_landuse` (land-use polygons with population density attribute) and `lagos_hex_grid` (hex grid). 2. I need to find the top 10% of hex cells ranked by area-weighted mean population density. 3. Small polygons under 100 m² should be excluded (noise/slivers). 4. Outputs: - `hotspots.geoparquet` - GeoParquet with Polygon/MultiPolygon geometry, columns: hex_id, rank, area_weighted_density, geometry - in Nigeria's national grid CRS - `hotspot_ranking.parquet` - plain Parquet with no geometry, columns: hex_id, rank, area_weighted_density, n_overlap_polygons (int), n_slivers_filtered (int) 5. Both files share the same hex_id values (only top 10%) 6. rank is unique integer starting at 1, descending area_weighted_density Let me first explore the input files to understand their structure, attributes, and CRS.
0:04
0:04
Bash
ls -la /work/
total 4244 drwxr-xr-x 1 agent agent 86 Jun 18 02:00 . drwxr-xr-x 1 root root 14 Jun 18 02:00 .. -rw-r--r-- 1 agent agent 736774 Jun 18 02:00 lagos_hex_grid.geojson -rw-r--r-- 1 agent agent 3607119 Jun 18 02:00 lagos_landuse.geojson
0:08
Let me explore both files to understand their structure.
0:08
Bash
python3 -c "
python3 -c "
import geopandas as gpd
import pandas as pd
# Read landuse
lu = gpd.read_file('/work/lagos_landuse.geojson')
print('=== LANDUSE ===')
print('Columns:', lu.columns.tolist())
print('CRS:', lu.crs)
print('Shape:', lu.shape)
print('Dtypes:')
print(lu.dtypes)
print()
print('Area range (m²):', lu.area.min(), '-', lu.area.max())
print(lu.head(3))
print()
# Read hex grid
hg = gpd.read_file('/work/lagos_hex_grid.geojson')
print('=== HEX GRID ===')
print('Columns:', hg.columns.tolist())
print('CRS:', hg.crs)
print('Shape:', hg.shape)
print('Dtypes:')
print(hg.dtypes)
print(hg.head(3))
"
=== LANDUSE ===
Columns: ['id', 'class', 'pop_density', 'geometry']
CRS: EPSG:4326
Shape: (5542, 4)
Dtypes:
id str
class str
pop_density float64
geometry geometry
dtype: object
Area range (m²): 8.184750909128568e-11 - 0.002854502247539594
id ... geometry
0 004ca249-d2cf-3617-906e-dea1138e762b ... POLYGON ((3.46992 6.44615, 3.46993 6.44566, 3....
1 006a7dc6-2b5e-34fa-84a2-8a679620c9de ... POLYGON ((3.42136 6.43091, 3.42118 6.43051, 3....
2 00a334d1-7e8e-3e4c-8565-3501d4ed74df ... POLYGON ((3.45824 6.43361, 3.45872 6.43306, 3....
[3 rows x 4 columns]
=== HEX GRID ===
Columns: ['hex_id', 'geometry']
CRS: EPSG:4326
Shape: (1782, 2)
Dtypes:
hex_id str
geometry geometry
dtype: object
hex_id geometry
0 H000-000 POLYGON ((3.25491 6.40082, 3.2523 6.40535, 3.2...
1 H000-001 POLYGON ((3.25491 6.40987, 3.2523 6.41439, 3.2...
2 H000-002 POLYGON ((3.25492 6.41891, 3.25231 6.42344, 3....
<string>:14: UserWarning: Geometry is in a geographic CRS. Results from 'area' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.
0:18
Okay, I can see the data: 1. Landuse: 5542 polygons with id, class, pop_density. CRS: EPSG:4326 (geographic). The areas in degrees² are tiny because it's geographic. Need to reproject to Nigeria's national grid for area calculations. 2. Hex grid: 1782 hex cells with hex_id. CRS: EPSG:4326. Let me find out what Nigeria's national grid CRS is. The most common one is EPSG:26391 (Minna / Nigeria West Belt) or EPSG:26392 (Minna / Nigeria East Belt). But for all of Lagos, Nigeria West Belt (EPSG:26391) is appropriate. Let me check - actually, Lagos falls in the UTM zone 31N. Let me check which CRS is Nigeria's national grid. Nigeria's national grid is based on the "Minna" datum. The two belt projections: - EPSG:26391 - Minna / Nigeria West Belt (zone 2) - covers Lagos area - EPSG:26392 - Minna / Nigeria East Belt Lagos is at around longitude 3.4°E, so it's in the West Belt (EPSG:26391). Actually, another option is EPSG:26331 (Minna / UTM zone 31N) or EPSG:26332 (Minna / UTM zone 32N). Let me search for "Nigeria national grid CRS" or just check the CRS on the datasets. Let me check the area of the landuse polygons to understand what "100 m²" means. Actually, the areas shown are in degrees squared (since CRS is 4326). Let me reproject to a suitable projected CRS and check areas. Let me first figure out which projected CRS to use. EPSG:26391 (Minna / Nigeria West Belt) covers Lagos.
0:18
Bash
python3 -c "
python3 -c "
import geopandas as gpd
# Check what projected CRS would be suitable
lu = gpd.read_file('/work/lagos_landuse.geojson')
print('Bounds:', lu.total_bounds)
# Try EPSG:26391 - Minna / Nigeria West Belt
lu_proj = lu.to_crs('EPSG:26391')
print('Proj bounds:', lu_proj.total_bounds)
print('Area range (m²):', lu_proj.area.min(), '-', lu_proj.area.max())
print('Total area:', lu_proj.area.sum())
# Count polygons under 100 m²
small_count = (lu_proj.area < 100).sum()
print(f'Polygons < 100 m²: {small_count} out of {len(lu_proj)}')
"
Bounds: [3.25026042 6.4009372 3.6701081 6.6980539 ] Proj bounds: [ 92614.05688451 265487.763265 139052.06293866 298364.62844424] Area range (m²): 1.0007638007186197 - 34909797.88823449 Total area: 515536457.3676117 Polygons < 100 m²: 3022 out of 5542
0:44
Interesting - most polygons are under 100 m². Let me check the distribution more carefully. The area in EPSG:4326 was tiny (as expected) but now in EPSG:26391 (Minna / Nigeria West Belt) we get more meaningful areas. 3022 out of 5542 polygons are under 100 m². That's a lot of "noise" polygons to filter out. Let me now plan the approach: 1. Load both datasets 2. Filter landuse polygons with area >= 100 m² (need to project first) 3. For each hex cell, find overlapping landuse polygons 4. Calculate area-weighted mean pop_density 5. Sort, assign ranks, take top 10% 6. Output both files with Nigeria's national grid CRS Wait, "Nigeria's national grid" - what CRS is that? Let me check. Nigeria uses several Transverse Mercator projections - typically "Minna / Nigeria West Belt" (EPSG:26391) for the western part (where Lagos is) or "Minna / Nigeria Mid Belt" (EPSG:26392) or "Minna / Nigeria East Belt" (EPSG:26393). Since Lagos is in the west, EPSG:26391 is the right one. Actually, let me think about the "Nigeria's national grid" more carefully. The "Nigeria National Grid" could also refer to EPSG:26381 (Minna / Nigeria West) or simply EPSG:26391. Let me just use EPSG:26391 which is Minna / Nigeria West Belt covering Lagos. Let me write the full script. For area-weighted mean population density: - For each hex cell, find all landuse polygons that intersect with it - For each overlapping pair, compute the intersection area - area_weighted_density = sum(pop_density_i * intersection_area_i) / sum(intersection_area_i) "n_overlap_polygons" = number of landuse polygons that overlap with this hex cell (after filtering slivers) "n_slivers_filtered" = number of landuse polygons that were excluded because they are under 100 m² Wait, "n_slivers_filtered" - this is per hex cell. How many small polygons (< 100 m²) were filtered out that overlap with this hex cell? Let me re-read: "Tiny polygons under 100 m² are noise and should be excluded." And the column `n_slivers_filtered` (integer) - I think it means for each hex cell, the number of landuse polygons that overlap with it but were filtered out because they're under 100 m². Let me be more careful. Let me count per hex cell: - All landuse polygons that overlap the hex - Among those, how many are < 100 m² (filtered/slivers) - Among those, how many are >= 100 m² (used for the calculation) Actually, let me think again. "Tiny polygons under 100 m² are noise and should be excluded." This means we filter out all polygons < 100 m² from the dataset before the overlap calculation. Then `n_slivers_filtered` is, for each hex cell, the count of polygons < 100 m² that overlapped the hex cell. And `n_overlap_polygons` is the count of polygons >= 100 m² that overlap the hex cell. Let me write the script. Now for the top 10% - we have 1782 hex cells, so top 10% = 178.2, ceil to 178? Or floor? Let's say ceil = 179. Actually, top 10% of hex cells... Let me use ceil (round up) to be safe. So ceil(1782 * 0.1) = ceil(178.2) = 179. Wait, let me check: "the top 10 % of hex cells" - typically this means the top 10% when ranked by density. So if there are 1782 cells, we take the top 179 cells (rounding up). Actually, looking at it more carefully - "the top 10 % of hex cells ranked by area-weighted mean population density". So we rank all hex cells by area_weighted_density descending, and take the top 10% (ceiling). Let me code this up.
0:44
Write
/work/solve.py
content · 143 lines · py
import geopandas as gpd
import pandas as pd
import numpy as np
# ------------------------------------------------------------
# 1. Load data
# ------------------------------------------------------------
lu = gpd.read_file('/work/lagos_landuse.geojson')
hg = gpd.read_file('/work/lagos_hex_grid.geojson')
# ------------------------------------------------------------
# 2. Project to Nigeria national grid (Minna / Nigeria West Belt, EPSG:26391)
# Lagos falls in the west belt.
# ------------------------------------------------------------
NIGERIA_GRID = 'EPSG:26391'
lu_proj = lu.to_crs(NIGERIA_GRID)
hg_proj = hg.to_crs(NIGERIA_GRID)
# ------------------------------------------------------------
# 3. Compute polygon areas and split into "valid" (>= 100 m²) and "slivers" (< 100 m²)
# ------------------------------------------------------------
lu_proj['poly_area'] = lu_proj.area # area in m²
lu_valid = lu_proj[lu_proj['poly_area'] >= 100].copy()
lu_sliver = lu_proj[lu_proj['poly_area'] < 100].copy()
print(f"Valid polygons (>= 100 m²): {len(lu_valid)}")
print(f"Sliver polygons (< 100 m²): {len(lu_sliver)}")
# ------------------------------------------------------------
# 4. Spatial join: hex grid × valid landuse polygons
# For each hex cell we need:
# - intersecting valid polygons (for area-weighted density)
# - count of intersecting slivers (n_slivers_filtered)
# ------------------------------------------------------------
# 4a. Valid overlaps
valid_join = gpd.sjoin(
hg_proj[['hex_id', 'geometry']],
lu_valid[['pop_density', 'geometry', 'poly_area']],
how='inner',
predicate='intersects'
)
# Compute intersection area for each hex-polygon pair
# We loop over rows or use overlay
print("Computing intersection areas for valid overlaps...")
# Use overlay 'intersection' to get exact intersection geometries
overlay_valid = gpd.overlay(
hg_proj[['hex_id', 'geometry']],
lu_valid[['pop_density', 'geometry']],
how='intersection'
)
overlay_valid['intersect_area'] = overlay_valid.area
# Area-weighted density per hex cell
agg_valid = (
overlay_valid
.groupby('hex_id')
.apply(
lambda grp: pd.Series({
'area_weighted_density': (
(grp['pop_density'] * grp['intersect_area']).sum()
/ grp['intersect_area'].sum()
),
'n_overlap_polygons': len(grp)
})
)
.reset_index()
)
print(f"Hex cells with valid overlaps: {len(agg_valid)}")
# 4b. Sliver overlaps (count only)
sliver_join = gpd.sjoin(
hg_proj[['hex_id', 'geometry']],
lu_sliver[['geometry']],
how='inner',
predicate='intersects'
)
n_slivers = (
sliver_join
.groupby('hex_id')
.size()
.reset_index(name='n_slivers_filtered')
)
# ------------------------------------------------------------
# 5. Merge results
# ------------------------------------------------------------
results = hg_proj[['hex_id', 'geometry']].merge(
agg_valid, on='hex_id', how='left'
).merge(
n_slivers, on='hex_id', how='left'
)
# Fill NA: hex cells with no valid overlaps get 0 density, etc.
results['area_weighted_density'] = results['area_weighted_density'].fillna(0.0)
results['n_overlap_polygons'] = results['n_overlap_polygons'].fillna(0).astype(int)
results['n_slivers_filtered'] = results['n_slivers_filtered'].fillna(0).astype(int)
# ------------------------------------------------------------
# 6. Rank all hex cells (descending density)
# ------------------------------------------------------------
results = results.sort_values('area_weighted_density', ascending=False).reset_index(drop=True)
results['rank'] = range(1, len(results) + 1)
# ------------------------------------------------------------
# 7. Top 10 % (ceiling)
# ------------------------------------------------------------
n_top = int(np.ceil(len(results) * 0.10))
hotspots = results.head(n_top).copy()
print(f"Total hex cells: {len(results)}")
print(f"Top 10% (ceiling): {n_top}")
# ------------------------------------------------------------
# 8. Output files
# hotspots.geoparquet: geometry in Nigeria grid, cols: hex_id, rank, area_weighted_density, geometry
# hotspot_ranking.parquet: no geometry, cols: hex_id, rank, area_weighted_density, n_overlap_polygons, n_slivers_filtered
# ------------------------------------------------------------
# 8a. GeoParquet
hotspots_out = hotspots[['hex_id', 'rank', 'area_weighted_density', 'geometry']].copy()
hotspots_out.to_parquet(
'/work/hotspots.geoparquet',
index=False,
schema_version='1.0.0',
version='2.6'
)
# 8b. Plain Parquet (no geometry)
rank_out = hotspots[['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered']].copy()
rank_out.to_parquet(
'/work/hotspot_ranking.parquet',
index=False,
version='2.6'
)
print("Done. Files written.")
print("\nhotspots.geoparquet columns:", hotspots_out.columns.tolist())
print("hotspot_ranking.parquet columns:", rank_out.columns.tolist())
# Verify CRS
print(f"\nhotspots CRS: {hotspots_out.crs}")Wrote /work/solve.py (5076 chars).
1:03
Wait, I used `.apply` with a lambda on a grouped object which can be problematic. Let me use a more robust approach. Also, `overlay_valid` needs the polygon area from the original valid polygons since the overlay will carry attributes. Let me rethink the approach. Actually, let me think about this more carefully. For the area-weighted mean density per hex cell: - For each hex cell, we find all valid land-use polygons that intersect it. - For each intersecting pair, we compute the intersection area. - The area-weighted density = sum(pop_density_i * intersection_area_i) / sum(intersection_area_i) The issue with my approach using `gpd.overlay` is that after intersection, each resulting polygon has the attributes from both sides. But I also need the original polygon area to compute... wait, no. I need the intersection area, not the original polygon area. So `gpd.overlay` with `how='intersection'` works fine - the `area` of the resulting geometry is exactly the intersection area. But there's a subtlety: if a polygon partially overlaps a hex, only the part inside the hex contributes. That's exactly what we want. But wait - the `.apply` approach can be slow. Let me simplify with groupby aggregations. Actually, `groupby.apply` returning a Series works fine in modern pandas. Let me fix the approach to be cleaner. Let me also think about edge cases: 1. What if a hex cell has no valid overlapping polygons? Then area_weighted_density = 0, n_overlap_polygons = 0. 2. The ranking is across ALL hex cells, and we take top 10%. Let me rewrite with cleaner code.
1:03
Write
content · 128 lines
import geopandas as gpd
import pandas as pd
import numpy as np
# ------------------------------------------------------------
# 1. Load data
# ------------------------------------------------------------
print("Loading data...")
lu = gpd.read_file('/work/lagos_landuse.geojson')
hg = gpd.read_file('/work/lagos_hex_grid.geojson')
# ------------------------------------------------------------
# 2. Project to Nigeria national grid (Minna / Nigeria West Belt, EPSG:26391)
# ------------------------------------------------------------
NIGERIA_GRID = 'EPSG:26391'
lu_proj = lu.to_crs(NIGERIA_GRID)
hg_proj = hg.to_crs(NIGERIA_GRID)
# ------------------------------------------------------------
# 3. Compute polygon areas and split into valid vs slivers
# ------------------------------------------------------------
lu_proj['poly_area'] = lu_proj.area
lu_valid = lu_proj[lu_proj['poly_area'] >= 100].copy()
lu_sliver = lu_proj[lu_proj['poly_area'] < 100].copy()
print(f"Valid polygons (>= 100 m²): {len(lu_valid)}")
print(f"Sliver polygons (< 100 m²): {len(lu_sliver)}")
# ------------------------------------------------------------
# 4. Valid overlaps — area-weighted mean density per hex cell
# ------------------------------------------------------------
print("Overlaying hex grid with valid polygons...")
overlay = gpd.overlay(
hg_proj[['hex_id', 'geometry']],
lu_valid[['pop_density', 'geometry']],
how='intersection'
)
overlay['intersect_area'] = overlay.area
# Per hex cell: weighted sum / total intersect area
agg = overlay.groupby('hex_id').agg(
weighted_sum=('area_weighted_density', lambda s: (s * overlay.loc[s.index, 'intersect_area']).sum()),
# Actually easier: compute directly
)
# Better approach: compute numerator and denominator separately
# Weighted sum of density * intersect_area
num = (overlay['pop_density'] * overlay['intersect_area']).groupby(overlay['hex_id']).sum()
den = overlay['intersect_area'].groupby(overlay['hex_id']).sum()
area_weighted_density = num / den
# Number of distinct valid polygons overlapping
n_overlap = overlay.groupby('hex_id').size()
valid_stats = pd.DataFrame({
'area_weighted_density': area_weighted_density,
'n_overlap_polygons': n_overlap
}).reset_index()
print(f"Hex cells with valid overlaps: {len(valid_stats)}")
# ------------------------------------------------------------
# 5. Sliver overlaps — count only
# ------------------------------------------------------------
print("Counting sliver overlaps...")
sliver_join = gpd.sjoin(
hg_proj[['hex_id', 'geometry']],
lu_sliver[['geometry']],
how='inner',
predicate='intersects'
)
n_slivers = (
sliver_join
.groupby('hex_id')
.size()
.reset_index(name='n_slivers_filtered')
)
print(f"Hex cells with sliver overlaps: {len(n_slivers)}")
# ------------------------------------------------------------
# 6. Merge everything to full hex grid
# ------------------------------------------------------------
results = hg_proj[['hex_id', 'geometry']].copy()
results = results.merge(valid_stats, on='hex_id', how='left')
results = results.merge(n_slivers, on='hex_id', how='left')
results['area_weighted_density'] = results['area_weighted_density'].fillna(0.0)
results['n_overlap_polygons'] = results['n_overlap_polygons'].fillna(0).astype(int)
results['n_slivers_filtered'] = results['n_slivers_filtered'].fillna(0).astype(int)
# ------------------------------------------------------------
# 7. Rank all hex cells by descending area_weighted_density
# ------------------------------------------------------------
results = results.sort_values('area_weighted_density', ascending=False).reset_index(drop=True)
results['rank'] = range(1, len(results) + 1)
# ------------------------------------------------------------
# 8. Top 10 % (ceiling)
# ------------------------------------------------------------
n_top = int(np.ceil(len(results) * 0.10))
hotspots = results.head(n_top).copy().reset_index(drop=True)
print(f"\nTotal hex cells: {len(results)}")
print(f"Top 10% (ceiling of {len(results)*0.10:.1f}): {n_top} cells")
print(f"Rank range: {hotspots['rank'].min()} – {hotspots['rank'].max()}")
# ------------------------------------------------------------
# 9. Write output
# ------------------------------------------------------------
# 9a. GeoParquet — Polygon geometry, Nigeria grid
hotspots_geo = hotspots[['hex_id', 'rank', 'area_weighted_density', 'geometry']].copy()
hotspots_geo.to_parquet('/work/hotspots.geoparquet', index=False, schema_version='1.0.0')
# 9b. Plain Parquet — no geometry
hotspots_tab = hotspots[['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered']].copy()
hotspots_tab.to_parquet('/work/hotspot_ranking.parquet', index=False)
print("\n=== Output files ===")
print(f"hotspots.geoparquet — CRS: {hotspots_geo.crs}, cols: {hotspots_geo.columns.tolist()}, rows: {len(hotspots_geo)}")
print(f"hotspot_ranking.parquet — cols: {hotspots_tab.columns.tolist()}, rows: {len(hotspots_tab)}")
# Verify identical hex_id sets
s1 = set(hotspots_geo['hex_id'])
s2 = set(hotspots_tab['hex_id'])
assert s1 == s2, "hex_id sets differ!"
print("✓ Both files share identical hex_id sets")missing or non-string argument: 'file_path'
1:19
I need to provide a file_path parameter. Let me write this to /work/solve.py.
1:19
Write
/work/solve.py
content · 121 lines · py
import geopandas as gpd
import pandas as pd
import numpy as np
# ------------------------------------------------------------
# 1. Load data
# ------------------------------------------------------------
print("Loading data...")
lu = gpd.read_file('/work/lagos_landuse.geojson')
hg = gpd.read_file('/work/lagos_hex_grid.geojson')
# ------------------------------------------------------------
# 2. Project to Nigeria national grid (Minna / Nigeria West Belt, EPSG:26391)
# ------------------------------------------------------------
NIGERIA_GRID = 'EPSG:26391'
lu_proj = lu.to_crs(NIGERIA_GRID)
hg_proj = hg.to_crs(NIGERIA_GRID)
# ------------------------------------------------------------
# 3. Compute polygon areas and split into valid vs slivers
# ------------------------------------------------------------
lu_proj['poly_area'] = lu_proj.area
lu_valid = lu_proj[lu_proj['poly_area'] >= 100].copy()
lu_sliver = lu_proj[lu_proj['poly_area'] < 100].copy()
print(f"Valid polygons (>= 100 m²): {len(lu_valid)}")
print(f"Sliver polygons (< 100 m²): {len(lu_sliver)}")
# ------------------------------------------------------------
# 4. Valid overlaps — area-weighted mean density per hex cell
# ------------------------------------------------------------
print("Overlaying hex grid with valid polygons...")
overlay = gpd.overlay(
hg_proj[['hex_id', 'geometry']],
lu_valid[['pop_density', 'geometry']],
how='intersection'
)
overlay['intersect_area'] = overlay.area
# Weighted sum of density * intersect_area
num = (overlay['pop_density'] * overlay['intersect_area']).groupby(overlay['hex_id']).sum()
den = overlay['intersect_area'].groupby(overlay['hex_id']).sum()
area_weighted_density = num / den
# Number of distinct valid polygons overlapping
n_overlap = overlay.groupby('hex_id').size()
valid_stats = pd.DataFrame({
'area_weighted_density': area_weighted_density,
'n_overlap_polygons': n_overlap
}).reset_index()
print(f"Hex cells with valid overlaps: {len(valid_stats)}")
# ------------------------------------------------------------
# 5. Sliver overlaps — count only
# ------------------------------------------------------------
print("Counting sliver overlaps...")
sliver_join = gpd.sjoin(
hg_proj[['hex_id', 'geometry']],
lu_sliver[['geometry']],
how='inner',
predicate='intersects'
)
n_slivers = (
sliver_join
.groupby('hex_id')
.size()
.reset_index(name='n_slivers_filtered')
)
print(f"Hex cells with sliver overlaps: {len(n_slivers)}")
# ------------------------------------------------------------
# 6. Merge everything to full hex grid
# ------------------------------------------------------------
results = hg_proj[['hex_id', 'geometry']].copy()
results = results.merge(valid_stats, on='hex_id', how='left')
results = results.merge(n_slivers, on='hex_id', how='left')
results['area_weighted_density'] = results['area_weighted_density'].fillna(0.0)
results['n_overlap_polygons'] = results['n_overlap_polygons'].fillna(0).astype(int)
results['n_slivers_filtered'] = results['n_slivers_filtered'].fillna(0).astype(int)
# ------------------------------------------------------------
# 7. Rank all hex cells by descending area_weighted_density
# ------------------------------------------------------------
results = results.sort_values('area_weighted_density', ascending=False).reset_index(drop=True)
results['rank'] = range(1, len(results) + 1)
# ------------------------------------------------------------
# 8. Top 10 % (ceiling)
# ------------------------------------------------------------
n_top = int(np.ceil(len(results) * 0.10))
hotspots = results.head(n_top).copy().reset_index(drop=True)
print(f"\nTotal hex cells: {len(results)}")
print(f"Top 10% (ceiling of {len(results)*0.10:.1f}): {n_top} cells")
print(f"Rank range: {hotspots['rank'].min()} – {hotspots['rank'].max()}")
# ------------------------------------------------------------
# 9. Write output
# ------------------------------------------------------------
# 9a. GeoParquet — Polygon geometry, Nigeria grid
hotspots_geo = hotspots[['hex_id', 'rank', 'area_weighted_density', 'geometry']].copy()
hotspots_geo.to_parquet('/work/hotspots.geoparquet', index=False, schema_version='1.0.0')
# 9b. Plain Parquet — no geometry
hotspots_tab = hotspots[['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered']].copy()
hotspots_tab.to_parquet('/work/hotspot_ranking.parquet', index=False)
print("\n=== Output files ===")
print(f"hotspots.geoparquet — CRS: {hotspots_geo.crs}, cols: {hotspots_geo.columns.tolist()}, rows: {len(hotspots_geo)}")
print(f"hotspot_ranking.parquet — cols: {hotspots_tab.columns.tolist()}, rows: {len(hotspots_tab)}")
# Verify identical hex_id sets
s1 = set(hotspots_geo['hex_id'])
s2 = set(hotspots_tab['hex_id'])
assert s1 == s2, "hex_id sets differ!"
print("✓ Both files share identical hex_id sets")Wrote /work/solve.py (4915 chars).
1:23
Let me run the script now to see if it works.
1:23
Bash
cd /work && python solve.py 2>&1
Loading data... Valid polygons (>= 100 m²): 2520 Sliver polygons (< 100 m²): 3022 Overlaying hex grid with valid polygons... Hex cells with valid overlaps: 1032 Counting sliver overlaps... Hex cells with sliver overlaps: 1436 Total hex cells: 1782 Top 10% (ceiling of 178.2): 179 cells Rank range: 1 – 179 === Output files === hotspots.geoparquet — CRS: EPSG:26391, cols: ['hex_id', 'rank', 'area_weighted_density', 'geometry'], rows: 179 hotspot_ranking.parquet — cols: ['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered'], rows: 179 ✓ Both files share identical hex_id sets
1:31
The script ran successfully. Let me verify the outputs by reading them back and checking all the specifications.
1:31
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import pandas as pd
import geopandas as gpd
# Check hotspots.geoparquet
h = gpd.read_parquet('/work/hotspots.geoparquet')
print('=== hotspots.geoparquet ===')
print('Shape:', h.shape)
print('Columns:', h.columns.tolist())
print('CRS:', h.crs)
print('dtypes:')
print(h.dtypes)
print('head:')
print(h.head(10))
print()
# Check hotspot_ranking.parquet
r = pd.read_parquet('/work/hotspot_ranking.parquet')
print('=== hotspot_ranking.parquet ===')
print('Shape:', r.shape)
print('Columns:', r.columns.tolist())
print('dtypes:')
print(r.dtypes)
print('head:')
print(r.head(10))
print()
# Verify sets match
s1 = set(h['hex_id'])
s2 = set(r['hex_id'])
print(f'hex_id sets match: {s1 == s2}')
print(f'Count of unique hex_id in geo: {len(s1)}, in tab: {len(s2)}')
# Verify ranks 1..N unique
print(f'Rank range: {r[\"rank\"].min()} – {r[\"rank\"].max()}')
print(f'Ranks are unique: {r[\"rank\"].is_unique}')
print(f'Rank order is ascending density-descending: {list(r[\"area_weighted_density\"].is_monotonic_decreasing)}')
"
=== hotspots.geoparquet ===
Shape: (179, 4)
Columns: ['hex_id', 'rank', 'area_weighted_density', 'geometry']
CRS: {"$schema": "https://proj.org/schemas/v0.7/projjson.schema.json", "type": "ProjectedCRS", "name": "Minna / Nigeria West Belt", "base_crs": {"name": "Minna", "datum": {"type": "GeodeticReferenceFrame", "name": "Minna", "ellipsoid": {"name": "Clarke 1880 (RGS)", "semi_major_axis": 6378249.145, "inverse_flattening": 293.465}}, "coordinate_system": {"subtype": "ellipsoidal", "axis": [{"name": "Geodetic latitude", "abbreviation": "Lat", "direction": "north", "unit": "degree"}, {"name": "Geodetic longitude", "abbreviation": "Lon", "direction": "east", "unit": "degree"}]}, "id": {"authority": "EPSG", "code": 4263}}, "conversion": {"name": "Nigeria West Belt", "method": {"name": "Transverse Mercator", "id": {"authority": "EPSG", "code": 9807}}, "parameters": [{"name": "Latitude of natural origin", "value": 4, "unit": "degree", "id": {"authority": "EPSG", "code": 8801}}, {"name": "Longitude of natural origin", "value": 4.5, "unit": "degree", "id": {"authority": "EPSG", "code": 8802}}, {"name": "Scale factor at natural origin", "value": 0.99975, "unit": "unity", "id": {"authority": "EPSG", "code": 8805}}, {"name": "False easting", "value": 230738.26, "unit": "metre", "id": {"authority": "EPSG", "code": 8806}}, {"name": "False northing", "value": 0, "unit": "metre", "id": {"authority": "EPSG", "code": 8807}}]}, "coordinate_system": {"subtype": "Cartesian", "axis": [{"name": "Easting", "abbreviation": "E", "direction": "east", "unit": "metre"}, {"name": "Northing", "abbreviation": "N", "direction": "north", "unit": "metre"}]}, "scope": "Engineering survey, topographic mapping.", "area": "Nigeria - onshore west of 6\u00b030'E, onshore and offshore shelf.", "bbox": {"south_latitude": 3.57, "west_longitude": 2.69, "north_latitude": 13.9, "east_longitude": 6.5}, "id": {"authority": "EPSG", "code": 26391}}
dtypes:
hex_id str
rank int64
area_weighted_density float64
geometry geometry
dtype: object
head:
hex_id ... geometry
0 H011-017 ... POLYGON ((102683.536 282990.674, 102396.262 28...
1 H049-008 ... POLYGON ((135574.415 273890.473, 135287.161 27...
2 H031-006 ... POLYGON ((119976.505 271936.003, 119689.228 27...
3 H035-007 ... POLYGON ((123444.43 272926.044, 123157.16 2734...
4 H039-014 ... POLYGON ((126930.109 279917.444, 126642.857 28...
5 H039-016 ... POLYGON ((126936.075 281917.907, 126648.827 28...
6 H038-016 ... POLYGON ((126068.355 281420.376, 125781.104 28...
7 H038-015 ... POLYGON ((126065.372 280420.142, 125778.12 280...
8 H039-015 ... POLYGON ((126933.09 280917.676, 126645.84 2814...
9 H031-005 ... POLYGON ((119973.561 270935.743, 119686.282 27...
[10 rows x 4 columns]
=== hotspot_ranking.parquet ===
Shape: (179, 5)
Columns: ['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered']
dtypes:
hex_id str
rank int64
area_weighted_density float64
n_overlap_polygons int64
n_slivers_filtered int64
dtype: object
head:
hex_id rank ... n_overlap_polygons n_slivers_filtered
0 H011-017 1 ... 2 1
1 H049-008 2 ... 1 2
2 H031-006 3 ... 1 3
3 H035-007 4 ... 1 2
4 H039-014 5 ... 1 0
5 H039-016 6 ... 1 1
6 H038-016 7 ... 1 1
7 H038-015 8 ... 3 5
8 H039-015 9 ... 2 1
9 H031-005 10 ... 2 2
[10 rows x 5 columns]
hex_id sets match: True
Count of unique hex_id in geo: 179, in tab: 179
Rank range: 1 – 179
Ranks are unique: True
Traceback (most recent call last):
File "<string>", line 37, in <module>
print(f'Rank order is ascending density-descending: {list(r["area_weighted_density"].is_monotonic_decreasing)}')
~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: 'bool' object is not iterable
1:43
The script ran successfully. Let me fix the verification script quickly and also check that the CRS is correct. The CRS is PROJJSON for EPSG:26391 - Minna / Nigeria West Belt, which is Nigeria's national grid. Let me also verify the coordinate magnitudes make sense (should be ~10^5 meters for metric CRS). The coordinate values are like `102683.536, 282990.674` which are ~100k-300k meters - that's correct for EPSG:26391 in the Lagos area. Everything looks correct. Let me just fix the verification script's bug and verify one more time.
1:43
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import pandas as pd
import geopandas as gpd
h = gpd.read_parquet('/work/hotspots.geoparquet')
r = pd.read_parquet('/work/hotspot_ranking.parquet')
print('=== hotspots.geoparquet ===')
print('Shape:', h.shape)
print('Columns:', h.columns.tolist())
print('CRS:', h.crs)
print('dtypes:', dict(h.dtypes))
print()
print('=== hotspot_ranking.parquet ===')
print('Shape:', r.shape)
print('Columns:', r.columns.tolist())
print('dtypes:', dict(r.dtypes))
print()
# hex_id sets match
s1 = set(h['hex_id'])
s2 = set(r['hex_id'])
print(f'hex_id sets match: {s1 == s2}')
print()
# rank check
print(f'Rank range: {r[\"rank\"].min()} - {r[\"rank\"].max()}')
print(f'Ranks unique: {r[\"rank\"].is_unique}')
# check density is descending
d = r['area_weighted_density'].values
print(f'Density descending: {all(d[i] >= d[i+1] for i in range(len(d)-1))}')
print()
# coordinate magnitude check
bounds = h.total_bounds
print(f'X range: {bounds[0]:.1f} - {bounds[2]:.1f}')
print(f'Y range: {bounds[1]:.1f} - {bounds[3]:.1f}')
"
=== hotspots.geoparquet ===
Shape: (179, 4)
Columns: ['hex_id', 'rank', 'area_weighted_density', 'geometry']
CRS: {"$schema": "https://proj.org/schemas/v0.7/projjson.schema.json", "type": "ProjectedCRS", "name": "Minna / Nigeria West Belt", "base_crs": {"name": "Minna", "datum": {"type": "GeodeticReferenceFrame", "name": "Minna", "ellipsoid": {"name": "Clarke 1880 (RGS)", "semi_major_axis": 6378249.145, "inverse_flattening": 293.465}}, "coordinate_system": {"subtype": "ellipsoidal", "axis": [{"name": "Geodetic latitude", "abbreviation": "Lat", "direction": "north", "unit": "degree"}, {"name": "Geodetic longitude", "abbreviation": "Lon", "direction": "east", "unit": "degree"}]}, "id": {"authority": "EPSG", "code": 4263}}, "conversion": {"name": "Nigeria West Belt", "method": {"name": "Transverse Mercator", "id": {"authority": "EPSG", "code": 9807}}, "parameters": [{"name": "Latitude of natural origin", "value": 4, "unit": "degree", "id": {"authority": "EPSG", "code": 8801}}, {"name": "Longitude of natural origin", "value": 4.5, "unit": "degree", "id": {"authority": "EPSG", "code": 8802}}, {"name": "Scale factor at natural origin", "value": 0.99975, "unit": "unity", "id": {"authority": "EPSG", "code": 8805}}, {"name": "False easting", "value": 230738.26, "unit": "metre", "id": {"authority": "EPSG", "code": 8806}}, {"name": "False northing", "value": 0, "unit": "metre", "id": {"authority": "EPSG", "code": 8807}}]}, "coordinate_system": {"subtype": "Cartesian", "axis": [{"name": "Easting", "abbreviation": "E", "direction": "east", "unit": "metre"}, {"name": "Northing", "abbreviation": "N", "direction": "north", "unit": "metre"}]}, "scope": "Engineering survey, topographic mapping.", "area": "Nigeria - onshore west of 6\u00b030'E, onshore and offshore shelf.", "bbox": {"south_latitude": 3.57, "west_longitude": 2.69, "north_latitude": 13.9, "east_longitude": 6.5}, "id": {"authority": "EPSG", "code": 26391}}
dtypes: {'hex_id': <StringDtype(na_value=nan)>, 'rank': dtype('int64'), 'area_weighted_density': dtype('float64'), 'geometry': <geopandas.array.GeometryDtype object at 0x7f8955e281a0>}
=== hotspot_ranking.parquet ===
Shape: (179, 5)
Columns: ['hex_id', 'rank', 'area_weighted_density', 'n_overlap_polygons', 'n_slivers_filtered']
dtypes: {'hex_id': <StringDtype(na_value=nan)>, 'rank': dtype('int64'), 'area_weighted_density': dtype('float64'), 'n_overlap_polygons': dtype('int64'), 'n_slivers_filtered': dtype('int64')}
hex_id sets match: True
Rank range: 1 - 179
Ranks unique: True
Density descending: True
X range: 92829.5 - 139039.2
Y range: 264967.9 - 293983.8
1:51
Everything checks out: 1. Both files have exactly 179 rows (top 10% ceiling of 1782 = 179). 2. Both files share identical hex_id sets. 3. `hotspots.geoparquet` has columns: `hex_id`, `rank`, `area_weighted_density`, `geometry` — in EPSG:26391 (Nigeria West Belt / Minna). 4. `hotspot_ranking.parquet` has columns: `hex_id`, `rank`, `area_weighted_density`, `n_overlap_polygons`, `n_slivers_filtered` — no geometry. 5. Rank is unique, starting at 1, ascending with descending density. 6. Coordinate magnitudes are in meters (90k-140k E, 265k-294k N) which is correct for EPSG:26391 Nigeria West Belt. 7. Density is monotonically decreasing. All good! Let me also verify the files exist on disk.
1:51
1:51
Bash
ls -lh /work/hotspots.geoparquet /work/hotspot_ranking.parquet
-rw-r--r-- 1 agent agent 7.2K Jun 18 02:01 /work/hotspot_ranking.parquet -rw-r--r-- 1 agent agent 29K Jun 18 02:01 /work/hotspots.geoparquet
1:58
Everything looks good. Both output files are ready and verified.
1:58