dc-l2-cairo-invalid-dedup

Rerun

import geopandas as gpd

def analyze_geojson(file_path):
    try:
        gdf = gpd.read_file(file_path)
        print(f"Total features: {len(gdf)}")
        print("Columns:", gdf.columns.tolist())
        print("CRS:", gdf.crs)
        print("\nFirst few rows:")
        print(gdf.head())
        print("\nGeometry types:")
        print(gdf.geometry.type.value_counts())
        print("\nNull geometries:")
        print(gdf.geometry.isna().sum())
        
        # Check for invalid geometries
        invalid_count = gdf.geometry.is_valid.value_counts()
        print("\nValidity check:")
        print(invalid_count)

    except Exception as e:
        print(f"Error: {e}")

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

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

def process_parcels(input_path, output_path):
    # Load the data
    gdf = gpd.read_file(input_path)
    print(f"Initial features: {len(gdf)}")

    # 1. Fix invalid geometries
    # Using buffer(0) is a common trick to fix invalid geometries
    gdf['geometry'] = gdf.geometry.buffer(0)
    
    # Drop any rows that are still invalid or have null geometry
    gdf = gdf[gdf.geometry.is_valid & gdf.geometry.notna()]
    print(f"After fixing invalid/null geometries: {len(gdf)}")

    # 2. Handle duplicates: Keep record with earliest record_seq per parcel_id
    # Sort by parcel_id and record_seq (ascending)
    gdf = gdf.sort_values(by=['parcel_id', 'record_seq'])
    # Keep first record for each parcel_id
    gdf = gdf.drop_duplicates(subset=['parcel_id'], keep='first')
    print(f"After removing duplicate parcel_ids: {len(gdf)}")

    # 3. Remove polygon fragments smaller than 1 m²
    # Note: The input is in EPSG:22992 which is likely meters (Gauss-Kruger or similar)
    # We calculate area from the geometry itself.
    gdf['calc_area'] = gdf.geometry.area
    gdf = gdf[gdf['calc_area'] >= 1.0]
    print(f"After removing fragments <<  1m2: {len(gdf)}")

    # 4. Ensure every feature is a MultiPolygon
    def ensure_multipolygon(geom):
        if geom.geom_type == 'Polygon':
            return MultiPolygon([geom])
        return geom

    gdf['geometry'] = gdf.geometry.apply(ensure_multipolygon)
    print(f"After converting to MultiPolygon: {len(gdf)}")

    # 5. Remove duplicate geometries
    # Note: We drop duplicates based on the geometry itself to remove overlapping/identical artifacts
    gdf = gdf.drop_duplicates(subset=['geometry'])
    print(f"After removing duplicate geometries: {len(gdf)}")

    # 6. Final Cleanup: Select requested columns and ensure no nulls in them
    # Requested columns: parcel_id, parcel_class, district, area_m2
    # Note: area_m2 was in the original, but we should probably use our calculated one 
    # if the task implies 'area_m2' column in the output should be the actual area.
    # The prompt says "including these columns... area_m2". 
    # I will use the original column if it exists, but I'll verify.
    
    # Let's make sure we use the recalculated area for the area_m2 column to be accurate.
    gdf['area_m2'] = gdf['calc_area']
    
    output_cols = ['parcel_id', 'parcel_class', 'district', 'area_m2', 'geometry']
    gdf = gdf[output_cols]

    # Remove any remaining rows with nulls in key columns
    gdf = gdf.dropna(subset=['parcel_id', 'parcel_class', 'district'])

    # 7. Save to GeoParquet
    gdf.to_parquet(output_path)
    print(f"Final feature count: {len(gdf)}")
    print(f"Saved to {output_path}")

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

import geopandas as gpd

def verify_output(file_path):
    try:
        gdf = gpd.read_parquet(file_path)
        print(f"Total features: {len(gdf)}")
        print("Columns:", gdf.columns.tolist())
        print("CRS:", gdf.crs)
        
        # Check MultiPolygon
        geom_types = gdf.geometry.type.value_counts()
        print("\nGeometry types:")
        print(geom_types)
        
        # Check for null/empty geometries
        print("\nNull geometries:", gdf.geometry.isna().sum())
        print("Empty geometries:", gdf.geometry.is_empty.sum())
        
        # Check validity
        print("Invalid geometries:", (not gdf.geometry.is_valid).sum())
        
        # Check duplicates in parcel_id
        print("Duplicate parcel_id:", gdf['parcel_id'].duplicated().sum())
        
        # Check duplicate geometries
        print("Duplicate geometries:", gdf.geometry.duplicated().sum())
        
        # Check area
        print("Min area:", gdf.geometry.area.min())
        print("Max area:", gdf.geometry.area.max())
        
        # Check requested columns
        requested = ['parcel_id', 'parcel_class', 'district', 'area_m2']
        missing = [c for c in requested if c not in gdf.columns]
        if missing:
            print(f"Missing columns: {missing}")
        else:
            print("All requested columns present.")

        # Check for nulls in requested columns
        for col in requested:
            print(f"Nulls in {col}: {gdf[col].isna().sum()}")

    except Exception as e:
        print(f"Error: {e}")

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

import geopandas as gpd

def verify_output(file_path):
    try:
        gdf = gpd.read_parquet(file_path)
        print(f"Total features: {len(gdf)}")
        print("Columns:", gdf.columns.tolist())
        print("CRS:", gdf.crs)
        
        # Check MultiPolygon
        geom_types = gdf.geometry.type.value_counts()
        print("\nGeometry types:")
        print(geom_types)
        
        # Check for null/empty geometries
        print("\nNull geometries:", gdf.geometry.isna().sum())
        print("Empty geometries:", gdf.geometry.is_empty.sum())
        
        # Check validity
        print("Invalid geometries:", (gdf.geometry.is_valid == False).sum())
        
        # Check duplicates in parcel_id
        print("Duplicate parcel_id:", gdf['parcel_id'].duplicated().sum())
        
        # Check duplicate geometries
        print("Duplicate geometries:", gdf.geometry.duplicated().sum())
        
        # Check area
        print("Min area:", gdf.geometry.area.min())
        print("Max area:", gdf.geometry.area.max())
        
        # Check requested columns
        requested = ['parcel_id', 'parcel_class', 'district', 'area_m2']
        missing = [c for c in requested if c not in gdf.columns]
        if missing:
            print(f"Missing columns: {missing}")
        else:
            print("All requested columns present.")

        # Check for nulls in requested columns
        for col in requested:
            print(f"Nulls in {col}: {gdf[col].isna().sum()}")

    except Exception as e:
        print(f"Error: {e}")

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