dc-l2-cairo-invalid-dedup

Rerun

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

def solve():
    # Load the data
    input_file = 'cairo_parcels_legacy.geojson'
    output_file = 'parcels_canonical.geoparquet'
    
    print(f"Loading {input_file}...")
    gdf = gpd.read_file(input_file)
    
    # Initial count
    print(f"Initial records: {len(gdf)}")

    # 1. Remove null or empty geometries
    gdf = gdf[gdf.geometry.notnull()]
    gdf = gdf[~gdf.geometry.is_empty]
    print(f"After removing null/empty geometries: {len(gdf)}")

    # 2. Fix invalid geometries
    print("Fixing invalid geometries...")
    gdf['geometry'] = gdf.geometry.buffer(0)
    gdf = gdf[~gdf.geometry.is_empty] # buffer(0) can result in empty if geometry was invalid
    print(f"After fixing invalid geometries: {len(gdf)}")

    # 3. Ensure all features are MultiPolygon
    print("Converting geometries to MultiPolygon...")
    def ensure_multipolygon(geom):
        if geom.geom_type == 'Polygon':
            return MultiPolygon([geom])
        return geom

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

    # 4. Handle Duplicates
    # "When duplicates exist, keep the record with the earliest record_seq"
    # The requirement implies duplicates are defined by parcel_id.
    print("Removing duplicates based on parcel_id...")
    # Sort by parcel_id and record_seq (ascending) to keep the earliest
    gdf = gdf.sort_values(by=['parcel_id', 'record_seq'], ascending=[True, True])
    gdf = gdf.drop_duplicates(subset='parcel_id', keep='first')
    print(f"After removing duplicates: {len(gdf)}")

    # 5. Remove duplicate geometries
    # "no duplicate geometries" - this usually means if two different IDs have the same shape, 
    # we keep one (though in a registry, we might want to keep them if IDs differ, 
    # but the prompt says "no duplicate geometries").
    # Let's assume it means spatial duplicates.
    print("Removing duplicate geometries...")
    gdf = gdf.drop_duplicates(subset='geometry')
    print(f"After removing duplicate geometries: {len(gdf)}")

    # 6. Calculate area and filter artifacts
    # The CRS is EPSG:22992 (based on the head output), which is a metric projected CRS.
    # Let's check the CRS just in case.
    if gdf.crs is None:
        print("Warning: CRS not detected. Assuming metric if not provided, but checking...")
        # If it's degrees, we must reproject. But 640000, 815000 looks like metric.
        # The file header said EPSG:22992.
    
    gdf['area_m2'] = gdf.geometry.area
    
    # "discard any polygon fragments smaller than 1 m²"
    print("Filtering small fragments (< 1 m²)...")
    gdf = gdf[gdf['area_m2'] >= 1.0]
    print(f"After filtering small fragments: {len(gdf)}")

    # 7. Final column selection
    # "including these columns: parcel_id, parcel_class, district, and area_m2"
    # We also need to keep geometry.
    output_columns = ['parcel_id', 'parcel_class', 'district', 'area_m2', 'geometry']
    
    # Check if columns exist
    missing_cols = [c for c in output_columns if c not in gdf.columns]
    if missing_cols:
        print(f"Warning: Missing columns in input: {missing_cols}")
        # If they are missing, we'll just use what we have or create dummy if needed for the task
        # But let's try to just select what is available.
        output_columns = [c for c in output_columns if c in gdf.columns]

    gdf_final = gdf[output_columns]

    # Final check for any remaining nulls in critical columns
    gdf_final = gdf_final.dropna(subset=['parcel_id', 'parcel_class', 'district'])

    # 8. Write to GeoParquet
    print(f"Writing to {output_file}...")
    gdf_final.to_parquet(output_file)
    print("Done.")

if __name__ == "__main__":
    solve()