hypernets_example_band_integration.py

import matplotlib.pyplot as plt
import xarray as xr
import numpy as np
from matheo.band_integration import band_integration
from matheo.utils import function_def as fd
import obsarray

ds_HYP = xr.open_dataset("HYPERNETS_L_GHNA_L2A_REF_20231103T0901_20240124T2246_v2.0.nc")  # read digital effects table
reflectance_HYP=ds_HYP["reflectance"].values
u_reflectance_HYP=ds_HYP.unc["reflectance"].total_unc()
wavelength_HYP=ds_HYP["wavelength"].values

refl_S2, band_centres_S2 = band_integration.spectral_band_int_sensor(
    d=reflectance_HYP,
    wl=wavelength_HYP,
    platform_name="Sentinel-2A",
    sensor_name="MSI",
)

wav_sat = np.arange(400,1600,10)
width_sat = 10*np.ones_like(wav_sat)

refl_band = band_integration.pixel_int(
   d=reflectance_HYP,
   x=wavelength_HYP,
   x_pixel=wav_sat,
   width_pixel=width_sat,
   band_shape="triangle",
)

wav_SRF = np.arange(390,1610,0.1)
r_SRF = np.array([fd.f_triangle(wav_SRF, sat_wav_i, 10) for sat_wav_i in wav_sat])
refl_band2 = band_integration.band_int(reflectance_HYP, wavelength_HYP, r_SRF, wav_SRF)

plt.plot(wavelength_HYP,reflectance_HYP[:,30],label="HYPERNETS reflectance for index 30")
plt.fill_between(wavelength_HYP,reflectance_HYP[:,30]-u_reflectance_HYP[:,30],reflectance_HYP[:,30]+u_reflectance_HYP[:,30], alpha=0.3, label="HYPERNETS reflectance uncertainty")
plt.plot(wav_sat,refl_band[:,30],label="pixel_int with 10nm sampling and 10 nm width")
plt.plot(wav_sat,refl_band2[:,30],label="band_int with 10nm sampling and 10 nm width")
plt.plot(band_centres_S2,refl_S2[:,30],"s",label="spectral_band_int_sensor with S2A SRF")
plt.legend()
plt.ylim([0,0.6])
plt.show()