#!/usr/bin/env python """ Test script to verify shapefile overlay functionality """ import geopandas as gpd import numpy as np from pathlib import Path # Test shapefile path shapefile_path = "ChauThanh/HienTrang/ChauThanh_kiemke.shp" print("=" * 70) print("TESTING SHAPEFILE OVERLAY") print("=" * 70) # Check if file exists shp = Path(shapefile_path) print(f"\n1. Checking file existence:") print(f" Path: {shp}") print(f" Exists: {shp.exists()}") print(f" Absolute: {shp.absolute()}") if shp.exists(): # Read shapefile print(f"\n2. Reading shapefile...") gdf = gpd.read_file(str(shp)) print(f" Features: {len(gdf)}") print(f" CRS: {gdf.crs}") print(f" Bounds: {gdf.total_bounds}") print(f" Columns: {list(gdf.columns)}") # Check geometries print(f"\n3. Checking geometries...") valid_count = sum(1 for geom in gdf.geometry if geom is not None and geom.is_valid) print(f" Valid geometries: {valid_count} / {len(gdf)}") # Sample geometry bounds if len(gdf) > 0: sample_geom = gdf.geometry.iloc[0] print(f" Sample geometry type: {sample_geom.geom_type}") print(f" Sample geometry bounds: {sample_geom.bounds}") # Test reprojection to EPSG:4326 print(f"\n4. Testing reprojection to EPSG:4326...") try: gdf_4326 = gdf.to_crs("EPSG:4326") print(f" Success!") print(f" New bounds: {gdf_4326.total_bounds}") except Exception as e: print(f" ERROR: {e}") # Test boundary extraction print(f"\n5. Testing boundary extraction...") boundaries = [] for geom in gdf.geometry: if geom is not None and geom.is_valid: boundary = geom.boundary if boundary is not None: boundaries.append(boundary) print(f" Extracted boundaries: {len(boundaries)}") # Test buffering print(f"\n6. Testing buffer...") buffer_size = 0.001 # degrees or meters depending on CRS buffered = [] for boundary in boundaries[:10]: # Test first 10 try: buf = boundary.buffer(buffer_size) buffered.append(buf) except Exception as e: print(f" Buffer error: {e}") print(f" Successfully buffered: {len(buffered)} / 10") else: print(" ERROR: Shapefile not found!") print("\n" + "=" * 70) print("TEST COMPLETE") print("=" * 70)