fio-l2-capetown-landuse-dissolve

Rerun

import geopandas as gpd

def main():
    # Load the data
    gdf = gpd.read_file('capetown_landuse.fgb')

    # Ensure we are in a projected CRS to calculate area in m2 correctly
    # If the data is already projected (UTM/Lo), we'll use it.
    # If it's EPSG:4326, we must reproject.
    if gdf.crs.is_geographic:
        print("Reprojecting to UTM Zone 34S (EPSG:32734) for accurate area calculation...")
        gdf = gdf.to_crs(epsg=32734)
    else:
        print(f"Using existing CRS: {gdf.crs}")

    # Group by 'class'
    # We need: 
    # 1. Dissolved geometry (MultiPolygon)
    # 2. Total area (sum of areas)
    # 3. Count of original parcels

    # First, calculate area of each individual parcel in the projected CRS
    gdf['area_m2_parcel'] = gdf.geometry.area

    # Grouping and aggregating
    # dissolved: we'll use unary_union via groupby.apply or dissolve()
    # total area: sum of parcel areas
    # parcel count: count of rows per group

    print("Dissolving landuse classes...")
    # The dissolve() method is highly efficient for this
    dissolved = gdf.dissolve(by='class', aggfunc={
        'area_m2_parcel': 'sum',
        'id': 'count'
    })

    # Rename columns to match requirements
    dissolved = dissolved.rename(columns={
        'area_m2_parcel': 'area_m2',
        'id': 'parcel_count'
    })

    # Reset index so 'class' becomes a column
    dissolved = dissolved.reset_index()

    # Select/Order columns as requested
    # requested: 'class', 'area_m2', 'parcel_count'
    result = dissolved[['class', 'area_m2', 'parcel_count', 'geometry']]

    # Write to GeoParquet
    print(f"Writing result to landuse_dissolved.geoparquet...")
    result.to_parquet('landuse_dissolved.geoparquet')
    print("Done.")

if __name__ == "__main__":
    main()