geo-l2-bangkok-landuse-intersect

Rerun

import geopandas as gpd
import pandas as pd

# Check landcover columns
lc = gpd.read_parquet('bangkok_landcover.parquet')
print("Landcover columns:", lc.columns.tolist())
print("Landcover head:\n", lc.head())

# Check study area columns
sa = gpd.read_file('bma_study_area.geojson')
print("\nStudy area columns:", sa.columns.tolist())
print("Study area head:\n", sa.head())

import geopandas as gpd
from shapely.geometry import MultiPolygon

def solve():
    # 1. Load data
    print("Loading data...")
    lc = gpd.read_parquet('bangkok_landcover.parquet')
    sa = gpd.read_file('bma_study_area.geojson')

    # Ensure same CRS
    if lc.crs != sa.crs:
        print(f"Reprojecting study area from {sa.crs} to {lc.crs}")
        sa = sa.to_crs(lc.crs)

    # 2. Intersect landcover with study area
    print("Intersecting landcover with study area...")
    # Use overlay for intersection to get clipped geometries
    # Note: overlay returns a GeoDataFrame of the intersection
    intersected = gpd.overlay(lc, sa, how='intersection')

    if intersected.empty:
        print("No intersection found.")
        return

    # 3. Calculate area in m2
    # Since it's already in a metric CRS (based on coordinate values in preview)
    # but let's double check. The preview shows 647454 1506832 which looks like UTM/metric.
    print("Calculating area...")
    intersected['area_m2'] = intersected.geometry.area

    # 4. Simplify geometries to reduce file size
    # We need to pick a tolerance. Since coordinates are in meters, 
    # a tolerance of 1-5 meters should be invisible for a browser preview
    # but significantly reduces vertex count.
    print("Simplifying geometries...")
    tolerance = 1.0 
    intersected['geometry'] = intersected.geometry.simplify(tolerance, preserve_topology=True)

    # 5. Convert everything to MultiPolygon as requested
    print("Converting to MultiPolygons...")
    def ensure_multipolygon(geom):
        if geom.geom_type == 'Polygon':
            return MultiPolygon([geom])
        return geom

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

    # 6. Select only required columns
    # "with the original 'class' string and a per-feature 'area_m2' in square metres"
    # The overlay might have kept 'name' from sa, so we filter.
    final_gdf = intersected[['class', 'area_m2', 'geometry']]

    # 7. Write to GeoJSON
    print(f"Writing to bma_landcover_intersect.geojson...")
    final_gdf.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')
    print("Done!")

if __name__ == "__main__":
    solve()

import geopandas as gpd
from shapely.geometry import MultiPolygon

def solve():
    # 1. Load data
    print("Loading data...")
    lc = gpd.read_parquet('bangkok_landcover.parquet')
    sa = gpd.read_file('bma_study_area.geojson')

    # Ensure same CRS
    if lc.crs != sa.crs:
        print(f"Reprojecting study area from {sa.crs} to {lc.crs}")
        sa = sa.to_crs(lc.crs)

    # 2. Intersect landcover with study area
    print("Intersecting landcover with study area...")
    intersected = gpd.overlay(lc, sa, how='intersection')

    if intersected.empty:
        print("No intersection found.")
        return

    # 3. Calculate area in m2
    print("Calculating area...")
    intersected['area_m2'] = intersected.geometry.area

    # 4. Simplify geometries
    # Increase tolerance to 5.0m to ensure small file size for browser preview
    print("Simplifying geometries (tolerance=5.0)...")
    tolerance = 5.0 
    intersected['geometry'] = intersected.geometry.simplify(tolerance, preserve_topology=True)

    # 5. Convert everything to MultiPolygon
    print("Converting to MultiPolygons...")
    def ensure_multipolygon(geom):
        if geom.is_empty:
            return None
        if geom.geom_type == 'Polygon':
            return MultiPolygon([geom])
        return geom

    intersected['geometry'] = intersected.geometry.apply(ensure_multipolygon)
    # Remove any None entries from invalid/empty simplifications
    intersected = intersected[intersected.geometry.notnull()]

    # 6. Select only required columns
    final_gdf = intersected[['class', 'area_m2', 'geometry']]

    # 7. Write to GeoJSON
    print(f"Writing to bma_landcover_intersect.geojson...")
    final_gdf.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')
    print("Done!")

