Files
remote-sensing/test_shapefile_overlay.py
T

80 lines
2.4 KiB
Python
Executable File

#!/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)