# Imagen satelital (Esri World Imagery, zoom 17 ≈ 1,2 m/px) reproyectada a la grilla UTM 19S del estudio (2 m/px) para el visor 3D.
# Solo visualización (no entra en el cálculo). Atribución: Esri, Maxar, Earthstar Geographics.
import math, io, numpy as np, urllib.request
from concurrent.futures import ThreadPoolExecutor
from PIL import Image
from pyproj import Transformer
from scipy import ndimage as nd
T = np.load('terreno.npz'); X0, Y0, DX, NX, NY = float(T['X0']), float(T['Y0']), float(T['DX']), int(T['NX']), int(T['NY'])
Z = 17; URL = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}'
inv = Transformer.from_crs('EPSG:32719', 'EPSG:4326', always_xy=True)
R = 2.0; W, H = int(NX * DX / R), int(NY * DX / R)
gx = X0 + (np.arange(W) + 0.5) * R; gy = Y0 + NY * DX - (np.arange(H) + 0.5) * R      # fila 0 = norte
GX, GY = np.meshgrid(gx, gy); LON, LAT = inv.transform(GX, GY)
n = 2 ** Z
fx = (LON + 180) / 360 * n; fy = (1 - np.log(np.tan(np.radians(LAT)) + 1 / np.cos(np.radians(LAT))) / math.pi) / 2 * n
tx0, tx1, ty0, ty1 = int(fx.min()), int(fx.max()), int(fy.min()), int(fy.max())
def grab(t):
    tx, ty = t
    req = urllib.request.Request(URL.format(z=Z, x=tx, y=ty), headers={'User-Agent': 'estudio-aluvion/1.0'})
    try: return t, Image.open(io.BytesIO(urllib.request.urlopen(req, timeout=30).read())).convert('RGB')
    except Exception: return t, Image.new('RGB', (256, 256), (90, 90, 90))
mos = Image.new('RGB', ((tx1 - tx0 + 1) * 256, (ty1 - ty0 + 1) * 256))
tiles = [(x, y) for y in range(ty0, ty1 + 1) for x in range(tx0, tx1 + 1)]
with ThreadPoolExecutor(12) as ex:
    for (x, y), im in ex.map(grab, tiles): mos.paste(im, ((x - tx0) * 256, (y - ty0) * 256))
M = np.asarray(mos).astype(np.float32); px = (fx - tx0) * 256 - 0.5; py = (fy - ty0) * 256 - 0.5
out = np.stack([nd.map_coordinates(M[..., c], [py, px], order=1) for c in range(3)], -1).clip(0, 255).astype(np.uint8)
Image.fromarray(out).save('web/satelite.jpg', quality=86); print('satélite', out.shape, 'teselas', len(tiles))