if __name__ == "__main__":
    solve()

import geopandas as gpd
from shapely.geometry import MultiPolygon

def solve():
    # 1. Load data
    print("Loading data...")
    lc = gpd.read_parquet('bangkok_landcover.parquet')
    sa = gpd.read_file('bma_study_area.geojson')

    # Ensure same CRS
    if lc.crs != sa.crs:
        print(f"Reprojecting study area from {sa.crs} to {lc.crs}")
        sa = sa.to_crs(lc.crs)

    # 2. Intersect landcover with study area
    print("Intersecting landcover with study area...")
    # Using clip is faster and often more robust for simple intersection
    intersected = lc.clip(sa)

    if intersected.empty:
        print("No intersection found.")
        return

    # 3. Calculate area in m2
    print("Calculating area...")
    intersected['area_m2'] = intersected.geometry.area

    # 4. Simplify geometries
    # Try an even more aggressive simplification for the 'preview' requirement
    print("Simplifying geometries (tolerance=20.0)...")
    tolerance = 20.0 
    intersected['geometry'] = intersected.geometry.simplify(tolerance, preserve_topology=True)

    # 5. Convert everything to MultiPolygon
    print("Converting to MultiPolygons...")
    def ensure_multipolygon(geom):
        if geom.is_empty:
            return None
        if geom.geom_type == 'Polygon':
            return MultiPolygon([geom])
        return geom

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

    # 6. Select only required columns
    final_gdf = intersected[['class', 'area_m2', 'geometry']]

    # 7. Write to GeoJSON
    print(f"Writing to bma_landcover_intersect.geojson...")
    final_gdf.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')
    print("Done!")

if __name__ == "__main__":
    solve()

import geopandas as gpd
from shapely.geometry import MultiPolygon

def solve():
    # 1. Load data
    print("Loading data...")
    lc = gpd.read_parquet('bangkok_landcover.parquet')
    sa = gpd.read_file('bma_study_area.geojson')

    # Ensure same CRS
    if lc.crs != sa.crs:
        print(f"Reprojecting study area from {sa.crs} to {lc.crs}")
        sa = sa.to_crs(lc.crs)

    # Fix invalid geometries in both
    print("Fixing invalid geometries...")
    lc['geometry'] = lc.geometry.buffer(0)
    sa['geometry'] = sa.geometry.buffer(0)

    # 2. Intersect landcover with study area
    print("Intersecting landcover with study area...")
    # Using overlay instead of clip to avoid some GEOS issues, and we know overlay handles it
    # But we'll use the safe way with overlay
    intersected = gpd.overlay(lc, sa, how='intersection')

    if intersected.empty:
        print("No intersection found.")
        return

    # 3. Calculate area in m2
    print("Calculating area...")
    intersected['area_m2'] = intersected.geometry.area

    # 4. Simplify geometries
    # Use a tolerance that balances size and detail. 
    # 10m is a good middle ground for a 'preview'
    print("Simplifying geometries (tolerance=10.0)...")
    tolerance = 10.0 
    intersected['geometry'] = intersected.geometry.simplify(tolerance, preserve_topology=True)

    # 5. Convert everything to MultiPolygon
    print("Converting to MultiPolygons...")
    def ensure_multipolygon(geom):
        if geom.is_empty:
            return None
        if geom.geom_type == 'Polygon':
            return MultiPolygon([geom])
        return geom

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

    # 6. Select only required columns
    final_gdf = intersected[['class', 'area_m2', 'geometry']]

    # 7. Write to GeoJSON
    print(f"Writing to bma_landcover_intersect.geojson...")
    final_gdf.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')
    print("Done!")

if __name__ == "__main__":
    solve()

import geopandas as gpd
from shapely.geometry import MultiPolygon

