geo-l2-bangkok-landuse-intersect

Rerun

import geopandas as gpd

# Load the datasets
landcover = gpd.read_parquet('bangkok_landcover.parquet')
study_area = gpd.read_file('bma_study_area.geojson')

print("Landcover columns:", landcover.columns)
print("Study area columns:", study_area.columns)
print("Landcover CRS:", landcover.crs)
print("Study area CRS:", study_area.crs)
print("Landcover head:\n", landcover.head())

import geopandas as gpd
from shapely.geometry import MultiPolygon

def solve():
    # Load the datasets
    landcover = gpd.read_parquet('bangkok_landcover.parquet')
    study_area = gpd.read_file('bma_study_area.geojson')

    # Both are already in EPSG:32647 (UTM zone 47N) based on previous check
    # Ensure they match just in case
    if landcover.crs != study_area.crs:
        study_area = study_area.to_crs(landcover.crs)

    # 1. Intersection: Get landcover within the study area
    # We use overlay to handle the intersection of polygons
    intersected = gpd.overlay(landcover, study_area, how='intersection')

    # 2. Calculate area in m2
    # Since the CRS is UTM (meters), we can calculate area directly
    intersected['area_m2'] = intersected.geometry.area

    # 3. Simplify the geometry
    # A reasonable simplification tolerance for a "small file for policy lead" 
    # might be 1 to 5 meters in UTM. Let's try 1.0 meter first.
    # Also, ensure everything is a MultiPolygon.
    
    # Simplify
    simplified = intersected.simplify(tolerance=1.0, preserve_topology=True)

    # 4. Convert all geometries to MultiPolygon
    def ensure_multipolygon(geom):
        if geom.geom_type == 'Polygon':
            return MultiPolygon([geom])
        return geom

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

    # 5. Select and reorder columns
    # We need 'class' and 'area_m2'
    result = simplified[['class', 'area_m2', 'geometry']]

    # 6. Convert to WGS84 (EPSG:4326) for GeoJSON standard/web browsing
    result = result.to_crs(epsg=4326)

    # 7. Save to GeoJSON
    result.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')
    
    print(f"Successfully created bma_landcover_intersect.geojson")
    print(f"Number of features: {len(result)}")
    print(f"Columns: {result.columns.tolist()}")

if __name__ == "__main__":
    solve()

import geopandas as gpd
from shapely.geometry import MultiPolygon

def solve():
    # Load the datasets
    landcover = gpd.read_parquet('bangkok_landcover.parquet')
    study_area = gpd.read_file('bma_study_area.geojson')

    # Both are already in EPSG:32647 (UTM zone 47N)
    if landcover.crs != study_area.crs:
        study_area = study_area.to_crs(landcover.crs)

    # 1. Intersection: Get landcover within the study area
    # overlay returns a GeoDataFrame
    intersected = gpd.overlay(landcover, study_area, how='intersection')

    # 2. Calculate area in m2
    intersected['area_m2'] = intersected.geometry.area

    # 3. Simplify the geometry
    # Use a tolerance that makes the file smaller. 
    # Let's try a slightly higher tolerance if needed, but start with 1.0.
    simplified = intersected.simplify(tolerance=1.0, preserve_topology=True)

    # 4. Convert all geometries to MultiPolygon
    # It is safer to create a new GeoSeries or just update the geometry column
    def ensure_multipolygon(geom):
        if geom.geom_type == 'Polygon':
            return MultiPolygon([geom])
        return geom

    new_geoms = simplified.geometry.apply(ensure_multipolygon)
    
    # Re-construct the GeoDataFrame to avoid issues with setting the geometry column on a potentially "broken" object
    result = gpd.GeoDataFrame(simplified.drop(columns='geometry').copy(), geometry=new_geoms, crs=simplified.crs)

    # 5. Select and reorder columns
    # We need 'class' and 'area_m2'
    result = result[['class', 'area_m2', 'geometry']]

    # 6. Convert to WGS84 (EPSG:4326) for GeoJSON standard/web browsing
    result = result.to_crs(epsg=4326)

    # 7. Save to GeoJSON
    result.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')
    
    print(f"Successfully created bma_landcover_intersect.geojson")
    print(f"Number of features: {len(result)}")
    print(f"Columns: {result.columns.tolist()}")

