"""Tutorial 09 — Multispectral and Hyperspectral Simulation

RGB, Sentinel-2A, hyperspectral simulation of cornfield scene, NDVI calculation."""

import less
import numpy as np

# ── Constructing a corn field ───────────────────────────────────────────
scene = less.Scene()
scene.size = 10.0
scene.repetitive = False

scene.terrain = less.Terrain(property=less.Lambertian(reflectance=0.15))
scene.illumination = less.Illumination(source=less.Sun(zenith=30, azimuth=150), atmosphere=less.NoAtmosphere())

maize = less.Object("maize", mesh=less.examples.asset_path("maize.obj"))
maize.set_property(less.Prospect(cab=45, car=10, cw=0.012, cm=0.006, N=1.55))

rng = np.random.RandomState(42)
margin, row_sp, plant_sp = 0.5, 0.75, 0.40
n_rows  = int((10.0 - 2 * margin) / row_sp) + 1
n_plant = int((10.0 - 2 * margin) / plant_sp) + 1
xs = margin + np.arange(n_rows) * row_sp
ys = margin + np.arange(n_plant) * plant_sp
gx, gy = np.meshgrid(xs, ys, indexing='ij')
positions = np.column_stack([
    gx.ravel() + rng.uniform(-0.04, 0.04, gx.size),
    gy.ravel() + rng.uniform(-0.04, 0.04, gy.size),
    np.zeros(gx.size),
])
scene.add(maize, positions=positions,
          scales=rng.uniform(0.9, 1.1, len(positions)),
          rotations=rng.uniform(0, 360, len(positions)))
img_rgb = scene.simulate(less.OpticalImager(
    less.Orthographic(image_size=512),
    bands=[650, 550, 450], quality=128,
))
img_rgb.save("09_rgb.png")
print("RGB image saved")

# ── Sentinel-2A ──────────────────────────────────────────────
img_s2 = scene.simulate(less.OpticalImager(
    less.Orthographic(image_size=512),
    bands=less.Sentinel2A(), spectral_resolution=5, quality=128,
))
img_s2.save("09_sentinel2a.tif")

data_s2 = img_s2.data
ndvi_s2 = (data_s2[:,:,7] - data_s2[:,:,3]) / (data_s2[:,:,7] + data_s2[:,:,3] + 1e-10)
print(f"Sentinel-2A NDVI mean: {np.nanmean(ndvi_s2):.3f}")

# ── Hyperspectral ──────────────────────────────────────────────
img_hyper = scene.simulate(less.OpticalImager(
    less.Orthographic(image_size=256),
    bands=less.Hyperspectral(start=400, stop=1000, step=5),
    quality=64,
))
img_hyper.save("09_hyperspectral.tif")

mean_spec = np.mean(img_hyper.data, axis=(0, 1))
print(f"Hyperspectral bands: {mean_spec.shape[0]}")

print("Tutorial 09 Complete!")