geo-l2-bangkok-landuse-intersect
Map (reference ← swipe → agent)
0:00
Working on a flood-mitigation green-cover briefing. `landcover` has land-cover polygons for the Bangkok metro; `study_area` defines the BMA boundary. I need the land-cover within the study area, simplified enough that the file is small for the policy lead to preview in a browser, with the original `class` string and a per-feature `area_m2` in square metres.
Write the result to `bma_landcover_intersect.geojson` as a GeoJSON FeatureCollection with every feature stored as a MultiPolygon.
0:00
System info
{
"adapter": "openrouter",
"model": "deepseek/deepseek-v4-flash",
"base_url": "https://openrouter.ai/api/v1"
}
0:03
Let me start by reading the input files to understand their structure, CRS, and contents.
0:03
Bash
ls -la /work/
total 13024 drwxr-xr-x 1 agent agent 94 Jun 18 01:21 . drwxr-xr-x 1 root root 14 Jun 18 01:21 .. -rw-r--r-- 1 agent agent 13328966 Jun 18 01:21 bangkok_landcover.parquet -rw-r--r-- 1 agent agent 1350 Jun 18 01:21 bma_study_area.geojson
0:11
Let me first inspect the input files to understand their structure.
0:11
Bash
cd /work && python -c "
cd /work && python -c "
import geopandas as gpd
# Read the study area
study = gpd.read_file('bma_study_area.geojson')
print('=== STUDY AREA ===')
print(study.head())
print('CRS:', study.crs)
print('Geometry types:', study.geometry.geom_type.unique())
print('---')
# Read landcover parquet
lc = gpd.read_parquet('bangkok_landcover.parquet')
print('=== LANDCOVER ===')
print(lc.head())
print('Columns:', lc.columns.tolist())
print('CRS:', lc.crs)
print('Shape:', lc.shape)
print('Geometry types:', lc.geometry.geom_type.unique())
print('Classes:', lc['class'].unique() if 'class' in lc.columns else 'N/A')
"
=== STUDY AREA ===
name geometry
0 BMA study area POLYGON ((653609.601 1506149.152, 662313.016 1...
CRS: EPSG:32647
Geometry types: <ArrowStringArray>
['Polygon']
Length: 1, dtype: str
---
=== LANDCOVER ===
id ... geometry
0 00023927-8760-5550-a47d-48ae7f087574 ... POLYGON ((647454.402 1506832.283, 647738.841 1...
1 0002ad8a-093a-539c-8f6c-0fb0b4f2c3c1 ... POLYGON ((643806.222 1536031.975, 643880.924 1...
2 001206dc-28bb-593e-87bd-72efd2f32c65 ... POLYGON ((695139.782 1532663.1, 695231.345 153...
3 001748f2-f294-53b7-be74-29acd35b3af6 ... POLYGON ((650585.19 1544580.727, 650660.692 15...
4 001ca192-b61e-56b8-a9b1-8f6d429f6fc0 ... POLYGON ((652703.522 1505402.689, 652740.294 1...
[5 rows x 3 columns]
Columns: ['id', 'class', 'geometry']
CRS: {"$schema": "https://proj.org/schemas/v0.7/projjson.schema.json", "type": "ProjectedCRS", "name": "WGS 84 / UTM zone 47N", "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 47N", "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": 99, "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": 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": "Navigation and medium accuracy spatial referencing.", "area": "Between 96\u00b0E and 102\u00b0E, northern hemisphere between equator and 84\u00b0N, onshore and offshore. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Myanmar (Burma). Russian Federation. Thailand.", "bbox": {"south_latitude": 0, "west_longitude": 96, "north_latitude": 84, "east_longitude": 102}, "id": {"authority": "EPSG", "code": 32647}}
Shape: (21660, 3)
Geometry types: <ArrowStringArray>
['Polygon']
Length: 1, dtype: str
Classes: <ArrowStringArray>
['crop', 'urban', 'forest', 'shrub', 'barren', 'wetland', 'grass', 'mangrove']
Length: 8, dtype: str
0:28
Good. Let me understand the data: 1. `bma_study_area.geojson` — 1 polygon, CRS EPSG:32647 (UTM zone 47N, metres) 2. `bangkok_landcover.parquet` — 21660 polygons, same CRS, with columns `id`, `class`, `geometry`. Classes: crop, urban, forest, shrub, barren, wetland, grass, mangrove. All Polygons. Tasks: 1. Intersect land-cover polygons with the study area boundary 2. Simplify enough to keep file small for browser preview 3. Keep original `class` string 4. Add `area_m2` in square metres (per feature) 5. Store every feature as MultiPolygon 6. Write to `bma_landcover_intersect.geojson` as GeoJSON FeatureCollection For simplification — the output needs to be small enough for browser preview. I'll use a reasonable simplification tolerance. Since the CRS is in metres (UTM 47N), something like 1m or 5m would be reasonable. Let me think... Actually, for a policy lead to preview, it just needs to not be enormous. The input is 21660 polygons. After intersection with the BMA boundary, it'll be fewer. I'll use a gentle simplification, maybe 1 metre tolerance, to reduce coordinate precision while keeping shape fidelity. Or maybe 0.5m. Let me use 1m — that's still very accurate for a metropolitan-scale map. Also need to ensure every geometry is a MultiPolygon (not Polygon). Let me also reduce coordinate precision. GeoJSON with lots of decimal places inflates file size. Since the CRS is in metres, 1cm precision (2 decimal places) is fine for this scale. Let me write the script.
0:28
Write
/work/solve.py
content · 56 lines · py
#!/usr/bin/env python3
"""Intersect landcover with BMA study area, simplify, compute area_m2, output GeoJSON."""
import geopandas as gpd
import pandas as pd
from shapely.geometry import MultiPolygon
# 1. Read inputs
study = gpd.read_file("/work/bma_study_area.geojson")
lc = gpd.read_parquet("/work/bangkok_landcover.parquet")
# Both are EPSG:32647 (metres) — confirm
assert study.crs.to_epsg() == 32647, f"Study CRS: {study.crs}"
assert lc.crs.to_epsg() == 32647, f"Landcover CRS: {lc.crs}"
study_area_geom = study.geometry.iloc[0]
# 2. Clip (spatial intersect) — faster than pairwise overlay
lc_clipped = lc.clip(study_area_geom, keep_geom_type=True)
print(f"After clip: {len(lc_clipped)} features")
# 3. Simplify (tolerance in metres — 1 m is negligible for a metro map)
# Use preserve_topology=True to avoid self-intersections
tol = 1.0
lc_clipped["geometry"] = lc_clipped.geometry.simplify(tol, preserve_topology=True)
# 4. Compute area in square metres (CRS is metric)
lc_clipped["area_m2"] = lc_clipped.geometry.area.round(1)
# 5. Ensure every geometry is MultiPolygon
def to_multipolygon(geom):
if geom is None:
return None
if geom.geom_type == "Polygon":
return MultiPolygon([geom])
if geom.geom_type == "MultiPolygon":
return geom
# Shouldn't happen with keep_geom_type=True, but handle gracefully
raise ValueError(f"Unexpected geometry type: {geom.geom_type}")
lc_clipped["geometry"] = lc_clipped.geometry.apply(to_multipolygon)
# 6. Keep only needed columns, drop id to keep file smaller
result = lc_clipped[["class", "area_m2", "geometry"]].copy()
# 7. Write GeoJSON — reduce coordinate precision to 2 decimals (cm-level)
# The GeoJSON driver writes with configurable precision
result.to_file(
"/work/bma_landcover_intersect.geojson",
driver="GeoJSON",
engine="pyogrio",
**{"COORDINATE_PRECISION": 2},
)
print("Done — wrote bma_landcover_intersect.geojson")Wrote /work/solve.py (1949 chars).
0:33
Let me run the script and check the output.
0:33
Bash
cd /work && python solve.py
Traceback (most recent call last):
File "/work/solve.py", line 19, in <module>
lc_clipped = lc.clip(study_area_geom, keep_geom_type=True)
File "/usr/local/lib/python3.14/site-packages/geopandas/geodataframe.py", line 2776, in clip
return geopandas.clip(self, mask=mask, keep_geom_type=keep_geom_type, sort=sort)
~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/local/lib/python3.14/site-packages/geopandas/tools/clip.py", line 218, in clip
clipped = _clip_gdf_with_mask(gdf, combined_mask, sort=sort)
File "/usr/local/lib/python3.14/site-packages/geopandas/tools/clip.py", line 93, in _clip_gdf_with_mask
gdf_sub.geometry.values[non_point_mask].intersection(mask)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^
File "/usr/local/lib/python3.14/site-packages/geopandas/array.py", line 858, in intersection
self._binary_method("intersection", self, other), crs=self.crs
~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/local/lib/python3.14/site-packages/geopandas/array.py", line 795, in _binary_method
return getattr(shapely, op)(left._data, right, **kwargs)
~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/local/lib/python3.14/site-packages/shapely/decorators.py", line 173, in wrapper
result = func(*args, **kwargs)
File "/usr/local/lib/python3.14/site-packages/shapely/decorators.py", line 88, in wrapped
return func(*args, **kwargs)
File "/usr/local/lib/python3.14/site-packages/shapely/set_operations.py", line 168, in intersection
return lib.intersection(a, b, **kwargs)
~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^
shapely.errors.GEOSException: TopologyException: side location conflict at 671567.85823834711 1501411.7438156258. This can occur if the input geometry is invalid.
0:39
The clip operation failed due to topology issues. Let me try with `keep_geom_type=False` or use an overlay approach instead. Let me first check if the study area geometry is valid.
0:39
Bash
cd /work && python -c "
cd /work && python -c "
import geopandas as gpd
study = gpd.read_file('bma_study_area.geojson')
print('Study area valid:', study.geometry.iloc[0].is_valid)
if not study.geometry.iloc[0].is_valid:
buf = study.geometry.iloc[0].buffer(0)
print('Buffered valid:', buf.is_valid)
study.at[0, 'geometry'] = buf
study.to_file('bma_study_area_fixed.geojson', driver='GeoJSON')
print('Fixed saved')
"
Study area valid: True
0:44
The study area is valid. The topology error is likely from the landcover data having invalid geometries during intersection. Let me check the landcover data and fix invalid geometries before clipping.
0:44
Bash
cd /work && python -c "
cd /work && python -c "
import geopandas as gpd
lc = gpd.read_parquet('bangkok_landcover.parquet')
print('Total:', len(lc))
print('Invalid:', (~lc.geometry.is_valid).sum())
# Check a few invalid
invalid = lc[~lc.geometry.is_valid]
if len(invalid) > 0:
print(invalid.head())
"
Total: 21660
Invalid: 25
id ... geometry
0 00023927-8760-5550-a47d-48ae7f087574 ... POLYGON ((647454.402 1506832.283, 647738.841 1...
866 0a425bd9-2aa5-5978-a12a-c0125b5d370e ... POLYGON ((698782.177 1496736.328, 699047.867 1...
1732 151cbf21-86ed-58f9-82ad-6cc6f74bf01d ... POLYGON ((643203.004 1530416.541, 643335.782 1...
2598 1fa2fc51-db29-51c9-a6e0-a3efc442fe12 ... POLYGON ((700485.954 1527549.694, 700754.336 1...
3464 2978c0d4-37c5-5b10-b8c4-6851e43f769a ... POLYGON ((650489.329 1502756.332, 650629.523 1...
[5 rows x 3 columns]
0:48
There are 25 invalid geometries in the landcover data. I should fix those with `buffer(0)` before the clip operation. Let me update the script.
0:48
Read
/work/solve.py
1 #!/usr/bin/env python3
2 """Intersect landcover with BMA study area, simplify, compute area_m2, output GeoJSON."""
3
4 import geopandas as gpd
5 import pandas as pd
6 from shapely.geometry import MultiPolygon
7
8 # 1. Read inputs
9 study = gpd.read_file("/work/bma_study_area.geojson")
10 lc = gpd.read_parquet("/work/bangkok_landcover.parquet")
11
12 # Both are EPSG:32647 (metres) — confirm
13 assert study.crs.to_epsg() == 32647, f"Study CRS: {study.crs}"
14 assert lc.crs.to_epsg() == 32647, f"Landcover CRS: {lc.crs}"
15
16 study_area_geom = study.geometry.iloc[0]
17
18 # 2. Clip (spatial intersect) — faster than pairwise overlay
19 lc_clipped = lc.clip(study_area_geom, keep_geom_type=True)
20
21 print(f"After clip: {len(lc_clipped)} features")
22
23 # 3. Simplify (tolerance in metres — 1 m is negligible for a metro map)
24 # Use preserve_topology=True to avoid self-intersections
25 tol = 1.0
26 lc_clipped["geometry"] = lc_clipped.geometry.simplify(tol, preserve_topology=True)
27
28 # 4. Compute area in square metres (CRS is metric)
29 lc_clipped["area_m2"] = lc_clipped.geometry.area.round(1)
30
31 # 5. Ensure every geometry is MultiPolygon
32 def to_multipolygon(geom):
33 if geom is None:
34 return None
35 if geom.geom_type == "Polygon":
36 return MultiPolygon([geom])
37 if geom.geom_type == "MultiPolygon":
38 return geom
39 # Shouldn't happen with keep_geom_type=True, but handle gracefully
40 raise ValueError(f"Unexpected geometry type: {geom.geom_type}")
41
42 lc_clipped["geometry"] = lc_clipped.geometry.apply(to_multipolygon)
43
44 # 6. Keep only needed columns, drop id to keep file smaller
45 result = lc_clipped[["class", "area_m2", "geometry"]].copy()
46
47 # 7. Write GeoJSON — reduce coordinate precision to 2 decimals (cm-level)
48 # The GeoJSON driver writes with configurable precision
49 result.to_file(
50 "/work/bma_landcover_intersect.geojson",
51 driver="GeoJSON",
52 engine="pyogrio",
53 **{"COORDINATE_PRECISION": 2},
54 )
55
56 print("Done — wrote bma_landcover_intersect.geojson")
0:57
The clip failed because of invalid geometries in the landcover data. I need to fix invalid geometries first, then clip. Let me use `buffer(0)` to fix them, and also catch any topology issues during intersection. I'll do a manual spatial intersection using overlay or a more robust approach.
0:57
Write
/work/solve.py
content · 63 lines · py
#!/usr/bin/env python3
"""Intersect landcover with BMA study area, simplify, compute area_m2, output GeoJSON."""
import geopandas as gpd
from shapely.geometry import MultiPolygon
# 1. Read inputs
study = gpd.read_file("/work/bma_study_area.geojson")
lc = gpd.read_parquet("/work/bangkok_landcover.parquet")
# Both are EPSG:32647 (metres) — confirm
assert study.crs.to_epsg() == 32647
assert lc.crs.to_epsg() == 32647
study_area_geom = study.geometry.iloc[0]
# 2. Fix any invalid geometries in landcover (buffer(0))
bad = ~lc.geometry.is_valid
if bad.any():
print(f"Fixing {bad.sum()} invalid geometries...")
lc.loc[bad, "geometry"] = lc.loc[bad, "geometry"].buffer(0)
# 3. Spatial subset first using sjoin (bbox filter) to make clip more robust
# Then use overlay='intersection' which handles topology better
bbox = study_area_geom.envelope
lc_in_bbox = lc[lc.intersects(bbox)].copy()
print(f"Features in study bbox: {len(lc_in_bbox)}")
# 4. Intersect via overlay (more robust than clip)
lc_inter = gpd.overlay(lc_in_bbox, study, how="intersection", keep_geom_type=True)
print(f"After intersection: {len(lc_inter)} features")
# 5. Simplify (tolerance in metres)
tol = 1.0
lc_inter["geometry"] = lc_inter.geometry.simplify(tol, preserve_topology=True)
# 6. Compute area in square metres (CRS is metric)
lc_inter["area_m2"] = lc_inter.geometry.area.round(1)
# 7. Ensure every geometry is MultiPolygon
def to_multipolygon(geom):
if geom is None:
return None
if geom.geom_type == "Polygon":
return MultiPolygon([geom])
if geom.geom_type == "MultiPolygon":
return geom
raise ValueError(f"Unexpected geometry type: {geom.geom_type}")
lc_inter["geometry"] = lc_inter.geometry.apply(to_multipolygon)
# 8. Keep only needed columns, drop id to keep file smaller
result = lc_inter[["class", "area_m2", "geometry"]].copy()
# 9. Write GeoJSON — reduce coordinate precision to 2 decimals (cm-level)
result.to_file(
"/work/bma_landcover_intersect.geojson",
driver="GeoJSON",
engine="pyogrio",
**{"COORDINATE_PRECISION": 2},
)
print("Done — wrote bma_landcover_intersect.geojson")Wrote /work/solve.py (2145 chars).
1:02
Wait, `overlay` with a single polygon study area might also have topology issues. Let me think about a more robust approach. The issue is that some landcover polygons are invalid and when doing intersection with the study area, the invalid geometries cause topology exceptions. Actually, I already fixed the invalid geometries. Let me also consider that the `overlay` function might still have issues. Let me try a different approach: use `intersection` directly on each geometry after fixing, which is what `clip` does but with pre-fixed geometries. Actually, let me also consider that the issue could be with the study area being very large and complex, and the clip/overlay operations having numerical precision issues at polygon boundaries. A more robust approach: 1. Fix invalid geometries in the landcover dataset 2. Use a row-by-row intersection with `shapely.intersection` wrapped in try/except 3. Or use the `clip` function after fixing both datasets Let me try a simpler approach - just fix invalid geometries first and then clip.
1:02
Bash
cd /work && python solve.py
Fixing 25 invalid geometries... Features in study bbox: 6470 After intersection: 3453 features Done — wrote bma_landcover_intersect.geojson
1:06
The script ran successfully. Let me now verify the output file.
1:06
Bash
ls -lh /work/bma_landcover_intersect.geojson
-rw-r--r-- 1 agent agent 3.2M Jun 18 01:22 /work/bma_landcover_intersect.geojson
1:10
The file is 3.2MB. Let me verify the content is as expected - check geometry types, columns, etc.
1:10
Bash
cd /work && python -c "
cd /work && python -c "
import geopandas as gpd
result = gpd.read_file('bma_landcover_intersect.geojson')
print('CRS:', result.crs)
print('Columns:', result.columns.tolist())
print('Shape:', result.shape)
print('Geometry types:', result.geometry.geom_type.unique())
print('Classes:', result['class'].unique())
print('Area range:', result['area_m2'].min(), '-', result['area_m2'].max())
print('Total area (m2):', result['area_m2'].sum())
print()
print(result.head())
"
CRS: EPSG:32647
Columns: ['class', 'area_m2', 'geometry']
Shape: (3453, 3)
Geometry types: <ArrowStringArray>
['MultiPolygon']
Length: 1, dtype: str
Classes: <ArrowStringArray>
['forest', 'barren', 'crop', 'shrub', 'urban', 'wetland', 'mangrove', 'grass']
Length: 8, dtype: str
Area range: 0.2 - 628183762.9
Total area (m2): 979954270.0
class area_m2 geometry
0 forest 15735.7 MULTIPOLYGON (((660072.92 1501476.79, 660051.3...
1 forest 32598.5 MULTIPOLYGON (((662162.93 1503314.93, 662117.9...
2 barren 17000.2 MULTIPOLYGON (((687414.87 1509653.4, 687417.38...
3 crop 106792.7 MULTIPOLYGON (((688024.32 1521450.52, 688025.6...
4 forest 18200.7 MULTIPOLYGON (((684185.57 1508479.55, 684215.4...
1:18
Everything looks good. Let me verify the key requirements: 1. ✅ Land-cover within the study area (clip/intersection was done) 2. ✅ Simplified enough for small file (3.2M with simplification tolerance 1m and precision 2) 3. ✅ Original `class` string preserved 4. ✅ Per-feature `area_m2` in square metres 5. ✅ Output file named `bma_landcover_intersect.geojson` 6. ✅ All geometries are MultiPolygon 7. ✅ GeoJSON FeatureCollection All checks pass.
1:18