if __name__ == "__main__":
    solve()

import geopandas as gpd
from shapely.geometry import MultiPolygon

def solve():
    # Load the datasets
    landcover = gpd.read_parquet('bangkok_landcover.parquet')
    study_area = gpd.read_file('bma_study_area.geojson')

    # Both are already in EPSG:32647 (UTM zone 47N)
    if landcover.crs != study_area.crs:
        study_area = study_area.to_crs(landcover.crs)

    # 1. Intersection: Get landcover within the study area
    intersected = gpd.overlay(landcover, study_area, how='intersection')

    # 2. Calculate area in m2
    intersected['area_m2'] = intersected.geometry.area

    # 3. Simplify the geometry
    # Using a tolerance that makes the file smaller.
    simplified = intersected.simplify(tolerance=1.0, preserve_topology=True)

    # 4. Convert all geometries to MultiPolygon
    def ensure_multipolygon(geom):
        if geom.is_empty:
            return geom
        if geom.geom_type == 'Polygon':
            return MultiPolygon([geom])
        return geom

    new_geoms = simplified.geometry.apply(ensure_multipolygon)
    
    # 5. Create the final GeoDataFrame with desired columns
    # We need 'class' from the original landcover and 'area_m2'
    # We must ensure 'class' exists in 'simplified'
    # In the previous attempt, it seems 'class' was lost or renamed. 
    # Let's check columns.
    
    # Note: gpd.overlay might name columns differently if they collide.
    # Let's check the columns of 'intersected' first.
    
    # For now, let's just re-create it properly.
    result_df = simplified.copy()
    result_df['geometry'] = new_geoms
    
    # 6. Select columns and ensure they exist
    final_cols = ['class', 'area_m2', 'geometry']
    # Verify columns
    for col in ['class', 'area_m2']:
        if col not in result_df.columns:
             raise KeyError(f"Column {col} missing. Available: {result_df.columns.tolist()}")

    result = result_df[final_cols]

    # 7. Convert to WGS84 (EPSG:4326) for GeoJSON standard/web browsing
    result = result.to_crs(epsg=4326)

    # 8. Save to GeoJSON
    result.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')
    
    print(f"Successfully created bma_landcover_intersect.geojson")
    print(f"Number of features: {len(result)}")
    print(f"Columns: {result.columns.tolist()}")

if __name__ == "__main__":
    solve()

import geopandas as gpd
from shapely.geometry import MultiPolygon

def solve():
    # Load the datasets
    landcover = gpd.read_parquet('bangkok_landcover.parquet')
    study_area = gpd.read_file('bma_study_area.geojson')

    # Both are already in EPSG:32647 (UTM zone 47N)
    if landcover.crs != study_area.crs:
        study_area = study_area.to_crs(landcover.crs)

    # 1. Intersection: Get landcover within the study area
    # overlay returns a GeoDataFrame
    intersected = gpd.overlay(landcover, study_area, how='intersection')

    # 2. Calculate area in m2
    intersected['area_m2'] = intersected.geometry.area

    # 3. Simplify the geometry
    # Using a tolerance that makes the file smaller.
    simplified = intersected.simplify(tolerance=1.0, preserve_topology=True)

    # 4. Convert all geometries to MultiPolygon
    def ensure_multipolygon(geom):
        if geom.is_empty:
            return geom
        if geom.geom_type == 'Polygon':
            return MultiPolygon([geom])
        return geom

    # Instead of trying to assign back to a GeoDataFrame which is causing issues
    # We will create a new GeoDataFrame from scratch using the attributes and the new geometries.
    
    new_geoms = simplified.geometry.apply(ensure_multipolygon)
    
    # Create a new GeoDataFrame to avoid any geometry assignment issues
    # We take the attributes we need from the simplified GeoDataFrame
    result = gpd.GeoDataFrame(
        {
            'class': simplified['class'],
            'area_m2': simplified['area_m2']
        },
        geometry=new_geoms,
        crs=simplified.crs
    )

    # 5. Convert to WGS84 (EPSG:4326) for GeoJSON standard/web browsing
    result = result.to_crs(epsg=4326)

    # 6. Save to GeoJSON
    result.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')
    
    print(f"Successfully created bma_landcover_intersect.geojson")
    print(f"Number of features: {len(result)}")
    print(f"Columns: {result.columns.tolist()}")

