Commit 0a070687 authored by Niklas Bohn's avatar Niklas Bohn
Browse files

Implemented option to choose between two solar irradiance models: 'new_kurucz' and 'fontenla'.


Signed-off-by: Niklas Bohn's avatarnbohn <nbohn@gfz-potsdam.de>
parent 96bb8f83
......@@ -99,7 +99,8 @@ setup(
"sicor/sensors/S2MSI/data/S2A_SNR_model.xlsx",
"sicor/AC/data/k_liquid_water_ice.xlsx",
"sicor/AC/data/newkur_EnMAP.dat",
"sicor/AC/data/solar_irradiances_400_2500_1.dill"
"sicor/AC/data/fontenla_EnMAP.dat",
"sicor/AC/data/conv_fac_fon_lut.dill"
])],
keywords=["SICOR", "EnMAP", "EnMAP-Box", "hyperspectral", "remote sensing", "satellite", "atmospheric correction"],
include_package_data=True,
......
......@@ -41,8 +41,8 @@ import platform
from py_tools_ds.processing.progress_mon import ProgressBar
from sicor.Tools.EnMAP.metadata import varsol
from sicor.Tools.EnMAP.LUT import get_data_file, read_lut_enmap_formatted, download_LUT, interpol_lut_c
from sicor.Tools.EnMAP.LUT import reduce_lut, interpol_lut_c_red
from sicor.Tools.EnMAP.LUT import get_data_file, read_lut_enmap_formatted, download_LUT, interpol_lut
from sicor.Tools.EnMAP.LUT import reduce_lut, interpol_lut_red
from sicor.Tools.EnMAP.conversion import generate_filter, table_to_array
from sicor.Tools.EnMAP.optimal_estimation import invert_function
from sicor.Tools.EnMAP.multiprocessing import SharedNdarray, initializer, mp_progress_bar
......@@ -284,14 +284,8 @@ class FoGen(object):
# LUT interpolation
vtest = np.asarray([pt[0], pt[1], pt[2], pt[3], pt[4], xx[0]])
f_int = interpol_lut_c(lut1=self.lut1_fit,
lut2=self.lut2_fit,
xnodes=self.xnodes,
nm_nodes=self.nm_nodes,
ndim=self.ndim,
x_cell=self.x_cell,
vtest=vtest,
intp_wvl=self.wvl_sel)
f_int = interpol_lut(lut1=self.lut1_fit, lut2=self.lut2_fit, xnodes=self.xnodes, nm_nodes=self.nm_nodes,
ndim=self.ndim, x_cell=self.x_cell, vtest=vtest, intp_wvl=self.wvl_sel)
f_int_l0 = f_int[0, :] * 1.e+3
f_int_edir = f_int[1, :] * 1.e+3
......@@ -326,23 +320,11 @@ class FoGen(object):
vtest = np.asarray([pt[0], pt[1], pt[2], pt[3], pt[4], xx[0]])
if mode == "full":
f_int = interpol_lut_c(lut1=self.lut1_full,
lut2=self.lut2_full,
xnodes=self.xnodes,
nm_nodes=self.nm_nodes,
ndim=self.ndim,
x_cell=self.x_cell,
vtest=vtest,
intp_wvl=self.wvl)
f_int = interpol_lut(lut1=self.lut1_full, lut2=self.lut2_full, xnodes=self.xnodes, nm_nodes=self.nm_nodes,
ndim=self.ndim, x_cell=self.x_cell, vtest=vtest, intp_wvl=self.wvl)
else:
f_int = interpol_lut_c(lut1=self.lut1_fit,
lut2=self.lut2_fit,
xnodes=self.xnodes,
nm_nodes=self.nm_nodes,
ndim=self.ndim,
x_cell=self.x_cell,
vtest=vtest,
intp_wvl=self.wvl_sel)
f_int = interpol_lut(lut1=self.lut1_fit, lut2=self.lut2_fit, xnodes=self.xnodes, nm_nodes=self.nm_nodes,
ndim=self.ndim, x_cell=self.x_cell, vtest=vtest, intp_wvl=self.wvl_sel)
f_int_l0 = f_int[0, :] * 1.e+3 * self.fac
f_int_edir = f_int[1, :] * 1.e+3 * self.fac
......@@ -606,23 +588,43 @@ class Fo(object):
self.wvl_swir = np.array(enmap_l1b.swir.detector_meta.wvl_center[:])
self.fwhm_swir = np.array(enmap_l1b.swir.detector_meta.fwhm[:])
# get solar irradiances for absorption feature shoulders wavelengths
self.logger.info("Getting solar irradiances for absorption feature shoulders wavelengths...")
new_kur_fn = get_data_file(module_name="sicor", file_basename="newkur_EnMAP.dat")
new_kur = pd.read_csv(new_kur_fn, delim_whitespace=True)
freq_0 = 10e6 / self.wvl_sel[0]
freq_1 = 10e6 / self.wvl_sel[-1]
solar_0 = np.zeros(1)
solar_1 = np.zeros(1)
# load solar irradiance model
self.logger.info("Loading solar irradiance model...")
self.sol_model = options["retrieval"]["sol_model"]
for ii in range(1, new_kur.shape[0]):
if int(freq_0) - int(new_kur.at[ii, "FREQ"]) == 0:
solar_0 = new_kur.at[ii, "SOLAR"]
if int(freq_1) - int(new_kur.at[ii, "FREQ"]) == 0:
solar_1 = new_kur.at[ii, "SOLAR"]
try:
assert self.sol_model == "new_kurucz" or self.sol_model == "fontenla"
except AssertionError:
raise AssertionError("No valid solar model is provided. Please indicate either 'new_kurucz' or 'fontenla' "
"in the options file.")
self.s0 = float(solar_0) * (10e6 / self.wvl_sel[0]) ** 2
self.s1 = float(solar_1) * (10e6 / self.wvl_sel[-1]) ** 2
if self.sol_model == "new_kurucz":
new_kur_fn = get_data_file(module_name="sicor", file_basename="newkur_EnMAP.dat")
new_kur = pd.read_csv(new_kur_fn, delim_whitespace=True)
new_kur_wvl = 10e6 / np.asarray(new_kur["FREQ"][1:], dtype=float)
new_kur_sol = np.asarray(new_kur["SOLAR"][1:], dtype=float) * (10e6 / new_kur_wvl) ** 2
self.solar_lut = {"wvl": new_kur_wvl, "sol_irr": new_kur_sol}
# get solar irradiances for shoulders of 1140 nm water absorption feature
s_norm_new_kur_fit = generate_filter(wvl_m=new_kur_wvl, wvl=self.wvl_sel, wl_resol=self.fwhm_sel)
new_kur_sol_fit = new_kur_sol @ s_norm_new_kur_fit
self.s0 = new_kur_sol_fit[0]
self.s1 = new_kur_sol_fit[-1]
elif self.sol_model == "fontenla":
fon_fn = get_data_file(module_name="sicor", file_basename="fontenla_EnMAP.dat")
fon = pd.read_csv(fon_fn, delim_whitespace=True)
fon_wvl = np.asarray(fon["um=micron"], dtype=float) * 1000
fon_sol = np.asarray(fon["E0"], dtype=float) * 10
self.solar_lut = {"wvl": fon_wvl, "sol_irr": fon_sol}
# get solar irradiances for shoulders of 1140 nm water absorption feature
s_norm_fon_fit = generate_filter(wvl_m=fon_wvl, wvl=self.wvl_sel, wl_resol=self.fwhm_sel)
new_kur_sol_fit = fon_sol @ s_norm_fon_fit
self.s0 = new_kur_sol_fit[0]
self.s1 = new_kur_sol_fit[-1]
# load RT LUT
# check if LUT file is available
......@@ -655,6 +657,33 @@ class Fo(object):
self.ndim = ndim
self.x_cell = x_cell
# applying spectral conversion factor to LUT in case Fontenla solar model is to be used
if self.sol_model == "fontenla":
self.logger.info("Converting LUT to Fontenla solar model...")
path_conv_fac = get_data_file(module_name="sicor", file_basename="conv_fac_fon_lut.dill")
with open(path_conv_fac, "rb") as fl:
conv_fac_fon_lut = dill.load(fl)
lut1_fon = np.zeros(lut1.shape)
lut2_fon = np.zeros(lut2.shape)
for ii in range(5):
for jj in range(6):
for kk in range(4):
for ll in range(6):
for mm in range(7):
lut1_fon[ii, jj, kk, ll, mm, 0, :, 0] = lut1[ii, jj, kk, ll, mm, 0, :, 0] * \
conv_fac_fon_lut
for nn in range(2):
lut2_fon[ii, jj, kk, ll, 0, mm, :, nn] = lut2[ii, jj, kk, ll, 0, mm, :, nn] * \
conv_fac_fon_lut
for ii in range(2):
lut2_fon[:, :, :, :, :, :, :, ii + 2] = lut2[:, :, :, :, :, :, :, ii + 2]
lut1 = lut1_fon
lut2 = lut2_fon
# resampling LUT to EnMAP wavelengths
self.logger.info("Resampling LUT to EnMAP wavelengths...")
......@@ -715,35 +744,18 @@ class Fo(object):
self.abs_co_w = 4 * np.pi * self.kw / self.wvl_sel
self.abs_co_i = 4 * np.pi * self.ki / self.wvl_sel
# calculate mean solar exoatmospheric irradiances
self.logger.info("Calculating mean solar exoatmospheric irradiances...")
path_sol = get_data_file(module_name="sicor", file_basename="solar_irradiances_400_2500_1.dill")
with open(path_sol, "rb") as fl:
solar_lut = dill.load(fl)
self.solar_res = solar_lut @ self.s_norm_fit
# perform first guess water vapor retrieval based on a common band ratio using VNIR data
warnings.filterwarnings("ignore")
if self.use_prior_mean[0]:
self.cwv_fg = np.full(self.data_vnir.shape[:2], self.prior_mean[0])
else:
logger.info("Performing first guess water vapor retrieval based on a common band ratio using VNIR data...")
self.cwv_fg = wv_band_ratio(data=self.data_vnir,
water_msk=self.water_mask_vnir,
fn_table=self.fn_table,
vza=self.pt[:, :, 0].mean(),
sza=self.pt[:, :, 1].mean(),
dem=self.hsf["vnir"].mean(),
aot=self.pt[:, :, 3].mean(),
raa=self.pt[:, :, 4].mean(),
intp_wvl=self.wvl_vnir,
intp_fwhm=self.fwhm_vnir,
jday=self.jday,
month=self.month,
idx=[73, 76, 81],
disable=self.disable_progressbars)
self.cwv_fg = wv_band_ratio(data=self.data_vnir, water_msk=self.water_mask_vnir, fn_table=self.fn_table,
sol_model=self.solar_lut, vza=self.pt[:, :, 0].mean(),
sza=self.pt[:, :, 1].mean(), dem=self.hsf["vnir"].mean(),
aot=self.pt[:, :, 3].mean(), raa=self.pt[:, :, 4].mean(),
intp_wvl=self.wvl_vnir, intp_fwhm=self.fwhm_vnir, jday=self.jday,
month=self.month, idx=[73, 76, 81], disable=self.disable_progressbars)
self.cwv_fg[np.logical_or(np.isnan(self.cwv_fg), np.isinf(self.cwv_fg))] = self.prior_mean[0]
self.cwv_fg[self.water_mask_vnir != 1] = self.prior_mean[0]
......@@ -752,8 +764,7 @@ class Fo(object):
logger.info("Transforming CWV first guess map to SWIR sensor geometry to enable segmentation and 3 phases of "
"water retrieval...")
self.cwv_fg_trans = enmap_l1b.transform_vnir_to_swir_raster(array_vnirsensorgeo=self.cwv_fg,
resamp_alg=self.resamp_alg,
respect_keystone=False)
resamp_alg=self.resamp_alg, respect_keystone=False)
self.cwv_fg_trans[self.cwv_fg_trans == 0.0] = self.prior_mean[0]
......@@ -815,8 +826,7 @@ class Fo(object):
if self.segmentation:
self.logger.info("Segmenting SWIR L1B spectra to enhance processing speed...")
self.X, self.segs, self.labels = SLIC_segmentation(data_rad_all=self.data_swir,
n_pca=self.n_pca,
self.X, self.segs, self.labels = SLIC_segmentation(data_rad_all=self.data_swir, n_pca=self.n_pca,
segs=self.segs)
# prepare segmented SWIR L1B data cube
......@@ -886,9 +896,9 @@ class Fo(object):
# LUT interpolation
vtest = np.asarray([pt[2], xx[0]])
f_int = interpol_lut_c_red(lut1=self.lut_red_fit[:, :, :, 0], lut2=self.lut_red_fit[:, :, :, 1:],
xnodes=self.xnodes_red, nm_nodes=self.nm_nodes_red, ndim=self.ndim_red,
x_cell=self.x_cell_red, vtest=vtest, intp_wvl=self.wvl_sel)
f_int = interpol_lut_red(lut1=self.lut_red_fit[:, :, :, 0], lut2=self.lut_red_fit[:, :, :, 1:],
xnodes=self.xnodes_red, nm_nodes=self.nm_nodes_red, ndim=self.ndim_red,
x_cell=self.x_cell_red, vtest=vtest, intp_wvl=self.wvl_sel)
f_int_l0 = f_int[0, :] * 1.e+3
f_int_edir = f_int[1, :] * 1.e+3
......@@ -923,17 +933,17 @@ class Fo(object):
vtest = np.asarray([pt[2], xx])
if mode == "vnir":
f_int = interpol_lut_c_red(lut1=self.lut_red_vnir[:, :, :, 0], lut2=self.lut_red_vnir[:, :, :, 1:],
xnodes=self.xnodes_red, nm_nodes=self.nm_nodes_red, ndim=self.ndim_red,
x_cell=self.x_cell_red, vtest=vtest, intp_wvl=self.wvl_vnir)
f_int = interpol_lut_red(lut1=self.lut_red_vnir[:, :, :, 0], lut2=self.lut_red_vnir[:, :, :, 1:],
xnodes=self.xnodes_red, nm_nodes=self.nm_nodes_red, ndim=self.ndim_red,
x_cell=self.x_cell_red, vtest=vtest, intp_wvl=self.wvl_vnir)
elif mode == "swir":
f_int = interpol_lut_c_red(lut1=self.lut_red_swir[:, :, :, 0], lut2=self.lut_red_swir[:, :, :, 1:],
xnodes=self.xnodes_red, nm_nodes=self.nm_nodes_red, ndim=self.ndim_red,
x_cell=self.x_cell_red, vtest=vtest, intp_wvl=self.wvl_swir)
f_int = interpol_lut_red(lut1=self.lut_red_swir[:, :, :, 0], lut2=self.lut_red_swir[:, :, :, 1:],
xnodes=self.xnodes_red, nm_nodes=self.nm_nodes_red, ndim=self.ndim_red,
x_cell=self.x_cell_red, vtest=vtest, intp_wvl=self.wvl_swir)
else:
f_int = interpol_lut_c_red(lut1=self.lut_red_fit[:, :, :, 0], lut2=self.lut_red_fit[:, :, :, 1:],
xnodes=self.xnodes_red, nm_nodes=self.nm_nodes_red, ndim=self.ndim_red,
x_cell=self.x_cell_red, vtest=vtest, intp_wvl=self.wvl_sel)
f_int = interpol_lut_red(lut1=self.lut_red_fit[:, :, :, 0], lut2=self.lut_red_fit[:, :, :, 1:],
xnodes=self.xnodes_red, nm_nodes=self.nm_nodes_red, ndim=self.ndim_red,
x_cell=self.x_cell_red, vtest=vtest, intp_wvl=self.wvl_sel)
f_int_l0 = f_int[0, :] * 1.e+3 * self.fac
f_int_edir = f_int[1, :] * 1.e+3 * self.fac
......@@ -1353,25 +1363,11 @@ def __oe__(ii):
__retrieval_noise__[i1, i2, :, :] = np.full((__state_dim__, __state_dim__), np.nan)
__smoothing_error__[i1, i2, :, :] = np.full((__state_dim__, __state_dim__), np.nan)
else:
se = __calc_se__(xx=__xa__[i1, i2, :],
dt=__data__[i1, i2, :],
pt=__pt__[i1, i2, :],
sb=__sb__,
num_bd=len(__fit_wvl__),
snr=__snr__,
unknowns=__unknowns__)
res = __inv_func__(yy=__data__[i1, i2, :],
fparam=__pt__[i1, i2, :],
ll=__ll__,
ul=__ul__,
xa=__xa__[i1, i2, :],
sa=__sa__,
se=se,
gnform=__gnform__,
full=__full__,
maxiter=__maxiter__,
eps=__eps__)
se = __calc_se__(xx=__xa__[i1, i2, :], dt=__data__[i1, i2, :], pt=__pt__[i1, i2, :], sb=__sb__,
num_bd=len(__fit_wvl__), snr=__snr__, unknowns=__unknowns__)
res = __inv_func__(yy=__data__[i1, i2, :], fparam=__pt__[i1, i2, :], ll=__ll__, ul=__ul__, xa=__xa__[i1, i2, :],
sa=__sa__, se=se, gnform=__gnform__, full=__full__, maxiter=__maxiter__, eps=__eps__)
model = __forward__(xx=res[0], pt=__pt__[i1, i2, :], dt=__data__[i1, i2, :], model_output=True)
......
This diff is collapsed.
......@@ -238,19 +238,21 @@ def read_lut_enmap_formatted(file_lut):
return luts, axes_x, axes_y, wvl, lut1, lut2, xnodes, nm_nodes, ndim, x_cell
def interpol_lut_c(lut1, lut2, xnodes, nm_nodes, ndim, x_cell, vtest, intp_wvl):
def interpol_lut(lut1, lut2, xnodes, nm_nodes, ndim, x_cell, vtest, intp_wvl):
"""
Multidimensional LUT interpolation based on numba.
:param lut1: Path radiance LUT
:param lut2: Solar irradiance, spherical albedo and transmittance LUT
:param xnodes: Gridpoints for each LUT dimension
:param nm_nodes: Overall number of LUT gridpoints (2**ndim)
:param ndim: Dimension of LUT
:param x_cell: Interpolation grid (nm_nodes x n_dim)
:param vtest: Atmospheric state vector (interpolation point); contains VZA, SZA, HSF, AOT, RAA, CWV
Multidimensional LUT interpolation.
:param lut1: LUT containing path radiance
:param lut2: LUT containing downwelling direct and diffuse surface solar irradiance, spherical albedo and
transmittance
:param xnodes: gridpoints for each LUT dimension
:param nm_nodes: overall number of LUT gridpoints (2**ndim)
:param ndim: dimension of LUT
:param x_cell: interpolation grid (nm_nodes x n_dim)
:param vtest: atmospheric state vector (interpolation point); contains VZA, SZA, HSF, AOT, RAA, CWV
:param intp_wvl: wavelengths for which interpolation should be conducted
:return: Path radiance, solar irradiance, spherical albedo and transmittance interpolated at vtest
:return: path radiance, downwelling direct and diffuse surface solar irradiance, spherical albedo and
transmittance interpolated at vtest
"""
lim = np.zeros((2, ndim), np.int8)
......@@ -302,19 +304,21 @@ def interpol_lut_c(lut1, lut2, xnodes, nm_nodes, ndim, x_cell, vtest, intp_wvl):
return f_int
def interpol_lut_c_red(lut1, lut2, xnodes, nm_nodes, ndim, x_cell, vtest, intp_wvl):
def interpol_lut_red(lut1, lut2, xnodes, nm_nodes, ndim, x_cell, vtest, intp_wvl):
"""
LUT interpolation based on numba in a grduced grid.
:param lut1: Path radiance LUT
:param lut2: Solar irradiance, spherical albedo and transmittance LUT
:param xnodes: Gridpoints for each LUT dimension
:param nm_nodes: Overall number of LUT gridpoints (2**ndim)
:param ndim: Dimension of LUT
:param x_cell: Interpolation grid (nm_nodes x n_dim)
:param vtest: Atmospheric state vector (interpolation point); only surface elevation and CWV
LUT interpolation based in a reduced grid.
:param lut1: LUT containing path radiance
:param lut2: LUT containing downwelling direct and diffuse surface solar irradiance, spherical albedo and
transmittance
:param xnodes: gridpoints for each LUT dimension
:param nm_nodes: overall number of LUT gridpoints (2**ndim)
:param ndim: dimension of LUT
:param x_cell: interpolation grid (nm_nodes x n_dim)
:param vtest: atmospheric state vector (interpolation point); only surface elevation and CWV
:param intp_wvl: wavelengths for which interpolation should be conducted
:return: Path radiance, solar irradiance, spherical albedo and transmittance interpolated at vtest
:return: path radiance, downwelling direct and diffuse surface solar irradiance, spherical albedo and
transmittance interpolated at vtest
"""
lim = np.zeros((2, ndim), np.int8)
......@@ -358,17 +362,19 @@ def reduce_lut(lut1, lut2, xnodes, nm_nodes, ndim, x_cell, gp, intp_wvl):
"""
Reduce grid dimensionality of LUT based on fixed solar and view geometry as well as one single scene value of AOT.
:param lut1: Path radiance LUT
:param lut2: Solar irradiance, spherical albedo and transmittance LUT
:param xnodes: Gridpoints for each LUT dimension
:param nm_nodes: Overall number of LUT gridpoints (2**ndim)
:param ndim: Dimension of LUT
:param x_cell: Interpolation grid (nm_nodes x n_dim)
:param gp: State vector of fixed scene parameters (VZA, SZA, AOT, RAA)
:param lut1: LUT containing path radiance
:param lut2: LUT containing downwelling direct and diffuse surface solar irradiance, spherical albedo and
transmittance
:param xnodes: gridpoints for each LUT dimension
:param nm_nodes: overall number of LUT gridpoints (2**ndim)
:param ndim: dimension of LUT
:param x_cell: interpolation grid (nm_nodes x n_dim)
:param gp: state vector of fixed scene parameters (VZA, SZA, AOT, RAA)
:param intp_wvl: wavelengths for which interpolation should be conducted
:return: 2D LUT including path radiance, solar irradiance, spherical albedo and transmittance interpolated
at gp for varying surface elevation and CWV; reduced gridpoints for each LUT dimension; reduced
overall number of LUT gridpoints; reduced dimension of LUT; reduced interpolation grid
:return: 2D LUT including path radiance, downwelling direct and diffuse surface solar irradiance, spherical
albedo and transmittance interpolated at gp for varying surface elevation and CWV; reduced
gridpoints for each LUT dimension; reduced overall number of LUT gridpoints; reduced dimension of
LUT; reduced interpolation grid
"""
red_lut = np.zeros((len(xnodes[:4, 2]), len(xnodes[:, 5]), len(intp_wvl), lut2.shape[-1] + 1))
......@@ -376,8 +382,8 @@ def reduce_lut(lut1, lut2, xnodes, nm_nodes, ndim, x_cell, gp, intp_wvl):
for jj, cwv in enumerate(xnodes[:, 5]):
vtest = np.array([gp[0], gp[1], hsf, gp[2], gp[3], cwv])
pp = interpol_lut_c(lut1=lut1, lut2=lut2, xnodes=xnodes, nm_nodes=nm_nodes, ndim=ndim, x_cell=x_cell,
vtest=vtest, intp_wvl=intp_wvl)
pp = interpol_lut(lut1=lut1, lut2=lut2, xnodes=xnodes, nm_nodes=nm_nodes, ndim=ndim, x_cell=x_cell,
vtest=vtest, intp_wvl=intp_wvl)
red_lut[ii, jj, :, 0] = pp[0, :]
red_lut[ii, jj, :, 1] = pp[1, :]
......
......@@ -25,21 +25,21 @@
import numpy as np
import dill
from tqdm import tqdm
from sicor.Tools.EnMAP.LUT import read_lut_enmap_formatted, interpol_lut_c, get_data_file
from sicor.Tools.EnMAP.LUT import read_lut_enmap_formatted, interpol_lut
from sicor.Tools.EnMAP.conversion import generate_filter
from sicor.Tools.EnMAP.metadata import varsol
def wv_band_ratio(data, water_msk, fn_table, vza, sza, dem, aot, raa, intp_wvl, intp_fwhm, jday, month, idx,
def wv_band_ratio(data, water_msk, fn_table, sol_model, vza, sza, dem, aot, raa, intp_wvl, intp_fwhm, jday, month, idx,
disable=False):
"""Band ratio water vapor retrieval.
:param data: image dataset
:param water_msk: water mask
:param fn_table: path to radiative transfer LUT
:param sol_model: dictionary containing solar irradiance model ('wvl', 'sol_irr')
:param vza: viewing zenitg angle
:param sza: sun zenith angle
:param dem: digital elevation model, same shape as data
......@@ -67,14 +67,14 @@ def wv_band_ratio(data, water_msk, fn_table, vza, sza, dem, aot, raa, intp_wvl,
luts, axes_x, axes_y, wvl, lut1, lut2, xnodes, nm_nodes, ndim, x_cell = read_lut_enmap_formatted(file_lut=fn_table)
wvl_lut = wvl
s_norm = generate_filter(wvl_m=wvl_lut, wvl=intp_wvl, wl_resol=intp_fwhm)
s_norm_lut = generate_filter(wvl_m=wvl_lut, wvl=intp_wvl, wl_resol=intp_fwhm)
lut2_shape = np.array(lut2.shape)
lut2_shape[6] = len(intp_wvl)
lut2_res = np.zeros(lut2_shape)
lut1_res = lut1[:, :, :, :, :, :, :, 0] @ s_norm
lut1_res = lut1[:, :, :, :, :, :, :, 0] @ s_norm_lut
for ii in range(lut2.shape[-1]):
lut2_res[:, :, :, :, :, :, :, ii] = lut2[:, :, :, :, :, :, :, ii] @ s_norm
lut2_res[:, :, :, :, :, :, :, ii] = lut2[:, :, :, :, :, :, :, ii] @ s_norm_lut
dsol = varsol(jday, month)
dn2rad = dsol * dsol * 0.1
......@@ -89,8 +89,8 @@ def wv_band_ratio(data, water_msk, fn_table, vza, sza, dem, aot, raa, intp_wvl,
for jj, cwv in enumerate(cwvs):
vtest = np.asarray([vza, sza, hsf, aot, raa, cwv])
f_int = interpol_lut_c(lut1=lut1_res, lut2=lut2_res, xnodes=xnodes, nm_nodes=nm_nodes, ndim=ndim,
x_cell=x_cell, vtest=vtest, intp_wvl=intp_wvl)
f_int = interpol_lut(lut1=lut1_res, lut2=lut2_res, xnodes=xnodes, nm_nodes=nm_nodes, ndim=ndim,
x_cell=x_cell, vtest=vtest, intp_wvl=intp_wvl)
f_int_l0 = f_int[0, :] * 1.e+3
f_int_edir = f_int[1, :] * 1.e+3
......@@ -103,10 +103,8 @@ def wv_band_ratio(data, water_msk, fn_table, vza, sza, dem, aot, raa, intp_wvl,
l_toa = (f_int_l0 + f_int_ee * rho / np.pi / (1 - f_int_ss * rho)) * fac
l_toa_lut[ii, jj, kk, :] = l_toa
path_sol = get_data_file(module_name="sicor", file_basename="solar_irradiances_400_2500_1.dill")
with open(path_sol, "rb") as fl:
solar_lut = dill.load(fl)
solar_res = solar_lut @ s_norm
s_norm_sol = generate_filter(wvl_m=sol_model["wvl"], wvl=intp_wvl, wl_resol=intp_fwhm)
solar_res = sol_model["sol_irr"] @ s_norm_sol
rfl_1_img = data[:, :, idx[0]] / solar_res[idx[0]]
rfl_2_img = data[:, :, idx[1]] / solar_res[idx[1]]
......
......@@ -14,6 +14,7 @@
at coastlines); otherwise, all image pixels (land + water) are processed; default:
false*/
"fn_LUT": "",
"sol_model": "new_kurucz", /*available solar models: 'new_kurucz' (Kurucz 2005), 'fontenla' (Fontenla 2011)*/
"cpu": 32,
"disable_progressbars": false,
"segmentation": true,
......
......@@ -39,7 +39,7 @@ from enpt.options.config import EnPTConfig, config_for_testing
from sicor.options.options import get_options
from sicor.sicor_enmap import sicor_ac_enmap
from sicor.Tools.EnMAP.LUT import interpol_lut_c
from sicor.Tools.EnMAP.LUT import interpol_lut
from sicor.AC.RtFo_3_phases import Fo
......@@ -193,8 +193,8 @@ class TestSicor_EnMAP_interpol_LUT(TestCase):
print("################################################")
print("")
print("")
interpol_lut_c(lut1=self.lut1_fit, lut2=self.lut2_fit, xnodes=self.xnodes, nm_nodes=self.nm_nodes,
ndim=self.ndim, x_cell=self.x_cell, vtest=self.vtest, intp_wvl=self.wvl_sel)
interpol_lut(lut1=self.lut1_fit, lut2=self.lut2_fit, xnodes=self.xnodes, nm_nodes=self.nm_nodes,
ndim=self.ndim, x_cell=self.x_cell, vtest=self.vtest, intp_wvl=self.wvl_sel)
print("Done!")
print("")
print("")
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment