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 exploring the working directory to find any relevant data files, particularly landuse data.
0:04
Bash
ls /work/
capetown_landuse.fgb
0:08
Let me explore the data file to understand its structure.
0:08
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import geopandas as gpd
gdf = gpd.read_file('capetown_landuse.fgb')
print('Columns:', gdf.columns.tolist())
print('CRS:', gdf.crs)
print('Shape:', gdf.shape)
print('dtypes:')
print(gdf.dtypes)
print()
print('First rows:')
print(gdf.head())
"
Columns: ['id', 'class', 'subtype', 'geometry']
CRS: EPSG:32734
Shape: (30978, 4)
dtypes:
id str
class str
subtype str
geometry geometry
dtype: object
First rows:
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 b9d39db7-5dcc-379e-bfeb-c994628df5a7 ... MULTIPOLYGON (((254978.723 6280842.179, 254992...
4 0033302c-c1d6-3624-b71e-0e300db64963 ... MULTIPOLYGON (((264759.63 6280552.85, 264825.5...
[5 rows x 4 columns]
0:16
Let me understand the task: 1. Read the `capetown_landuse.fgb` file 2. Group by `class` column 3. For each class: - Dissolve (unify) all geometries into a single MultiPolygon - Calculate the area in m² (since CRS is EPSG:32734 which is UTM 34S with meters) - Count the number of source parcels 4. Write to `landuse_dissolved.geoparquet` with columns: `class`, `area_m2`, `parcel_count` Let me write a Python script to do this.
0:16
Write
/work/solve.py
content · 36 lines · py
import geopandas as gpd
import pandas as pd
from shapely.geometry import MultiPolygon
# Read the source file
gdf = gpd.read_file('/work/capetown_landuse.fgb')
print(f"Total parcels: {len(gdf)}")
print(f"Unique classes: {gdf['class'].nunique()}")
print(gdf['class'].value_counts())
# Group by 'class', dissolve geometries, count parcels
dissolved = gdf.dissolve(by='class', aggfunc='count')
# Rename the id count column to parcel_count
dissolved = dissolved.rename(columns={'id': 'parcel_count'})
# Only keep the columns we need
dissolved = dissolved[['parcel_count', 'geometry']]
# Calculate area in m² (CRS is EPSG:32734, units are metres)
dissolved['area_m2'] = dissolved.geometry.area
# Reset index so 'class' becomes a column
dissolved = dissolved.reset_index()
# Keep only desired columns in the right order
result = dissolved[['class', 'area_m2', 'parcel_count', 'geometry']]
print("\nResult summary:")
print(result[['class', 'area_m2', 'parcel_count']].to_string())
# Write to GeoParquet
result.to_parquet('/work/landuse_dissolved.geoparquet')
print("\nWritten to /work/landuse_dissolved.geoparquet")Wrote /work/solve.py (1118 chars).
0:27
Let me run the script.
0:27
Bash
cd /work && python3 solve.py
Total parcels: 30978
Unique classes: 72
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
Result summary:
class area_m2 parcel_count
0 airfield 2.444056e+06 1
1 allotments 1.715652e+04 215
2 animal_keeping 6.636804e+05 302
3 barracks 1.873189e+03 4
4 base 9.244053e+05 1
5 beach_resort 2.936896e+05 1
6 brownfield 1.418610e+05 13
7 bunker 6.712463e+04 638
8 camp_site 8.568784e+05 13
9 cemetery 4.515795e+06 65
10 clinic 9.552029e+04 19
11 college 1.178636e+07 37
12 commercial 7.190500e+06 327
13 construction 5.727714e+06 84
14 doctors 1.407718e+02 1
15 dog_park 1.021826e+05 11
16 driving_range 1.287389e+05 10
17 education 3.349939e+05 8
18 fairway 1.040348e+06 232
19 farmland 3.817290e+08 2289
20 farmyard 7.748543e+06 294
21 flowerbed 6.537666e+03 40
22 garages 6.635553e+02 3
23 garden 4.671889e+06 731
24 golf_course 1.353686e+07 50
25 grass 3.260907e+07 2208
26 grave_yard 1.920088e+05 17
27 green 1.316131e+05 310
28 greenfield 3.279033e+05 14
29 greenhouse_horticulture 8.850429e+04 5
30 hospital 4.673730e+06 55
31 industrial 6.496317e+07 321
32 kindergarten 1.540748e+05 65
33 landfill 4.550490e+06 18
34 lateral_water_hazard 6.541641e+04 12
35 marina 1.557197e+05 5
36 meadow 1.785986e+07 397
37 military 5.032478e+06 22
38 music_school 1.887399e+04 1
39 national_park 4.566999e+08 3
40 nature_reserve 1.700276e+09 42
41 naval_base 4.584897e+05 4
42 orchard 2.764034e+07 1052
43 park 2.179358e+07 1980
44 pedestrian 4.111786e+05 266
45 pitch 1.133306e+07 3393
46 plant_nursery 8.625895e+05 40
47 playground 5.245591e+05 506
48 plaza 2.755790e+04 3
49 protected_landscape_seascape 2.586956e+06 1
50 quarry 5.748884e+06 30
51 railway 4.382439e+06 19
52 recreation_ground 2.230580e+06 80
53 religious 1.284939e+06 358
54 residential 3.613121e+08 3537
55 resort 9.421589e+05 6
56 retail 8.692102e+06 414
57 rough 8.837839e+03 36
58 school 2.882922e+07 1046
59 species_management_area 1.758947e+07 3
60 stadium 4.357762e+05 12
61 tee 6.753886e+04 576
62 track 3.204643e+05 35
63 training_area 5.419677e+03 1
64 university 4.222409e+06 23
65 village_green 5.769054e+04 45
66 vineyard 3.329834e+08 8476
67 water_hazard 5.217910e+05 69
68 water_park 4.026974e+04 5
69 wilderness_area 3.064923e+08 5
70 works 2.503782e+06 68
71 zoo 4.314377e+05 5
Written to /work/landuse_dissolved.geoparquet
0:37
Let me verify the output file exists and is valid.
0:37
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import geopandas as gpd
result = gpd.read_parquet('landuse_dissolved.geoparquet')
print('Columns:', result.columns.tolist())
print('CRS:', result.crs)
print('Shape:', result.shape)
print()
print(result.head())
print()
# Verify all geometries are MultiPolygon
from shapely.geometry import MultiPolygon
all_mp = result.geometry.map(lambda g: isinstance(g, MultiPolygon)).all()
print(f'All MultiPolygons: {all_mp}')
"
Columns: ['class', 'area_m2', 'parcel_count', '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}}
Shape: (72, 4)
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]
All MultiPolygons: False
0:56
The output shows that some geometries are plain Polygon instead of MultiPolygon. The user asked for "geometry unified into a single MultiPolygon". So I should ensure all geometries are MultiPolygon. Let me update the script to convert any Polygon geometries to MultiPolygon.
0:56
0:56
Read
/work/solve.py
1 import geopandas as gpd
2 import pandas as pd
3 from shapely.geometry import MultiPolygon
4
5 # Read the source file
6 gdf = gpd.read_file('/work/capetown_landuse.fgb')
7
8 print(f"Total parcels: {len(gdf)}")
9 print(f"Unique classes: {gdf['class'].nunique()}")
10 print(gdf['class'].value_counts())
11
12 # Group by 'class', dissolve geometries, count parcels
13 dissolved = gdf.dissolve(by='class', aggfunc='count')
14
15 # Rename the id count column to parcel_count
16 dissolved = dissolved.rename(columns={'id': 'parcel_count'})
17
18 # Only keep the columns we need
19 dissolved = dissolved[['parcel_count', 'geometry']]
20
21 # Calculate area in m² (CRS is EPSG:32734, units are metres)
22 dissolved['area_m2'] = dissolved.geometry.area
23
24 # Reset index so 'class' becomes a column
25 dissolved = dissolved.reset_index()
26
27 # Keep only desired columns in the right order
28 result = dissolved[['class', 'area_m2', 'parcel_count', 'geometry']]
29
30 print("\nResult summary:")
31 print(result[['class', 'area_m2', 'parcel_count']].to_string())
32
33 # Write to GeoParquet
34 result.to_parquet('/work/landuse_dissolved.geoparquet')
35
36 print("\nWritten to /work/landuse_dissolved.geoparquet")
1:02
I need to enforce that all geometries are MultiPolygon. Some dissolve results for single-parcel classes (like `airfield`, `base`, etc.) come out as plain `Polygon`. I should convert them using `shapely.geometry.MultiPolygon` to wrap them.
1:02
Edit
/work/solve.py
Edited /work/solve.py.
1:04
Let me re-run the script to ensure all geometries are MultiPolygons.
1:04
Bash
cd /work && python3 solve.py
Total parcels: 30978
Unique classes: 72
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
Result summary:
class area_m2 parcel_count
0 airfield 2.444056e+06 1
1 allotments 1.715652e+04 215
2 animal_keeping 6.636804e+05 302
3 barracks 1.873189e+03 4
4 base 9.244053e+05 1
5 beach_resort 2.936896e+05 1
6 brownfield 1.418610e+05 13
7 bunker 6.712463e+04 638
8 camp_site 8.568784e+05 13
9 cemetery 4.515795e+06 65
10 clinic 9.552029e+04 19
11 college 1.178636e+07 37
12 commercial 7.190500e+06 327
13 construction 5.727714e+06 84
14 doctors 1.407718e+02 1
15 dog_park 1.021826e+05 11
16 driving_range 1.287389e+05 10
17 education 3.349939e+05 8
18 fairway 1.040348e+06 232
19 farmland 3.817290e+08 2289
20 farmyard 7.748543e+06 294
21 flowerbed 6.537666e+03 40
22 garages 6.635553e+02 3
23 garden 4.671889e+06 731
24 golf_course 1.353686e+07 50
25 grass 3.260907e+07 2208
26 grave_yard 1.920088e+05 17
27 green 1.316131e+05 310
28 greenfield 3.279033e+05 14
29 greenhouse_horticulture 8.850429e+04 5
30 hospital 4.673730e+06 55
31 industrial 6.496317e+07 321
32 kindergarten 1.540748e+05 65
33 landfill 4.550490e+06 18
34 lateral_water_hazard 6.541641e+04 12
35 marina 1.557197e+05 5
36 meadow 1.785986e+07 397
37 military 5.032478e+06 22
38 music_school 1.887399e+04 1
39 national_park 4.566999e+08 3
40 nature_reserve 1.700276e+09 42
41 naval_base 4.584897e+05 4
42 orchard 2.764034e+07 1052
43 park 2.179358e+07 1980
44 pedestrian 4.111786e+05 266
45 pitch 1.133306e+07 3393
46 plant_nursery 8.625895e+05 40
47 playground 5.245591e+05 506
48 plaza 2.755790e+04 3
49 protected_landscape_seascape 2.586956e+06 1
50 quarry 5.748884e+06 30
51 railway 4.382439e+06 19
52 recreation_ground 2.230580e+06 80
53 religious 1.284939e+06 358
54 residential 3.613121e+08 3537
55 resort 9.421589e+05 6
56 retail 8.692102e+06 414
57 rough 8.837839e+03 36
58 school 2.882922e+07 1046
59 species_management_area 1.758947e+07 3
60 stadium 4.357762e+05 12
61 tee 6.753886e+04 576
62 track 3.204643e+05 35
63 training_area 5.419677e+03 1
64 university 4.222409e+06 23
65 village_green 5.769054e+04 45
66 vineyard 3.329834e+08 8476
67 water_hazard 5.217910e+05 69
68 water_park 4.026974e+04 5
69 wilderness_area 3.064923e+08 5
70 works 2.503782e+06 68
71 zoo 4.314377e+05 5
Written to /work/landuse_dissolved.geoparquet
1:18
Let me verify that now all geometries are MultiPolygon.
1:18
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import geopandas as gpd
from shapely.geometry import MultiPolygon
result = gpd.read_parquet('landuse_dissolved.geoparquet')
print('Columns:', result.columns.tolist())
print('CRS:', result.crs)
print('Shape:', result.shape)
all_mp = result.geometry.map(lambda g: isinstance(g, MultiPolygon)).all()
print(f'All MultiPolygon geometries: {all_mp}')
# Quick summary stats
print()
print('Top 5 by area:')
print(result.nlargest(5, 'area_m2')[['class', 'area_m2', 'parcel_count']].to_string())
print()
print('Top 5 by parcel count:')
print(result.nlargest(5, 'parcel_count')[['class', 'area_m2', 'parcel_count']].to_string())
"
Columns: ['class', 'area_m2', 'parcel_count', '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}}
Shape: (72, 4)
All MultiPolygon geometries: True
Top 5 by area:
class area_m2 parcel_count
40 nature_reserve 1.700276e+09 42
39 national_park 4.566999e+08 3
19 farmland 3.817290e+08 2289
54 residential 3.613121e+08 3537
66 vineyard 3.329834e+08 8476
Top 5 by parcel count:
class area_m2 parcel_count
66 vineyard 3.329834e+08 8476
54 residential 3.613121e+08 3537
45 pitch 1.133306e+07 3393
19 farmland 3.817290e+08 2289
25 grass 3.260907e+07 2208
1:25
Everything looks good. The file `landuse_dissolved.geoparquet` is written with: - 72 rows (one per landuse class) - `class` column (string) - `area_m2` column (float, in square metres, calculated from the projected CRS EPSG:32734) - `parcel_count` column (integer, count of source parcels) - `geometry` column (all MultiPolygon type) - CRS preserved as EPSG:32734 Let me confirm the file exists.
1:25