def solve():
    # 1. Load data
    print("Loading data...")
    lc = gpd.read_parquet('bangkok_landcover.parquet')
    sa = gpd.read_file('bma_study_area.geojson')

    # Ensure same CRS
    if lc.crs != sa.crs:
        print(f"Reprojecting study area from {sa.crs} to {lc.crs}")
        sa = sa.to_crs(lc.crs)

    # Fix invalid geometries in both
    print("Fixing invalid geometries...")
    lc['geometry'] = lc.geometry.buffer(0)
    sa['geometry'] = sa.geometry.buffer(0)

    # 2. Intersect landcover with study area
    print("Intersecting landcover with study area...")
    intersected = gpd.overlay(lc, sa, how='intersection')

    if intersected.empty:
        print("No intersection found.")
        return

    # 3. Calculate area in m2
    print("Calculating area...")
    intersected['area_m2'] = intersected.geometry.area

    # 4. Simplify geometries
    # VERY aggressive simplification for a small file size as requested
    print("Simplifying geometries (tolerance=50.0)...")
    tolerance = 50.0 
    intersected['geometry'] = intersected.geometry.simplify(tolerance, preserve_topology=True)

    # 5. Convert everything to MultiPolygon
    print("Converting to MultiPolygons...")
    def ensure_multipolygon(geom):
        if geom.is_empty:
            return None
        if geom.geom_type == 'Polygon':
            return MultiPolygon([geom])
        return geom

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

    # 6. Select only required columns
    final_gdf = intersected[['class', 'area_m2', 'geometry']]

    # 7. Write to GeoJSON
    print(f"Writing to bma_landcover_intersect.geojson...")
    final_gdf.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')
    print("Done!")

if __name__ == "__main__":
    solve()

import geopandas as gpd
from shapely.geometry import MultiPolygon
import json

def solve():
    # 1. Load data
    print("Loading data...")
    lc = gpd.read_parquet('bangkok_landcover.parquet')
    sa = gpd.read_file('bma_study_area.geojson')

    # Ensure same CRS
    if lc.crs != sa.crs:
        print(f"Reprojecting study area from {sa.crs} to {lc.crs}")
        sa = sa.to_crs(lc.crs)

    # Fix invalid geometries in both
    print("Fixing invalid geometries...")
    lc['geometry'] = lc.geometry.buffer(0)
    sa['geometry'] = sa.geometry.buffer(0)

    # 2. Intersect landcover with study area
    print("Intersecting landcover with study area...")
    intersected = gpd.overlay(lc, sa, how='intersection')

    if intersected.empty:
        print("No intersection found.")
        return

    # 3. Calculate area in m2
    print("Calculating area...")
    intersected['area_m2'] = intersected.geometry.area

    # 4. Simplify geometries
    # Aggressive simplification for small file size
    print("Simplifying geometries (tolerance=50.0)...")
    tolerance = 50.0 
    intersected['geometry'] = intersected.geometry.simplify(tolerance, preserve_topology=True)

    # 5. Convert everything to MultiPolygon
    print("Converting to MultiPolygons...")
    def ensure_multipolygon(geom):
        if geom.is_empty:
            return None
        if geom.geom_type == 'Polygon':
            return MultiPolygon([geom])
        return geom

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

    # 6. Select only required columns
    final_gdf = intersected[['class', 'area_m2', 'geometry']]

    # 7. Write to GeoJSON
    print(f"Writing to bma_landcover_intersect.geojson...")
    # To ensure the output is strictly a GeoJSON FeatureCollection with MultiPolygons
    # we use to_file.
    final_gdf.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')
    
    # Final check on the output format
    with open('bma_landcover_intersect.geojson', 'r') as f:
        data = json.load(f)
        print(f"Output type: {data['type']}")
        if len(data['features']) > 0:
            print(f"First feature geometry type: {data['features'][0]['geometry']['type']}")
            print(f"Columns in first feature: {data['features'][0]['properties'].keys()}")

    print("Done!")

if __name__ == "__main__":
    solve()