if __name__ == "__main__":
    solve()

import geopandas as gpd
from shapely.geometry import MultiPolygon

def solve():
    # Load the datasets
    landcover = gpd.read_parquet('bangkok_landcover.parquet')
    study_area = gpd.read_file('bma_study_area.geojson')

    print("Landcover columns:", landcover.columns.tolist())
    print("Study area columns:", study_area.columns.tolist())

    # Both are already in EPSG:32647 (UTM zone 47N)
    if landcover.crs != study_area.crs:
        study_area = study_area.to_crs(landcover.crs)

    # 1. Intersection: Get landcover within the study area
    # overlay returns a GeoDataFrame
    intersected = gpd.overlay(landcover, study_area, how='intersection')

    print("Intersected columns:", intersected.columns.tolist())

    # 2. Calculate area in m2
    intersected['area_m2'] = intersected.geometry.area

    # 3. Simplify the geometry
    # Using a tolerance that makes the file smaller.
    simplified = intersected.simplify(tolerance=1.0, preserve_topology=True)

    # 4. Convert all geometries to MultiPolygon
    def ensure_multipolygon(geom):
        if geom.is_empty:
            return geom
        if geom.geom_type == 'Polygon':
            return MultiPolygon([geom])
        return geom

    new_geoms = simplified.geometry.apply(ensure_multipolygon)
    
    # 5. Create the final GeoDataFrame with desired columns
    # We'll check which columns actually exist to avoid KeyError
    available_cols = simplified.columns.tolist()
    target_cols = ['class', 'area_m2']
    existing_target_cols = [c for c in target_cols if c in available_cols]
    
    # If 'class' is missing, let's see what is there.
    if 'class' not in existing_target_cols:
        # Maybe it's 'Class' or something else?
        print("Warning: 'class' not found in columns. Checking available columns again.")

    # Construct the result
    result_data = {col: simplified[col].values for col in existing_target_cols}
    result = gpd.GeoDataFrame(result_data, geometry=new_geoms, crs=simplified.crs)

    # 6. Convert to WGS84 (EPSG:4326) for GeoJSON standard/web browsing
    result = result.to_crs(epsg=4326)

    # 7. Save to GeoJSON
    result.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')
    
    print(f"Successfully created bma_landcover_intersect.geojson")
    print(f"Number of features: {len(result)}")
    print(f"Columns: {result.columns.tolist()}")

if __name__ == "__main__":
    solve()

import geopandas as gpd
from shapely.geometry import MultiPolygon

