geo-l2-bangkok-landuse-intersect

Rerun

-rw-r--r-- 1 agent agent 1350 Jun  9 12:22 bma_study_area.geojson
total 13024
drwxr-xr-x 1 agent agent       94 Jun  9 12:22 .
drwxr-xr-x 1 root  root        14 Jun  9 12:22 ..
-rw-r--r-- 1 agent agent 13328966 Jun  9 12:22 bangkok_landcover.parquet
-rw-r--r-- 1 agent agent     1350 Jun  9 12:22 bma_study_area.geojson

cd /work && python3 -c "
import geopandas as gpd
# Read 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.type.unique())
print()
"

=== 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

cd /work && python3 -c "
import pandas as pd
df = pd.read_parquet('bangkok_landcover.parquet')
print('=== LANDCOVER (parquet) ===')
print('Shape:', df.shape)
print('Columns:', df.columns.tolist())
print('CRS:', df.attrs.get('crs', 'N/A'))
print(df.head(2))
print()
print('class values:', df['class'].value_counts().to_dict())
print()
# Check geometry type
print('Geometry types:', df.geometry.type.unique())
# Check if any Multi vs Single
print('MultiPolygon:', (df.geometry.type == 'MultiPolygon').sum())
print('Polygon:', (df.geometry.type == 'Polygon').sum())
"

=== LANDCOVER (parquet) ===
Shape: (21660, 3)
Columns: ['id', 'class', 'geometry']
CRS: N/A
                                     id  ...                                           geometry
0  00023927-8760-5550-a47d-48ae7f087574  ...  b'\x01\x03\x00\x00\x00\x01\x00\x00\x00\x05\x00...
1  0002ad8a-093a-539c-8f6c-0fb0b4f2c3c1  ...  b'\x01\x03\x00\x00\x00\x01\x00\x00\x00\r\x00\x...

[2 rows x 3 columns]

class values: {'barren': 5482, 'shrub': 5062, 'forest': 4509, 'crop': 3267, 'urban': 1697, 'wetland': 1433, 'mangrove': 115, 'grass': 95}

Traceback (most recent call last):
  File "<string>", line 13, in <module>
    print('Geometry types:', df.geometry.type.unique())
                             ^^^^^^^^^^^^^^^^
  File "/usr/local/lib/python3.14/site-packages/pandas/core/generic.py", line 6206, in __getattr__
    return object.__getattribute__(self, name)
           ~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^
AttributeError: 'Series' object has no attribute 'type'. Did you mean: 'dtype'?

cd /work && python3 -c "
import pandas as pd
import geopandas as gpd
from shapely import wkb

df = pd.read_parquet('bangkok_landcover.parquet')

# Convert WKB to geometry
df['geometry'] = df['geometry'].apply(wkb.loads)
gdf = gpd.GeoDataFrame(df, geometry='geometry', crs='EPSG:32647')

print('Shape:', gdf.shape)
print('CRS:', gdf.crs)
print('Geometry types:', gdf.geometry.type.unique())
print()
print('MultiPolygon count:', (gdf.geometry.type == 'MultiPolygon').sum())
print('Polygon count:', (gdf.geometry.type == 'Polygon').sum())
print('Other:', (gdf.geometry.type != 'MultiPolygon') & (gdf.geometry.type != 'Polygon').sum())
print()
print(gdf.head(3))
"

Shape: (21660, 3)
CRS: EPSG:32647
Geometry types: <ArrowStringArray>
['Polygon']
Length: 1, dtype: str

MultiPolygon count: 0
Polygon count: 21660
Other: 0        False
1        False
2        False
3        False
4        False
         ...  
21655    False
21656    False
21657    False
21658    False
21659    False
Length: 21660, dtype: bool

                                     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 rows x 3 columns]

cd /work && python3 -c "
import pandas as pd
import geopandas as gpd
from shapely import wkb

df = pd.read_parquet('bangkok_landcover.parquet')
df['geometry'] = df['geometry'].apply(wkb.loads)
landcover = gpd.GeoDataFrame(df, geometry='geometry', crs='EPSG:32647')

study = gpd.read_file('bma_study_area.geojson')

# Clip/intersect
intersected = gpd.overlay(landcover, study, how='intersection')
print('After intersection:', intersected.shape)
print(intersected.head())
print('Memory:', intersected.memory_usage(deep=True).sum() / 1e6, 'MB')
"

After intersection: (3453, 4)
                                     id  ...                                           geometry
