103 lines
3.5 KiB
Python
103 lines
3.5 KiB
Python
#!/usr/bin/env python
|
|
"""
|
|
Debug script to diagnose why field boundary masking produces no data
|
|
"""
|
|
|
|
import json
|
|
import rasterio
|
|
from rasterio.mask import mask
|
|
from pathlib import Path
|
|
import numpy as np
|
|
import shapely.geometry as shgeom
|
|
|
|
# Load a sample field boundary
|
|
geojson_path = Path("laravel_app/storage/app/angata/Data/pivot.geojson")
|
|
field_id = "79" # A field that had issues
|
|
|
|
print(f"Loading field {field_id} from GeoJSON...")
|
|
with open(geojson_path) as f:
|
|
geojson_data = json.load(f)
|
|
|
|
field_boundary = None
|
|
for feature in geojson_data.get('features', []):
|
|
props = feature.get('properties', {})
|
|
if str(props.get('field', '')) == str(field_id):
|
|
geometry = feature.get('geometry')
|
|
if geometry:
|
|
geom_type = geometry.get('type', '')
|
|
coordinates = geometry.get('coordinates', [])
|
|
|
|
if geom_type == 'MultiPolygon':
|
|
if coordinates and len(coordinates) > 0:
|
|
coords = coordinates[0][0]
|
|
field_boundary = shgeom.Polygon(coords)
|
|
elif geom_type == 'Polygon':
|
|
if coordinates and len(coordinates) > 0:
|
|
coords = coordinates[0]
|
|
field_boundary = shgeom.Polygon(coords)
|
|
break
|
|
|
|
if field_boundary is None:
|
|
print(f"Field {field_id} not found")
|
|
exit(1)
|
|
|
|
print(f"Field boundary bounds: {field_boundary.bounds}")
|
|
print(f"Field boundary area: {field_boundary.area}")
|
|
|
|
# Load a sample TIFF tile
|
|
tiff_dir = Path("laravel_app/storage/app/angata/merged_final_tif/5x5")
|
|
tile_file = None
|
|
for date_dir in sorted(tiff_dir.iterdir()):
|
|
if date_dir.is_dir():
|
|
for tif in date_dir.glob("*.tif"):
|
|
if tif.stat().st_size > 12e6:
|
|
tile_file = tif
|
|
break
|
|
if tile_file:
|
|
break
|
|
|
|
if not tile_file:
|
|
print("No suitable TIFF found")
|
|
exit(1)
|
|
|
|
print(f"\nTesting with TIFF: {tile_file.name}")
|
|
|
|
with rasterio.open(tile_file) as src:
|
|
print(f"TIFF Bounds: {src.bounds}")
|
|
print(f"TIFF CRS: {src.crs}")
|
|
|
|
# Check if field boundary is within tile bounds
|
|
tile_box = shgeom.box(*src.bounds)
|
|
intersects = field_boundary.intersects(tile_box)
|
|
print(f"Field boundary intersects tile: {intersects}")
|
|
|
|
if intersects:
|
|
intersection = field_boundary.intersection(tile_box)
|
|
print(f"Intersection area: {intersection.area}")
|
|
print(f"Intersection bounds: {intersection.bounds}")
|
|
|
|
# Try to mask and see what we get
|
|
print("\nAttempting to mask...")
|
|
geom = shgeom.mapping(field_boundary)
|
|
try:
|
|
masked_data, _ = mask(src, [geom], crop=True, indexes=[1, 2, 3])
|
|
print(f"Masked data shape: {masked_data.shape}")
|
|
print(f"Masked data dtype: {masked_data.dtype}")
|
|
|
|
# Check the data
|
|
for i, band_idx in enumerate([1, 2, 3]):
|
|
band_data = masked_data[i]
|
|
print(f"\nBand {band_idx}:")
|
|
print(f" min: {np.nanmin(band_data):.6f}")
|
|
print(f" max: {np.nanmax(band_data):.6f}")
|
|
print(f" mean: {np.nanmean(band_data):.6f}")
|
|
print(f" % valid (non-zero): {(band_data != 0).sum() / band_data.size * 100:.2f}%")
|
|
print(f" % NaN: {np.isnan(band_data).sum() / band_data.size * 100:.2f}%")
|
|
|
|
# Show sample values
|
|
valid_pixels = band_data[band_data != 0]
|
|
if len(valid_pixels) > 0:
|
|
print(f" Sample valid values: {valid_pixels[:10]}")
|
|
except ValueError as e:
|
|
print(f"Error during masking: {e}")
|