def solve():
    # Load the datasets
    landcover = gpd.read_parquet('bangkok_landcover.parquet')
    study_area = gpd.read_file('bma_study_area.geojson')

    # Both are already in EPSG:32647 (UTM zone 47N)
    if landcover.crs != study_area.crs:
        study_area = study_area.to_crs(landcover.crs)

    # 1. Intersection: Get landcover within the study area
    intersected = gpd.overlay(landcover, study_area, how='intersection')

    # 2. Calculate area in m2
    intersected['area_m2'] = intersected.geometry.area

    # 3. Simplify the geometry
    simplified = intersected.simplify(tolerance=1.0, preserve_topology=True)

    # 4. Convert all geometries to MultiPolygon
    def ensure_multipolygon(geom):
        if geom.is_empty:
            return geom
        if geom.geom_type == 'Polygon':
            return MultiPolygon([geom])
        return geom

    new_geoms = simplified.geometry.apply(ensure_multipolygon)
    
    # 5. Create the final GeoDataFrame
    # We want 'class' and 'area_m2'
    # Note: 'simplified' is a GeoDataFrame.
    
    result = gpd.GeoDataFrame(
        {'class': simplified['class'], 'area_m2': simplified['area_m2']},
        geometry=new_geoms,
        crs=simplified.crs
    )

    # 6. Convert to WGS84 (EPSG:4326) for GeoJSON standard/web browsing
    result = result.to_crs(epsg=4326)

    # 7. Save to GeoJSON
    result.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')
    
    print(f"Successfully created bma_landcover_intersect.geojson")
    print(f"Number of features: {len(result)}")
    print(f"Columns: {result.columns.tolist()}")

if __name__ == "__main__":
    solve()

import geopandas as gpd
from shapely.geometry import MultiPolygon

def solve():
    # Load the datasets
    landcover = gpd.read_parquet('bangkok_landcover.parquet')
    study_area = gpd.read_file('bma_study_area.geojson')

    # Both are already in EPSG:32647 (UTM zone 47N)
    if landcover.crs != study_area.crs:
        study_area = study_area.to_crs(landcover.crs)

    # 1. Intersection: Get landcover within the study area
    intersected = gpd.overlay(landcover, study_area, how='intersection')

    # 2. Calculate area in m2
    intersected['area_m2'] = intersected.geometry.area

    # 3. Simplify the geometry
    simplified = intersected.simplify(tolerance=1.0, preserve_topology=True)

    # 4. Convert all geometries to MultiPolygon
    def ensure_multipolygon(geom):
        if geom.is_empty:
            return geom
        if geom.geom_type == 'Polygon':
            return MultiPolygon([geom])
        return geom

    new_geoms = simplified.geometry.apply(ensure_multipolygon)
    
    # 5. Create the final GeoDataFrame
    # We'll use a more robust way to copy the columns.
    # Let's just create a new dict with the values.
    
    # Check if 'class' exists in simplified
    if 'class' in simplified.columns:
        classes = simplified['class']
    else:
        # This is strange, but let's handle it
        print("Error: 'class' not found in simplified.")
        print("Columns found:", simplified.columns.tolist())
        return

    result = gpd.GeoDataFrame(
        {'class': classes.values, 'area_m2': simplified['area_m2'].values},
        geometry=new_geoms.values,
        crs=simplified.crs
    )

    # 6. Convert to WGS84 (EPSG:4326) for GeoJSON standard/web browsing
    result = result.to_crs(epsg=4326)

    # 7. Save to GeoJSON
    result.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')
    
    print(f"Successfully created bma_landcover_intersect.geojson")
    print(f"Number of features: {len(result)}")
    print(f"Columns: {result.columns.tolist()}")

if __name__ == "__main__":
    solve()

import geopandas as gpd
from shapely.geometry import MultiPolygon

