dc-l2-cairo-invalid-dedup

Rerun

import geopandas as gpd

def inspect_data(file_path):
    gdf = gpd.read_file(file_path)
    print("--- Info ---")
    print(gdf.info())
    print("\n--- Head ---")
    print(gdf.head())
    print("\n--- Columns ---")
    print(gdf.columns)
    print("\n--- Geometry Types ---")
    print(gdf.geom_type.value_counts())
    print("\n--- Null Values ---")
    print(gdf.isnull().sum())

if __name__ == "__main__":
    inspect_data('cairo_parcels_legacy.geojson')

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

def clean_parcels(input_file, output_file):
    # 1. Load the data
    gdf = gpd.read_file(input_file)
    
    # 2. Remove null or empty geometries
    gdf = gdf[gdf.geometry.notnull() & ~gdf.geometry.is_empty]
    
    # 3. Fix invalid geometries
    gdf['geometry'] = gdf.geometry.make_valid()
    
    # 4. Discard any polygon fragments smaller than 1 m2
    # Note: A MultiPolygon might consist of several polygons. 
    # The requirement says "discard any polygon fragments smaller than 1 m2".
    # This usually implies we should check the components of the MultiPolygons.
    # However, if a MultiPolygon is considered a single feature, we should clarify.
    # Given the instruction "discard any polygon fragments", let's handle individual polygons.
    
    def filter_small_fragments(geom):
        if geom.is_empty:
            return None
        if isinstance(geom, Polygon):
            return geom if geom.area >= 1.0 else None
        elif isinstance(geom, MultiPolygon):
            # Keep only the parts of the MultiPolygon that are >= 1 m2
            parts = [p for p in geom.geoms if p.area >= 1.0]
            if not parts:
                return None
            return MultiPolygon(parts)
        return geom

    gdf['geometry'] = gdf['geometry'].apply(filter_small_fragments)
    gdf = gdf[gdf.geometry.notnull() & ~gdf.geometry.is_empty]

    # 5. Handle duplicates: one canonical record per parcel.
    # Keep the record with the earliest record_seq.
    gdf = gdf.sort_values('record_seq')
    gdf = gdf.drop_duplicates(subset=['parcel_id'], keep='first')
    
    # 6. Ensure every feature is a MultiPolygon
    def ensure_multipolygon(geom):
        if isinstance(geom, Polygon):
            return MultiPolygon([geom])
        return geom

    gdf['geometry'] = gdf['geometry'].apply(ensure_multipolygon)
    
    # 7. Remove duplicate geometries
    # We use a subset of geometry to identify duplicates
    gdf = gdf.drop_duplicates(subset=['geometry'])

    # 8. Final cleaning: remove any remaining empty/invalid (safety check)
    gdf = gdf[gdf.geometry.notnull() & ~gdf.geometry.is_empty]
    gdf['geometry'] = gdf.geometry.make_valid()
    
    # Ensure we only have MultiPolygons after make_valid (if make_valid returned Polygons)
    gdf['geometry'] = gdf['geometry'].apply(ensure_multipolygon)

    # 9. Select final columns
    # parcel_id, parcel_class, district, and area_m2
    # We should recalculate area_m2 to ensure it matches the cleaned geometry?
    # The instruction says "including these columns". It doesn't explicitly say 
    # to recalculate, but usually, in a registry, area should be accurate.
    # However, the original 'area_m2' might be useful. Let's see.
    # Let's recalculate to be safe and consistent with the cleaned geometry.
    gdf['area_m2'] = gdf.geometry.area

    final_gdf = gdf[['parcel_id', 'parcel_class', 'district', 'area_m2', 'geometry']]
    
    # 10. Save to GeoParquet
    final_gdf.to_parquet(output_file)
    print(f"Successfully saved {len(final_gdf)} records to {output_file}")

if __name__ == "__main__":
    clean_parcels('cairo_parcels_legacy.geojson', 'parcels_canonical.geoparquet')

import geopandas as gpd

def verify_output(file_path):
    gdf = gpd.read_parquet(file_path)
    
    print("--- Info ---")
    print(gdf.info())
    
    print("\n--- Head ---")
    print(gdf.head())
    
    print("\n--- Column Check ---")
    expected_cols = {'parcel_id', 'parcel_class', 'district', 'area_m2', 'geometry'}
    actual_cols = set(gdf.columns)
    print(f"Expected: {expected_cols}")
    print(f"Actual: {actual_cols}")
    assert expected_cols.issubset(actual_cols), "Missing columns!"
    
    print("\n--- Geometry Type Check ---")
    geom_types = gdf.geom_type.unique()
    print(f"Geometry types: {geom_types}")
    assert all(t == 'MultiPolygon' for t in geom_types), f"Unexpected geometry types: {geom_types}"
    
    print("\n--- Area Check ---")
    min_area = gdf.geometry.area.min()
    print(f"Minimum area: {min_area}")
    assert min_area >= 1.0, f"Found area smaller than 1m2: {min_area}"

    print("\n--- Null/Empty Check ---")
    print(f"Null geometries: {gdf.geometry.isnull().sum()}")
    print(f"Empty geometries: {gdf.geometry.is_empty.sum()}")
    assert gdf.geometry.isnull().sum() == 0, "Null geometries found!"
    assert gdf.geometry.is_empty.sum() == 0, "Empty geometries found!"

    print("\n--- Duplicate Geometry Check ---")
    # Note: This might fail if different parcels have the same geometry (unlikely but possible)
    # but the requirement says "no duplicate geometries".
    dupes = gdf.geometry.duplicated().sum()
    print(f"Duplicate geometries: {dupes}")
    assert dupes == 0, "Duplicate geometries found!"

    print("\n--- Duplicate Parcel ID Check ---")
    dupe_ids = gdf['parcel_id'].duplicated().sum()
    print(f"Duplicate parcel_ids: {dupe_ids}")
    assert dupe_ids == 0, "Duplicate parcel_ids found!"

    print("\nAll checks passed!")

if __name__ == "__main__":
    verify_output('parcels_canonical.geoparquet')