"""Tutorial 19 - RGB, Sentinel-2A and hyperspectral images of one maize field."""
import numpy as np
import less

scene = less.Scene()
scene.size = 10.0
scene.terrain = less.Terrain(property=less.Lambertian(reflectance=0.15))
scene.illumination = less.Illumination(
    source=less.Sun(zenith=30, azimuth=150), atmosphere=less.SimpleSpectralAtmosphere(turbidity=2.0))
maize = less.Object("maize", mesh=less.asset_path("maize.obj"))
maize.set_property(less.Prospect(cab=45, car=10, cw=0.012, cm=0.006, N=1.55))
positions = less.place.grid(scene, spacing=(0.75, 0.40), jitter=0.04, margin=0.3, seed=42)
rng = np.random.default_rng(42)
scene.add(maize, positions=positions, scales=rng.uniform(0.9, 1.1, len(positions)),
          rotations=rng.uniform(0, 360, len(positions)))

scene.simulate(less.OpticalImager(less.Orthographic(image_size=128), bands=[650, 550, 450],
                                  quality=16)).save("19_rgb.png")

bands = less.Sentinel2A()
s2 = scene.simulate(less.OpticalImager(less.Orthographic(image_size=128), bands=bands,
                                       spectral_resolution=20, quality=16))
s2.save("19_sentinel2a.tif")
brf = s2.to_brf().data
red, nir = brf[:, :, bands.band_names.index("B4")], brf[:, :, bands.band_names.index("B8")]
print("Mean NDVI:", round(float(((nir - red) / (nir + red)).mean()), 3))

edge_bands = less.Hyperspectral(start=650, stop=800, step=2)
edge = scene.simulate(less.OpticalImager(less.Orthographic(image_size=64), bands=edge_bands, quality=32))
wavelengths = edge_bands.get_band_centers()
spectrum = edge.to_brf().data.mean(axis=(0, 1))
print("Red-edge position:", wavelengths[np.argmax(np.gradient(spectrum, wavelengths))], "nm")