def solve():
    # Load the datasets
    landcover = gpd.read_parquet('bangkok_landcover.parquet')
    study_area = gpd.read_file('bma_study_area.geojson')

    # Both are already in EPSG:32647 (UTM zone 47N)
    if landcover.crs != study_area.crs:
        study_area = study_area.to_crs(landcover.crs)

    # 1. Intersection: Get landcover within the study area
    intersected = gpd.overlay(landcover, study_area, how='intersection')

    # 2. Calculate area in m2
    intersected['area_m2'] = intersected.geometry.area

    # 3. Simplify the geometry
    simplified = intersected.simplify(tolerance=1.0, preserve_topology=True)

    # 4. Convert all geometries to MultiPolygon
    def ensure_multipolygon(geom):
        if geom.is_empty:
            return geom
        if geom.geom_type == 'Polygon':
            return MultiPolygon([geom])
        return geom

    new_geoms = simplified.geometry.apply(ensure_multipolygon)
    
    # 5. Create the final GeoDataFrame
    # We'll use a more robust way to copy the columns.
    
    # Create a new GeoDataFrame from the simplified one
    # but only with the columns we want, and with the updated geometry
    
    # The issue was that 'simplified' was behaving like a GeoSeries in some contexts? 
    # Actually, simplified = intersected.simplify(...) returns a GeoSeries in some versions of geopandas if not careful.
    # But overlay returns a GeoDataFrame. Let's check type.
    
    if not isinstance(simplified, gpd.GeoDataFrame):
        print(f"Warning: simplified is {type(simplified)}, not GeoDataFrame")
        # If it's a GeoSeries, it's probably because simplify() was called on a GeoSeries.
        # But intersected is a GeoDataFrame.
        # Wait, if intersected is a GeoDataFrame, simplified = intersected.simplify() 
        # should return a GeoDataFrame. 
        # Let's check.
        pass

    # Let's try a different approach for step 5.
    result = simplified.copy()
    result['geometry'] = new_geoms
    
    # Filter to required columns
    # Using a list to ensure we only take what we need
    result = result[['class', 'area_m2', 'geometry']]

    # 6. Convert to WGS84 (EPSG:4326) for GeoJSON standard/web browsing
    result = result.to_crs(epsg=4326)

    # 7. Save to GeoJSON
    result.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')
    
    print(f"Successfully created bma_landcover_intersect.geojson")
    print(f"Number of features: {len(result)}")
    print(f"Columns: {result.columns.tolist()}")

if __name__ == "__main__":
    solve()

import geopandas as gpd
from shapely.geometry import MultiPolygon

def solve():
    # Load the datasets
    landcover = gpd.read_parquet('bangkok_landcover.parquet')
    study_area = gpd.read_file('bma_study_area.geojson')

    # Both are already in EPSG:32647 (UTM zone 47N)
    if landcover.crs != study_area.crs:
        study_area = study_area.to_crs(landcover.crs)

    # 1. Intersection: Get landcover within the study area
    # overlay returns a GeoDataFrame
    intersected = gpd.overlay(landcover, study_area, how='intersection')

    # 2. Calculate area in m2
    intersected['area_m2'] = intersected.geometry.area

    # 3. Simplify the geometry
    # Note: simplify() on a GeoDataFrame returns a GeoDataFrame.
    # However, the error message said simplified is a GeoSeries.
    # This happens if we do something like intersected.geometry.simplify().
    # Let's be explicit.
    
    # We want to simplify the geometries but keep the attributes.
    # The correct way to simplify a GeoDataFrame is intersected.simplify(tolerance=...)
    # BUT if it returns a GeoSeries, we need to handle it.
    
    # Let's try to avoid the ambiguity by simplifying the geometry column and then 
    # rebuilding the GeoDataFrame.
    
    new_geoms = intersected.geometry.simplify(tolerance=1.0, preserve_topology=True)
    
    # 4. Convert all geometries to MultiPolygon
    def ensure_multipolygon(geom):
        if geom.is_empty:
            return geom
        if geom.geom_type == 'Polygon':
            return MultiPolygon([geom])
        return geom

    new_geoms = new_geoms.apply(ensure_multipolygon)
    
    # 5. Create the final GeoDataFrame
    # Use the original 'intersected' for attributes and 'new_geoms' for geometry
    result = gpd.GeoDataFrame(
        intersected[['class', 'area_m2']].copy(),
        geometry=new_geoms,
        crs=intersected.crs
    )

    # 6. Convert to WGS84 (EPSG:4326) for GeoJSON standard/web browsing
    result = result.to_crs(epsg=4326)

    # 7. Save to GeoJSON
    result.to_file('bma_landcover_intersect.geojson', driver='GeoJSON')
    
    print(f"Successfully created bma_landcover_intersect.geojson")
    print(f"Number of features: {len(result)}")
    print(f"Columns: {result.columns.tolist()}")

if __name__ == "__main__":
    solve()