108 lines
3.9 KiB
Python
108 lines
3.9 KiB
Python
#!/usr/bin/env python
|
|
"""
|
|
Debug script to find all tiles for a date and check which overlap with field boundary
|
|
"""
|
|
|
|
import json
|
|
import rasterio
|
|
from rasterio.mask import mask
|
|
from pathlib import Path
|
|
import numpy as np
|
|
import shapely.geometry as shgeom
|
|
import pandas as pd
|
|
|
|
# Load field 79 boundary
|
|
geojson_path = Path("laravel_app/storage/app/angata/Data/pivot.geojson")
|
|
field_id = "79"
|
|
|
|
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}")
|
|
|
|
# Find a specific date directory
|
|
tiff_dir = Path("laravel_app/storage/app/angata/merged_final_tif/5x5")
|
|
target_date = pd.Timestamp("2026-01-15") # Use a recent date that exists
|
|
|
|
# Find tiles for that date
|
|
date_dirs = []
|
|
for date_dir in tiff_dir.iterdir():
|
|
if date_dir.is_dir():
|
|
try:
|
|
dir_name = date_dir.name
|
|
date_str = dir_name.split('_')[0]
|
|
tile_date = pd.Timestamp(date_str)
|
|
if tile_date == target_date:
|
|
date_dirs.append(date_dir)
|
|
except:
|
|
pass
|
|
|
|
if not date_dirs:
|
|
print(f"No tiles found for {target_date}")
|
|
exit(1)
|
|
|
|
print(f"\nFound {len(date_dirs)} date directory(ies) for {target_date}")
|
|
|
|
for date_dir in date_dirs:
|
|
print(f"\n=== Checking date directory: {date_dir.name} ===")
|
|
|
|
tiles = list(date_dir.glob("*.tif"))
|
|
print(f"Found {len(tiles)} tiles in this directory")
|
|
|
|
for tile_path in sorted(tiles):
|
|
try:
|
|
with rasterio.open(tile_path) as src:
|
|
tile_bounds = src.bounds
|
|
tile_geom = shgeom.box(*tile_bounds)
|
|
|
|
intersects = field_boundary.intersects(tile_geom)
|
|
intersection = field_boundary.intersection(tile_geom) if intersects else None
|
|
intersection_area = intersection.area if intersection else 0
|
|
|
|
print(f"\n{tile_path.name}")
|
|
print(f" Tile bounds: {tile_bounds}")
|
|
print(f" Intersects field: {intersects}")
|
|
if intersects:
|
|
print(f" Intersection area: {intersection_area:.8f}")
|
|
|
|
# Try to mask this tile
|
|
geom = shgeom.mapping(field_boundary)
|
|
try:
|
|
masked_data, _ = mask(src, [geom], crop=True, indexes=[1, 2, 3])
|
|
print(f" ✓ Successfully masked! Shape: {masked_data.shape}")
|
|
|
|
# Check the data in each band
|
|
for i, band_idx in enumerate([1, 2, 3]):
|
|
band_data = masked_data[i]
|
|
non_zero = (band_data != 0).sum()
|
|
print(f" Band {band_idx}: {non_zero} non-zero pixels out of {band_data.size}")
|
|
except ValueError as e:
|
|
print(f" ✗ Masking failed: {e}")
|
|
except Exception as e:
|
|
print(f" Error reading tile: {e}")
|