dc-l2-cairo-invalid-dedup

Rerun

import geopandas as gpd

# Load the data
gdf = gpd.read_file('cairo_parcels_legacy.geojson')

print("Initial info:")
print(gdf.info())
print("\nFirst few rows:")
print(gdf.head())
print("\nGeometry type counts:")
print(gdf.geom_type.value_counts())
print("\nColumns:")
print(gdf.columns)

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. Handle invalid geometries
    # Ensure geometries are valid. If not, try to fix them using buffer(0)
    gdf['geometry'] = gdf['geometry'].buffer(0)
    
    # 3. Remove null or empty geometries
    gdf = gdf[gdf.geometry.notnull() & ~gdf.geometry.is_empty]
    
    # 4. Handle duplicates based on parcel_id
    # Keep the record with the earliest record_seq
    gdf = gdf.sort_values('record_seq').drop_duplicates(subset='parcel_id', keep='first')
    
    # 5. Discard polygon fragments smaller than 1 m^2
    # Note: The requirement says "discard any polygon fragments smaller than 1 m^2".
    # Usually this implies that if a parcel is represented by multiple parts, 
    # we should filter the parts. However, the instructions also say "Every feature must be a MultiPolygon".
    # This suggests we should process the geometry components.
    
    def filter_small_parts(geom):
        if geom.is_empty:
            return geom
        
        # If it's a Polygon, check its area
        if isinstance(geom, Polygon):
            return geom if geom.area >= 1.0 else None
        
        # If it's a MultiPolygon, check each part
        if isinstance(geom, MultiPolygon):
            parts = [p for p in geom.geoms if p.area >= 1.0]
            if not parts:
                return None
            return MultiPolygon(parts)
        
        return geom

    # Apply filtering to geometries
    gdf['geometry'] = gdf['geometry'].apply(filter_small_parts)
    
    # Drop any rows that became empty/None after filtering parts
    gdf = gdf[gdf.geometry.notnull() & ~gdf.geometry.is_empty]
    
    # 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
    gdf = gdf.drop_duplicates(subset='geometry')

    # 8. Final column selection and cleanup
    # Required columns: parcel_id, parcel_class, district, area_m2
    # We'll also recalculate area_m2 based on the cleaned geometry to be sure, 
    # or just keep the column if it's expected to be the one from the record.
    # The prompt says "including these columns", implying they should be present.
    # Let's recalculate area_m2 to ensure it matches the cleaned geometry.
    gdf['area_m2'] = gdf['geometry'].area
    
    output_cols = ['parcel_id', 'parcel_class', 'district', 'area_m2', 'geometry']
    gdf = gdf[output_cols]

    # Save to geoparquet
    gdf.to_parquet(output_file)
    print(f"Successfully processed {len(gdf)} parcels.")

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

import geopandas as gpd

# Load the result
gdf = gpd.read_parquet('parcels_canonical.geoparquet')

print("Final info:")
print(gdf.info())
print("\nGeometry type counts:")
print(gdf.geom_type.value_counts())

# Check for duplicates in parcel_id
duplicates = gdf.duplicated(subset='parcel_id').sum()
print(f"\nDuplicate parcel_ids: {duplicates}")

# Check for null/empty geometries
null_empty = gdf.geometry.isna().sum() + gdf.geometry.is_empty.sum()
print(f"Null or empty geometries: {null_empty}")

# Check for invalid geometries
invalid = gdf.geometry.is_valid.sum()
print(f"Invalid geometries: {len(gdf) - invalid}")

# Check for duplicate geometries
geom_duplicates = gdf.duplicated(subset='geometry').sum()
print(f"Duplicate geometries: {geom_duplicates}")

# Check all are MultiPolygon
all_multipolygon = all(gdf.geom_type == 'MultiPolygon')
print(f"All are MultiPolygon: {all_multipolygon}")

# Check area_m2 matches geometry area
area_mismatch = (gdf['area_m2'] - gdf.geometry.area).abs().max()
print(f"Max area mismatch: {area_mismatch}")

# Check smallest area
print(f"Minimum area: {gdf['area_m2'].min()}")