#!/usr/bin/env python3
"""
diag_glo30.py
=============
Diagnose why glo30_remote.vrt returns 0 at a known coordinate.

Usage:
    python diag_glo30.py
"""
import rasterio

LON, LAT = -115.923944, 35.4355  # WHB502, known elevation ~1359 m

print(f"Query point: lon={LON}, lat={LAT}")
print()

# --- 1. The VRT itself -------------------------------------------------
print("=== glo30_remote.vrt ===")
with rasterio.open("glo30_remote.vrt") as ds:
    print("CRS:", ds.crs)
    print("Bounds:", ds.bounds)
    print("Nodata:", ds.nodata)
    print("Res:", ds.res)
    print("Shape:", ds.width, ds.height)
    in_bounds = ds.bounds.left <= LON <= ds.bounds.right and ds.bounds.bottom <= LAT <= ds.bounds.top
    print("Point in bounds:", in_bounds)
    val = list(ds.sample([(LON, LAT)]))[0][0]
    print("Sampled value:", val)
print()

# --- 2. The specific source tile directly ------------------------------
url = "/vsicurl/https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N35_00_W116_00_DEM/Copernicus_DSM_COG_10_N35_00_W116_00_DEM.tif"
print("=== Direct source tile ===")
print(url)
try:
    with rasterio.open(url) as ds:
        print("CRS:", ds.crs)
        print("Bounds:", ds.bounds)
        print("Nodata:", ds.nodata)
        print("Res:", ds.res)
        in_bounds = ds.bounds.left <= LON <= ds.bounds.right and ds.bounds.bottom <= LAT <= ds.bounds.top
        print("Point in bounds:", in_bounds)
        val = list(ds.sample([(LON, LAT)]))[0][0]
        print("Sampled value:", val)
except Exception as e:
    print("ERROR opening tile directly:", e)