0  001ff542-0ee6-58e9-b674-208f6a5e294b  ...  POLYGON ((660072.922 1501476.795, 660051.306 1...
1  002ac31b-101e-537d-a3b4-0667280f5fc3  ...  POLYGON ((662162.934 1503314.935, 662117.905 1...
2  00458581-0c65-5601-8552-c457223d99fe  ...  POLYGON ((687414.87 1509653.403, 687417.385 15...
3  004709eb-d772-5da4-b008-333639afed53  ...  POLYGON ((688024.32 1521450.524, 688025.62 152...
4  00bf0971-e164-5c82-9181-6aed665b9e81  ...  POLYGON ((684185.567 1508479.553, 684215.404 1...

[5 rows x 4 columns]
Memory: 0.303249 MB

Written 3453 features to bma_landcover_intersect.geojson
Classes: {'barren': 820, 'forest': 777, 'shrub': 728, 'crop': 698, 'wetland': 212, 'urban': 160, 'mangrove': 29, 'grass': 29}
Total area: 980026998 m²
File size: 5.94 MB

cd /work && python3 -c "
import geopandas as gpd
gdf = gpd.read_file('bma_landcover_intersect.geojson')
print('Bounds:', gdf.total_bounds)
# Calculate approximate extent
x_range = gdf.total_bounds[2] - gdf.total_bounds[0]
y_range = gdf.total_bounds[3] - gdf.total_bounds[1]
print(f'Extent: {x_range:.0f} x {y_range:.0f} m')
# Check if all geometries are MultiPolygon
print('All MultiPolygon?', (gdf.geometry.type == 'MultiPolygon').all())
"

Written 3453 features to bma_landcover_intersect.geojson
Classes: {'barren': 820, 'forest': 777, 'shrub': 728, 'crop': 698, 'wetland': 212, 'urban': 160, 'mangrove': 29, 'grass': 29}
Total area: 980026998 m²
File size: 4.86 MB

Written 3453 features to bma_landcover_intersect.geojson
Classes: {'barren': 820, 'forest': 777, 'shrub': 728, 'crop': 698, 'wetland': 212, 'urban': 160, 'mangrove': 29, 'grass': 29}
Total area: 980026998 m²
File size: 2.46 MB

cd /work && python3 -c "
import geopandas as gpd
gdf = gpd.read_file('bma_landcover_intersect.geojson')
print('Features:', len(gdf))
print('Columns:', gdf.columns.tolist())
print('All MultiPolygon:', (gdf.geometry.type == 'MultiPolygon').all())
print('CRS:', gdf.crs)
print('area_m2 range:', gdf['area_m2'].min(), '-', gdf['area_m2'].max())
print('class values:', sorted(gdf['class'].unique()))
print()
# Check the GeoJSON structure first few features
import json
with open('bma_landcover_intersect.geojson') as f:
    gj = json.load(f)
print('Type:', gj['type'])
print('Num features in JSON:', len(gj['features']))
print()
feat = gj['features'][0]
print('First feature keys:', feat.keys())
print('Geometry type:', feat['geometry']['type'])
print('Properties:', feat['properties'])
"

Features: 3453
Columns: ['class', 'area_m2', 'geometry']
All MultiPolygon: True
CRS: EPSG:32647
area_m2 range: 0.15202269192047457 - 628183755.0187986
class values: ['barren', 'crop', 'forest', 'grass', 'mangrove', 'shrub', 'urban', 'wetland']

Type: FeatureCollection
Num features in JSON: 3453

First feature keys: dict_keys(['type', 'properties', 'geometry'])
Geometry type: MultiPolygon
Properties: {'class': 'forest', 'area_m2': 15735.747013974595}

Written 3453 features to bma_landcover_intersect.geojson
Classes: {'barren': 820, 'forest': 777, 'shrub': 728, 'crop': 698, 'wetland': 212, 'urban': 160, 'mangrove': 29, 'grass': 29}
Total area: 980026998 m²
File size: 2.46 MB

cd /work && python3 -c "
import geopandas as gpd
import json

gdf = gpd.read_file('bma_landcover_intersect.geojson')

# Final checks
print('✅ Features:', len(gdf))
print('✅ Columns:', gdf.columns.tolist())
print('✅ All MultiPolygon:', (gdf.geometry.type == 'MultiPolygon').all())
print('✅ CRS:', gdf.crs)
print('✅ Classes:', sorted(gdf['class'].unique()))
print('✅ area_m2 range:', f'{gdf[\"area_m2\"].min():.1f} – {gdf[\"area_m2\"].max():.1f} m²')
print('✅ Total area:', f'{gdf[\"area_m2\"].sum() / 1e6:.1f} km²')

# Validate GeoJSON structure
with open('bma_landcover_intersect.geojson') as f:
    gj = json.load(f)
assert gj['type'] == 'FeatureCollection'
for feat in gj['features']:
    assert feat['geometry']['type'] == 'MultiPolygon'
    assert set(feat['properties'].keys()) == {'class', 'area_m2'}
print('✅ GeoJSON structure valid — all features are MultiPolygon with class + area_m2')
"

✅ Features: 3453
✅ Columns: ['class', 'area_m2', 'geometry']
✅ All MultiPolygon: True
✅ CRS: EPSG:32647
✅ Classes: ['barren', 'crop', 'forest', 'grass', 'mangrove', 'shrub', 'urban', 'wetland']
✅ area_m2 range: 0.2 – 628183755.0 m²
✅ Total area: 980.0 km²
✅ GeoJSON structure valid — all features are MultiPolygon with class + area_m2