fio-l2-capetown-landuse-dissolve

Rerun

import geopandas as gpd

def solve():
    # Load the landuse data
    input_file = 'capetown_landuse.fgb'
    output_file = 'landuse_dissolved.geoparquet'
    
    gdf = gpd.read_file(input_file)
    
    # Ensure we are working with a metric CRS for area calculations
    # The input is EPSG:32734, which is UTM zone 34S (metric)
    # We'll verify if it's already metric or if we need to project.
    # Since the print showed large numbers like 872000, it's likely metric.
    
    # Group by 'class'
    # We need:
    # 1. Unified geometry (dissolved)
    # 2. Total area in m2
    # 3. Count of source parcels (rows in the original gdf per class)
    
    # To get the count and area, we can aggregate first, then dissolve.
    # Or dissolve and then calculate area. 
    # But we need the count of original parcels.
    
    # Step 1: Aggregate to get count and area per class
    # Note: area must be calculated on individual geometries or the dissolved one.
    # Let's calculate area for each parcel first to be safe.
    gdf['area_m2_indiv'] = gdf.geometry.area
    
    summary = gdf.groupby('class').agg(
        area_m2=('area_m2_indiv', 'sum'),
        parcel_count=('id', 'count')
    ).reset_index()
    
    # Step 2: Dissolve geometries by class
    # We dissolve the original gdf to get the unified MultiPolygon
    dissolved_gdf = gdf.dissolve(by='class')
    
    # The dissolved_gdf index is 'class'. We need to merge the summary stats.
    # dissolved_gdf currently has only 'geometry' (and whatever else was in gdf)
    # Let's make it cleaner.
    
    # Re-dissolve to ensure we only have the geometry and the class index
    dissolved_gdf = gdf.dissolve(by='class')
    
    # Merge the summary stats back
    # summary has 'class', 'area_m2', 'parcel_count'
    # dissolved_gdf has index 'class'
    
    final_gdf = dissolved_gdf[['geometry']].merge(summary, left_index=True, right_on='class')
    
    # Reorder columns to match requested: class, area_m2, parcel_count, geometry
    # (GeoDataFrame will keep geometry as the active geometry column)
    final_gdf = final_gdf[['class', 'area_m2', 'parcel_count', 'geometry']]
    
    # Set geometry back correctly
    final_gdf = final_gdf.set_geometry('geometry')
    
    # Write to GeoParquet
    final_gdf.to_parquet(output_file)
    print(f"Successfully wrote {len(final_gdf)} rows to {output_file}")
    print(final_gdf.head())

if __name__ == "__main__":
    solve()