fio-l2-capetown-landuse-dissolve
Map (reference ← swipe → agent)
0:00
I'm prepping a transit-corridor study for the Cape Town metro, and I need a class-level summary of `landuse` so the team's spatial-SQL notebooks can join it against the bus-route table. Can you give me one row per landuse class with the geometry unified into a single MultiPolygon, the total area in m², and the count of source parcels that fed in? Please write the result to `landuse_dissolved.geoparquet` with `class`, `area_m2`, and `parcel_count` columns.
0:00
System info
{
"adapter": "openrouter",
"model": "deepseek/deepseek-v4-flash",
"base_url": "https://openrouter.ai/api/v1"
}
0:04
Let me start by looking at the working directory to see what input files are available.
0:04
Bash
ls /work/
capetown_landuse.fgb
0:10
Let's first inspect the input file to understand its schema and content.
0:10
Bash
cd /work && python -c "
cd /work && python -c "
import geopandas as gpd
gdf = gpd.read_file('capetown_landuse.fgb')
print('Columns:', gdf.columns.tolist())
print('CRS:', gdf.crs)
print('Geometry types:', gdf.geometry.type.unique())
print('Shape:', gdf.shape)
print('---')
print(gdf.head(3))
print('---')
print(gdf['class'].value_counts())
"
Columns: ['id', 'class', 'subtype', 'geometry']
CRS: EPSG:32734
Geometry types: <ArrowStringArray>
['MultiPolygon']
Length: 1, dtype: str
Shape: (30978, 4)
---
id ... geometry
0 f328c054-b040-3a7d-a015-53747d608fef ... MULTIPOLYGON (((872000.605 6252622.417, 871205...
1 f1943847-f785-390c-abb0-edc923f97997 ... MULTIPOLYGON (((253543.223 6283900.241, 254921...
2 c37e40aa-f202-3b15-b68d-89dafeeaa478 ... MULTIPOLYGON (((255047.308 6280773.284, 255024...
[3 rows x 4 columns]
---
class
vineyard 8476
residential 3537
pitch 3393
farmland 2289
grass 2208
...
airfield 1
music_school 1
doctors 1
training_area 1
protected_landscape_seascape 1
Name: count, Length: 72, dtype: int64
0:30
The data is in CRS EPSG:32734 (UTM zone 34S, which covers Cape Town). That's a metric CRS, so area calculations will be in m². Good. The task: 1. Dissolve by `class` - union all geometries per class into a single MultiPolygon 2. Compute total area in m² per class 3. Count of source parcels per class 4. Output to `landuse_dissolved.geoparquet` with columns `class`, `area_m2`, `parcel_count` Let me write a script.
0:30
Write
/work/solve.py
content · 32 lines · py
import geopandas as gpd
# Read the source FlatGeobuf
gdf = gpd.read_file("/work/capetown_landuse.fgb")
# CRS is EPSG:32734 (UTM zone 34S) — metric, so area will be in m²
crs = gdf.crs
print(f"Input CRS: {crs}")
print(f"Input rows: {len(gdf)}")
# Dissolve by class: union geometries, count parcels
dissolved = gdf.dissolve(by="class", aggfunc="count")
# Rename / compute columns
# After dissolve:
# - index is 'class'
# - geometry is the dissolved (unioned) geometry
# - id (or any column) holds the count since we used aggfunc="count"
dissolved = dissolved.rename(columns={"id": "parcel_count"})
dissolved["area_m2"] = dissolved.geometry.area
# Keep only the columns we need
result = dissolved[["parcel_count", "area_m2", "geometry"]].copy()
result.index.name = "class"
result = result.reset_index()
print(f"Output classes: {len(result)}")
print(result.head(10))
# Write to GeoParquet
result.to_parquet("/work/landuse_dissolved.geoparquet", index=False)
print("✅ Written to /work/landuse_dissolved.geoparquet")Wrote /work/solve.py (1023 chars).
0:33
Let me run the script.
0:33
Bash
cd /work && python solve.py
Input CRS: EPSG:32734
Input rows: 30978
Output classes: 72
class ... geometry
0 airfield ... POLYGON ((267829.758 6244743.665, 267820.919 6...
1 allotments ... MULTIPOLYGON (((263870.512 6232477.907, 263865...
2 animal_keeping ... MULTIPOLYGON (((264114.249 6261265.896, 264220...
3 barracks ... MULTIPOLYGON (((262288.681 6215543.331, 262280...
4 base ... POLYGON ((272069.73 6245748.538, 272064.541 62...
5 beach_resort ... POLYGON ((292732.371 6227028.582, 292737.841 6...
6 brownfield ... MULTIPOLYGON (((275011.734 6227089.673, 275021...
7 bunker ... MULTIPOLYGON (((262135.155 6221253.116, 262135...
8 camp_site ... MULTIPOLYGON (((300525.478 6221802.335, 300450...
9 cemetery ... MULTIPOLYGON (((264302.444 6213143.812, 264327...
[10 rows x 4 columns]
✅ Written to /work/landuse_dissolved.geoparquet
0:48
Let me verify the output to make sure it has the right columns and is valid.
0:48
Bash
cd /work && python -c "
cd /work && python -c "
import geopandas as gpd
gdf = gpd.read_parquet('landuse_dissolved.geoparquet')
print('Columns:', gdf.columns.tolist())
print('CRS:', gdf.crs)
print('Rows:', len(gdf))
print('---')
print(gdf.dtypes)
print('---')
# Check all geometries are MultiPolygon
types = gdf.geometry.type.unique()
print('Geometry types:', types)
# Show a couple rows
print(gdf.head(5))
"
Columns: ['class', 'parcel_count', 'area_m2', 'geometry']
CRS: {"$schema": "https://proj.org/schemas/v0.7/projjson.schema.json", "type": "ProjectedCRS", "name": "WGS 84 / UTM zone 34S", "base_crs": {"name": "WGS 84", "datum_ensemble": {"name": "World Geodetic System 1984 ensemble", "members": [{"name": "World Geodetic System 1984 (Transit)"}, {"name": "World Geodetic System 1984 (G730)"}, {"name": "World Geodetic System 1984 (G873)"}, {"name": "World Geodetic System 1984 (G1150)"}, {"name": "World Geodetic System 1984 (G1674)"}, {"name": "World Geodetic System 1984 (G1762)"}, {"name": "World Geodetic System 1984 (G2139)"}, {"name": "World Geodetic System 1984 (G2296)"}], "ellipsoid": {"name": "WGS 84", "semi_major_axis": 6378137, "inverse_flattening": 298.257223563}, "accuracy": "2.0", "id": {"authority": "EPSG", "code": 6326}}, "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": 4326}}, "conversion": {"name": "UTM zone 34S", "method": {"name": "Transverse Mercator", "id": {"authority": "EPSG", "code": 9807}}, "parameters": [{"name": "Latitude of natural origin", "value": 0, "unit": "degree", "id": {"authority": "EPSG", "code": 8801}}, {"name": "Longitude of natural origin", "value": 21, "unit": "degree", "id": {"authority": "EPSG", "code": 8802}}, {"name": "Scale factor at natural origin", "value": 0.9996, "unit": "unity", "id": {"authority": "EPSG", "code": 8805}}, {"name": "False easting", "value": 500000, "unit": "metre", "id": {"authority": "EPSG", "code": 8806}}, {"name": "False northing", "value": 10000000, "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": "Navigation and medium accuracy spatial referencing.", "area": "Between 18\u00b0E and 24\u00b0E, southern hemisphere between 80\u00b0S and equator, onshore and offshore. Angola. Botswana. Democratic Republic of the Congo (Zaire). Namibia. South Africa. Zambia.", "bbox": {"south_latitude": -80, "west_longitude": 18, "north_latitude": 0, "east_longitude": 24}, "id": {"authority": "EPSG", "code": 32734}}
Rows: 72
---
class str
parcel_count int64
area_m2 float64
geometry geometry
dtype: object
---
Geometry types: <ArrowStringArray>
['Polygon', 'MultiPolygon']
Length: 2, dtype: str
class ... geometry
0 airfield ... POLYGON ((267829.758 6244743.665, 267820.919 6...
1 allotments ... MULTIPOLYGON (((263870.512 6232477.907, 263865...
2 animal_keeping ... MULTIPOLYGON (((264114.249 6261265.896, 264220...
3 barracks ... MULTIPOLYGON (((262288.681 6215543.331, 262280...
4 base ... POLYGON ((272069.73 6245748.538, 272064.541 62...
[5 rows x 4 columns]
0:56
The output looks good, but I notice some geometries are single Polygon rather than MultiPolygon. The user asked for "the geometry unified into a single MultiPolygon". Let me fix that to ensure all geometries are MultiPolygon.
0:56
0:56
Edit
/work/solve.py
Edited /work/solve.py.
0:59
Wait, the import is wrong. Let me fix it more cleanly.
0:59
Edit
/work/solve.py
Edited /work/solve.py.
1:04
Let me check the current state of the script and then simplify the geometry handling.
1:04
Read
/work/solve.py
1 import geopandas as gpd
2 from shapely.geometry import MultiPolygon
3
4 # Read the source FlatGeobuf
5 gdf = gpd.read_file("/work/capetown_landuse.fgb")
6
7 # CRS is EPSG:32734 (UTM zone 34S) — metric, so area will be in m²
8 crs = gdf.crs
9 print(f"Input CRS: {crs}")
10 print(f"Input rows: {len(gdf)}")
11
12 # Dissolve by class: union geometries, count parcels
13 dissolved = gdf.dissolve(by="class", aggfunc="count")
14
15 # Rename / compute columns
16 # After dissolve:
17 # - index is 'class'
18 # - geometry is the dissolved (unioned) geometry
19 # - id (or any column) holds the count since we used aggfunc="count"
20 dissolved = dissolved.rename(columns={"id": "parcel_count"})
21 dissolved["area_m2"] = dissolved.geometry.area
22
23 # Ensure every geometry is MultiPolygon (some singletons stay as Polygon)
24 import shapely
25 result = dissolved[["parcel_count", "area_m2", "geometry"]].copy()
26 result["geometry"] = result["geometry"].apply(
27 lambda g: shapely.geometry.collection.GeometryCollection([g]) if g is None else
28 shapely.geometry.MultiPolygon([g]) if g.geom_type == "Polygon" else g
29 )
30 result.index.name = "class"
31 result = result.reset_index()
32
33 print(f"Output classes: {len(result)}")
34 print(result.head(10))
35
36 # Write to GeoParquet
37 result.to_parquet("/work/landuse_dissolved.geoparquet", index=False)
38 print("✅ Written to /work/landuse_dissolved.geoparquet")
1:12
The code uses `shapely.geometry.MultiPolygon` directly and imports `shapely` inside the apply. Let me clean that up - already imported MultiPolygon at the top. Let me use it properly.
1:12
Write
/work/solve.py
content · 36 lines · py
import geopandas as gpd
from shapely.geometry import MultiPolygon
# Read the source FlatGeobuf
gdf = gpd.read_file("/work/capetown_landuse.fgb")
# CRS is EPSG:32734 (UTM zone 34S) — metric, so area will be in m²
crs = gdf.crs
print(f"Input CRS: {crs}")
print(f"Input rows: {len(gdf)}")
# Dissolve by class: union geometries, count parcels
dissolved = gdf.dissolve(by="class", aggfunc="count")
# Rename / compute columns
# After dissolve:
# - index is 'class'
# - geometry is the dissolved (unioned) geometry
# - id (or any column) holds the count since we used aggfunc="count"
dissolved = dissolved.rename(columns={"id": "parcel_count"})
dissolved["area_m2"] = dissolved.geometry.area
# Ensure every geometry is MultiPolygon (single Polygons stay as Polygon)
result = dissolved[["parcel_count", "area_m2", "geometry"]].copy()
result["geometry"] = result["geometry"].apply(
lambda g: MultiPolygon([g]) if g.geom_type == "Polygon" else g
)
result.index.name = "class"
result = result.reset_index()
print(f"Output classes: {len(result)}")
print(result.head(10))
# Write to GeoParquet
result.to_parquet("/work/landuse_dissolved.geoparquet", index=False)
print("✅ Written to /work/landuse_dissolved.geoparquet")Wrote /work/solve.py (1223 chars).
1:17
Let me run the script and verify the output.
1:17
Bash
cd /work && python solve.py
Input CRS: EPSG:32734
Input rows: 30978
Output classes: 72
class ... geometry
0 airfield ... MULTIPOLYGON (((267829.758 6244743.665, 267820...
1 allotments ... MULTIPOLYGON (((263870.512 6232477.907, 263865...
2 animal_keeping ... MULTIPOLYGON (((264114.249 6261265.896, 264220...
3 barracks ... MULTIPOLYGON (((262288.681 6215543.331, 262280...
4 base ... MULTIPOLYGON (((272069.73 6245748.538, 272064....
5 beach_resort ... MULTIPOLYGON (((292732.371 6227028.582, 292737...
6 brownfield ... MULTIPOLYGON (((275011.734 6227089.673, 275021...
7 bunker ... MULTIPOLYGON (((262135.155 6221253.116, 262135...
8 camp_site ... MULTIPOLYGON (((300525.478 6221802.335, 300450...
9 cemetery ... MULTIPOLYGON (((264302.444 6213143.812, 264327...
[10 rows x 4 columns]
✅ Written to /work/landuse_dissolved.geoparquet
1:34
Now all geometries are MultiPolygon. Let me do a final verification pass.
1:34
Bash
cd /work && python -c "
cd /work && python -c "
import geopandas as gpd
gdf = gpd.read_parquet('landuse_dissolved.geoparquet')
print('Columns:', gdf.columns.tolist())
print('Rows:', len(gdf))
print('Geometry types:', gdf.geometry.type.unique())
print('CRS datum:', gdf.crs.is_projected)
print('---')
# Sum area sanity check
print(f'Total area: {gdf[\"area_m2\"].sum():,.0f} m²')
print(f'Total parcels: {gdf[\"parcel_count\"].sum()}')
print(f'Source rows: 30978 => matches? {gdf[\"parcel_count\"].sum() == 30978}')
print('---')
# Quick look at a few rows
print(gdf[['class','parcel_count','area_m2']].head(10))
# ls output
"
ls -lh /work/landuse_dissolved.geoparquet
Columns: ['class', 'parcel_count', 'area_m2', 'geometry']
Rows: 72
Geometry types: <ArrowStringArray>
['MultiPolygon']
Length: 1, dtype: str
CRS datum: True
---
Total area: 3,876,566,473 m²
Total parcels: 30978
Source rows: 30978 => matches? True
---
class parcel_count area_m2
0 airfield 1 2.444056e+06
1 allotments 215 1.715652e+04
2 animal_keeping 302 6.636804e+05
3 barracks 4 1.873189e+03
4 base 1 9.244053e+05
5 beach_resort 1 2.936896e+05
6 brownfield 13 1.418610e+05
7 bunker 638 6.712463e+04
8 camp_site 13 8.568784e+05
9 cemetery 65 4.515795e+06
-rw-r--r-- 1 agent agent 8.0M Jun 18 14:24 /work/landuse_dissolved.geoparquet
1:43
Everything checks out: 1. **Columns**: `class`, `parcel_count`, `area_m2`, `geometry` — exactly what was requested. 2. **All geometries are MultiPolygon** — the single Polygons got wrapped. 3. **CRS** is EPSG:32734 (projected, metric) — area values are in m², magnitude correct (~10⁵–10⁶ for UTM coords). 4. **Parcel count sums to 30,978** — matches the source row count exactly. 5. **72 classes** — one row per landuse class. 6. **File exists** at `/work/landuse_dissolved.geoparquet`.
1:43