-
Notifications
You must be signed in to change notification settings - Fork 0
/
hypernets_example_band_integration.py
42 lines (36 loc) · 1.59 KB
/
hypernets_example_band_integration.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
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()