geo-l2-bangkok-landuse-intersect

Rerun

import geopandas as gpd

# Read the input files
landcover = gpd.read_parquet('bangkok_landcover.parquet')
study_area = gpd.read_file('bma_study_area.geojson')

print("Landcover columns:", landcover.columns)
print("Landcover CRS:", landcover.crs)
print("Study Area CRS:", study_area.crs)

# Check if CRS match and transform if necessary
if landcover.crs != study_area.crs:
    print(f"Transforming landcover from {landcover.crs} to {study_area.crs}")
    landcover = landcover.to_crs(study_area.crs)

# Ensure study_area is a single polygon/multipolygon for intersection
# or at least handle it properly.
# We want landcover within study_area.
# Using overlay to perform the intersection.
intersected = gpd.overlay(landcover, study_area, how='intersection')

# Calculate area in square meters. 
# We should ensure we are in a projected CRS for accurate area calculation.
# If the input is in WGS84 (degrees), we must project to a metric CRS.
if intersected.crs.is_geographic:
    print("CRS is geographic. Projecting to EPSG:32647 (UTM zone 47N) for area calculation.")
    projected_crs = "EPSG:32647"
    intersected_projected = intersected.to_crs(projected_crs)
    intersected['area_m2'] = intersected_projected.geometry.area
else:
    intersected['area_m2'] = intersected.geometry.area

# The prompt asks for the original 'class' string and 'area_m2'.
# Let's verify 'class' exists.
if 'class' not in intersected.columns:
    # Check what the column name actually is
    print("Column 'class' not found. Available columns:", intersected.columns)
    # Maybe it's 'landcover_class' or something? 
    # Let's look at the first few rows of landcover to be sure.

# Simplify to keep file size small.
# Since we want to keep features somewhat meaningful, let's use a reasonable tolerance.
# A tolerance depends on the CRS units. If metric, maybe 1-5 meters?
# If we don't know the CRS units yet, let's check.
print("CRS info:", intersected.crs)

# Simplify geometry. 
# Note: simplify() can change geometry type. 
# We need to ensure everything is a MultiPolygon.
intersected['geometry'] = intersected.geometry.simplify(tolerance=1.0, preserve_topology=True)

# Convert all geometries to MultiPolygon
def ensure_multipolygon(geometry):
    if geometry.geom_type == 'Polygon':
        from shapely.geometry import MultiPolygon
        return MultiPolygon([geometry])
    return geometry

intersected['geometry'] = intersected.geometry.apply(ensure_multipolygon)

# Final selection of columns
# The prompt says: "original `class` string and a per-feature `area_m2`"
# We should also keep 'geometry'
output_gdf = intersected[['class', 'area_m2', 'geometry']]

# Write to GeoJSON
output_gdf.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')
print("Successfully wrote bma_landcover_intersect.geojson")

import geopandas as gpd
from shapely.geometry import MultiPolygon, Polygon

# Read the input files
landcover = gpd.read_parquet('bangkok_landcover.parquet')
study_area = gpd.read_file('bma_study_area.geojson')

# Check if CRS match and transform if necessary
if landcover.crs != study_area.crs:
    landcover = landcover.to_crs(study_area.crs)

# Perform intersection
intersected = gpd.overlay(landcover, study_area, how='intersection')

# Calculate area in square meters (input is already in EPSG:32647 which is metric)
intersected['area_m2'] = intersected.geometry.area

# Simplify geometry to reduce size
# Using a tolerance of 5 meters for more significant reduction while maintaining shape
intersected['geometry'] = intersected.geometry.simplify(tolerance=5.0, preserve_topology=True)

# Ensure every feature is a MultiPolygon as requested
def ensure_multipolygon(geometry):
    if geometry.is_empty:
        return geometry
    if geometry.geom_type == 'Polygon':
        return MultiPolygon([geometry])
    return geometry

intersected['geometry'] = intersected.geometry.apply(ensure_multipolygon)

# Final selection of columns
output_gdf = intersected[['class', 'area_m2', 'geometry']]

# Write to GeoJSON
# Note: GeoJSON standard uses WGS84 (degrees). 
# For the policy lead to preview in a browser, we should convert back to EPSG:4326.
output_gdf_wgs84 = output_gdf.to_crs("EPSG:4326")
output_gdf_wgs84.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')

print(f"Successfully wrote bma_landcover_intersect.geojson. Size: {output_gdf_wgs84.geometry.count} features.")

import geopandas as gpd
from shapely.geometry import MultiPolygon

# Read the input files
landcover = gpd.read_parquet('bangkok_landcover.parquet')
study_area = gpd.read_file('bma_study_area.geojson')

# Ensure CRS match
if landcover.crs != study_area.crs:
    landcover = landcover.to_crs(study_area.crs)

# Perform intersection
intersected = gpd.overlay(landcover, study_area, how='intersection')

# Calculate area in square meters using a metric projection (UTM 47N)
# The original data seems to be in EPSG:32647 based on previous logs
if intersected.crs.is_geographic:
    projected_crs = "EPSG:32647"
    intersected_projected = intersected.to_crs(projected_crs)
    intersected['area_m2'] = intersected_projected.geometry.area
else:
    intersected['area_m2'] = intersected.geometry.area

# Simplify geometry to reduce file size for web preview
# Tolerance is in CRS units (meters)
intersected['geometry'] = intersected.geometry.simplify(tolerance=10.0, preserve_topology=True)

# Ensure every feature is a MultiPolygon as requested
def ensure_multipolygon(geometry):
    if geometry.is_empty:
        return geometry
    if geometry.geom_type == 'Polygon':
        return MultiPolygon([geometry])
    return geometry

intersected['geometry'] = intersected.geometry.apply(ensure_multipolygon)

# Final selection of columns
output_gdf = intersected[['class', 'area_m2', 'geometry']]

# Convert back to WGS84 for GeoJSON standard (required for most web viewers)
output_gdf_wgs84 = output_gdf.to_crs("EPSG:4326")

# Write to GeoJSON
output_gdf_wgs84.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')

print(f"Successfully wrote bma_landcover_intersect.geojson")
print(f"Number of features: {len(output_gdf_wgs84)}")