Commit 6496a718 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Merge branch 'master' into enhancement/improve_ecmwf_crawler

parents ec52e2a8 84b1efd1
...@@ -6,10 +6,12 @@ History ...@@ -6,10 +6,12 @@ History
-------------------- --------------------
New features: New features:
* * Multiprocessing mode of SICOR available again on macOS.
* Option to choose between two solar irradiance models: 'new_kurucz' and 'fontenla'.
Bugfixes: Bugfixes:
* * Added missing initializer to multiprocessing pool in empirical line calculation and set multiprocessing start method to fork.
* Disabled water vapor first guess retrieval over water surfaces in case SICOR is running for land+water pixels.
0.16.4 (2021-06-18) 0.16.4 (2021-06-18)
......
...@@ -99,7 +99,8 @@ setup( ...@@ -99,7 +99,8 @@ setup(
"sicor/sensors/S2MSI/data/S2A_SNR_model.xlsx", "sicor/sensors/S2MSI/data/S2A_SNR_model.xlsx",
"sicor/AC/data/k_liquid_water_ice.xlsx", "sicor/AC/data/k_liquid_water_ice.xlsx",
"sicor/AC/data/newkur_EnMAP.dat", "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"], keywords=["SICOR", "EnMAP", "EnMAP-Box", "hyperspectral", "remote sensing", "satellite", "atmospheric correction"],
include_package_data=True, include_package_data=True,
......
...@@ -28,7 +28,7 @@ import logging ...@@ -28,7 +28,7 @@ import logging
import numpy as np import numpy as np
import pandas as pd import pandas as pd
import dill import dill
from multiprocessing import Pool import multiprocessing as mp
from itertools import product from itertools import product
from time import time from time import time
import warnings import warnings
...@@ -41,8 +41,8 @@ import platform ...@@ -41,8 +41,8 @@ import platform
from py_tools_ds.processing.progress_mon import ProgressBar from py_tools_ds.processing.progress_mon import ProgressBar
from sicor.Tools.EnMAP.metadata import varsol 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 get_data_file, read_lut_enmap_formatted, download_LUT, interpol_lut
from sicor.Tools.EnMAP.LUT import reduce_lut, interpol_lut_c_red 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.conversion import generate_filter, table_to_array
from sicor.Tools.EnMAP.optimal_estimation import invert_function from sicor.Tools.EnMAP.optimal_estimation import invert_function
from sicor.Tools.EnMAP.multiprocessing import SharedNdarray, initializer, mp_progress_bar from sicor.Tools.EnMAP.multiprocessing import SharedNdarray, initializer, mp_progress_bar
...@@ -284,14 +284,8 @@ class FoGen(object): ...@@ -284,14 +284,8 @@ class FoGen(object):
# LUT interpolation # LUT interpolation
vtest = np.asarray([pt[0], pt[1], pt[2], pt[3], pt[4], xx[0]]) vtest = np.asarray([pt[0], pt[1], pt[2], pt[3], pt[4], xx[0]])
f_int = interpol_lut_c(lut1=self.lut1_fit, f_int = interpol_lut(lut1=self.lut1_fit, lut2=self.lut2_fit, xnodes=self.xnodes, nm_nodes=self.nm_nodes,
lut2=self.lut2_fit, ndim=self.ndim, x_cell=self.x_cell, vtest=vtest, intp_wvl=self.wvl_sel)
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_l0 = f_int[0, :] * 1.e+3
f_int_edir = f_int[1, :] * 1.e+3 f_int_edir = f_int[1, :] * 1.e+3
...@@ -326,23 +320,11 @@ class FoGen(object): ...@@ -326,23 +320,11 @@ class FoGen(object):
vtest = np.asarray([pt[0], pt[1], pt[2], pt[3], pt[4], xx[0]]) vtest = np.asarray([pt[0], pt[1], pt[2], pt[3], pt[4], xx[0]])
if mode == "full": if mode == "full":
f_int = interpol_lut_c(lut1=self.lut1_full, f_int = interpol_lut(lut1=self.lut1_full, lut2=self.lut2_full, xnodes=self.xnodes, nm_nodes=self.nm_nodes,
lut2=self.lut2_full, ndim=self.ndim, x_cell=self.x_cell, vtest=vtest, intp_wvl=self.wvl)
xnodes=self.xnodes,
nm_nodes=self.nm_nodes,
ndim=self.ndim,
x_cell=self.x_cell,
vtest=vtest,
intp_wvl=self.wvl)
else: else:
f_int = interpol_lut_c(lut1=self.lut1_fit, f_int = interpol_lut(lut1=self.lut1_fit, lut2=self.lut2_fit, xnodes=self.xnodes, nm_nodes=self.nm_nodes,
lut2=self.lut2_fit, ndim=self.ndim, x_cell=self.x_cell, vtest=vtest, intp_wvl=self.wvl_sel)
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_l0 = f_int[0, :] * 1.e+3 * self.fac
f_int_edir = f_int[1, :] * 1.e+3 * self.fac f_int_edir = f_int[1, :] * 1.e+3 * self.fac
...@@ -606,23 +588,43 @@ class Fo(object): ...@@ -606,23 +588,43 @@ class Fo(object):
self.wvl_swir = np.array(enmap_l1b.swir.detector_meta.wvl_center[:]) self.wvl_swir = np.array(enmap_l1b.swir.detector_meta.wvl_center[:])
self.fwhm_swir = np.array(enmap_l1b.swir.detector_meta.fwhm[:]) self.fwhm_swir = np.array(enmap_l1b.swir.detector_meta.fwhm[:])
# get solar irradiances for absorption feature shoulders wavelengths # load solar irradiance model
self.logger.info("Getting solar irradiances for absorption feature shoulders wavelengths...") self.logger.info("Loading solar irradiance model...")
new_kur_fn = get_data_file(module_name="sicor", file_basename="newkur_EnMAP.dat") self.sol_model = options["retrieval"]["sol_model"]
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)
for ii in range(1, new_kur.shape[0]): try:
if int(freq_0) - int(new_kur.at[ii, "FREQ"]) == 0: assert self.sol_model == "new_kurucz" or self.sol_model == "fontenla"
solar_0 = new_kur.at[ii, "SOLAR"] except AssertionError:
if int(freq_1) - int(new_kur.at[ii, "FREQ"]) == 0: raise AssertionError("No valid solar model is provided. Please indicate either 'new_kurucz' or 'fontenla' "
solar_1 = new_kur.at[ii, "SOLAR"] "in the options file.")
self.s0 = float(solar_0) * (10e6 / self.wvl_sel[0]) ** 2 if self.sol_model == "new_kurucz":
self.s1 = float(solar_1) * (10e6 / self.wvl_sel[-1]) ** 2 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 # load RT LUT
# check if LUT file is available # check if LUT file is available
...@@ -655,6 +657,33 @@ class Fo(object): ...@@ -655,6 +657,33 @@ class Fo(object):
self.ndim = ndim self.ndim = ndim
self.x_cell = x_cell 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 # resampling LUT to EnMAP wavelengths
self.logger.info("Resampling LUT to EnMAP wavelengths...") self.logger.info("Resampling LUT to EnMAP wavelengths...")
...@@ -715,45 +744,27 @@ class Fo(object): ...@@ -715,45 +744,27 @@ class Fo(object):
self.abs_co_w = 4 * np.pi * self.kw / self.wvl_sel self.abs_co_w = 4 * np.pi * self.kw / self.wvl_sel
self.abs_co_i = 4 * np.pi * self.ki / 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 # perform first guess water vapor retrieval based on a common band ratio using VNIR data
warnings.filterwarnings("ignore") warnings.filterwarnings("ignore")
if self.use_prior_mean[0]: if self.use_prior_mean[0]:
self.cwv_fg = np.full(self.data_vnir.shape[:2], self.prior_mean[0]) self.cwv_fg = np.full(self.data_vnir.shape[:2], self.prior_mean[0])
else: else:
logger.info("Performing first guess water vapor retrieval based on a common band ratio using VNIR data...") 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, self.cwv_fg = wv_band_ratio(data=self.data_vnir, water_msk=self.water_mask_vnir, fn_table=self.fn_table,
water_msk=self.water_mask_vnir, sol_model=self.solar_lut, vza=self.pt[:, :, 0].mean(),
land_only=self.land_only, sza=self.pt[:, :, 1].mean(), dem=self.hsf["vnir"].mean(),
fn_table=self.fn_table, aot=self.pt[:, :, 3].mean(), raa=self.pt[:, :, 4].mean(),
vza=self.pt[:, :, 0].mean(), intp_wvl=self.wvl_vnir, intp_fwhm=self.fwhm_vnir, jday=self.jday,
sza=self.pt[:, :, 1].mean(), month=self.month, idx=[73, 76, 81], disable=self.disable_progressbars)
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[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]
# transform CWV first guess map to SWIR sensor geometry to enable segmentation and 3 phases of water retrieval # transform CWV first guess map to SWIR sensor geometry to enable segmentation and 3 phases of water retrieval
logger.info("Transforming CWV first guess map to SWIR sensor geometry to enable segmentation and 3 phases of " logger.info("Transforming CWV first guess map to SWIR sensor geometry to enable segmentation and 3 phases of "
"water retrieval...") "water retrieval...")
self.cwv_fg_trans = enmap_l1b.transform_vnir_to_swir_raster(array_vnirsensorgeo=self.cwv_fg, self.cwv_fg_trans = enmap_l1b.transform_vnir_to_swir_raster(array_vnirsensorgeo=self.cwv_fg,
resamp_alg=self.resamp_alg, resamp_alg=self.resamp_alg, respect_keystone=False)
respect_keystone=False)
self.cwv_fg_trans[self.cwv_fg_trans == 0.0] = self.prior_mean[0] self.cwv_fg_trans[self.cwv_fg_trans == 0.0] = self.prior_mean[0]
...@@ -815,8 +826,7 @@ class Fo(object): ...@@ -815,8 +826,7 @@ class Fo(object):
if self.segmentation: if self.segmentation:
self.logger.info("Segmenting SWIR L1B spectra to enhance processing speed...") 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, self.X, self.segs, self.labels = SLIC_segmentation(data_rad_all=self.data_swir, n_pca=self.n_pca,
n_pca=self.n_pca,
segs=self.segs) segs=self.segs)
# prepare segmented SWIR L1B data cube # prepare segmented SWIR L1B data cube
...@@ -886,9 +896,9 @@ class Fo(object): ...@@ -886,9 +896,9 @@ class Fo(object):
# LUT interpolation # LUT interpolation
vtest = np.asarray([pt[2], xx[0]]) 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:], 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, 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) x_cell=self.x_cell_red, vtest=vtest, intp_wvl=self.wvl_sel)
f_int_l0 = f_int[0, :] * 1.e+3 f_int_l0 = f_int[0, :] * 1.e+3
f_int_edir = f_int[1, :] * 1.e+3 f_int_edir = f_int[1, :] * 1.e+3
...@@ -923,17 +933,17 @@ class Fo(object): ...@@ -923,17 +933,17 @@ class Fo(object):
vtest = np.asarray([pt[2], xx]) vtest = np.asarray([pt[2], xx])
if mode == "vnir": if mode == "vnir":
f_int = interpol_lut_c_red(lut1=self.lut_red_vnir[:, :, :, 0], lut2=self.lut_red_vnir[:, :, :, 1:], 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, 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) x_cell=self.x_cell_red, vtest=vtest, intp_wvl=self.wvl_vnir)
elif mode == "swir": elif mode == "swir":
f_int = interpol_lut_c_red(lut1=self.lut_red_swir[:, :, :, 0], lut2=self.lut_red_swir[:, :, :, 1:], 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, 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) x_cell=self.x_cell_red, vtest=vtest, intp_wvl=self.wvl_swir)
else: else:
f_int = interpol_lut_c_red(lut1=self.lut_red_fit[:, :, :, 0], lut2=self.lut_red_fit[:, :, :, 1:], 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, 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) 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_l0 = f_int[0, :] * 1.e+3 * self.fac
f_int_edir = f_int[1, :] * 1.e+3 * self.fac f_int_edir = f_int[1, :] * 1.e+3 * self.fac
...@@ -1127,9 +1137,7 @@ def __minimize__(fo, opt_func, unknowns=False, logger=None): ...@@ -1127,9 +1137,7 @@ def __minimize__(fo, opt_func, unknowns=False, logger=None):
if platform.system() == "Windows" and processes > 1: if platform.system() == "Windows" and processes > 1:
logger.warning('Multiprocessing is currently not available on Windows.') logger.warning('Multiprocessing is currently not available on Windows.')
if platform.system() == "Darwin" and processes > 1: if platform.system() == "Windows" or processes == 1:
logger.warning('Multiprocessing is currently not available on macOS.')
if platform.system() == "Windows" or platform.system() == "Darwin" or processes == 1:
logger.info("Singleprocessing on 1 cpu") logger.info("Singleprocessing on 1 cpu")
else: else:
logger.info("Setting up multiprocessing...") logger.info("Setting up multiprocessing...")
...@@ -1234,11 +1242,12 @@ def __minimize__(fo, opt_func, unknowns=False, logger=None): ...@@ -1234,11 +1242,12 @@ def __minimize__(fo, opt_func, unknowns=False, logger=None):
# check if operating system is 'Windows'; in that case, multiprocessing is currently not working # check if operating system is 'Windows'; in that case, multiprocessing is currently not working
# TODO: enable Windows compatibility for multiprocessing # TODO: enable Windows compatibility for multiprocessing
if platform.system() == "Windows" or platform.system() == "Darwin" or processes == 1: if platform.system() == "Windows" or processes == 1:
initializer(globals(), globs) initializer(globals(), globs)
[mp_fun(ii) for ii in tqdm(rng, disable=fo.disable_progressbars)] [mp_fun(ii) for ii in tqdm(rng, disable=fo.disable_progressbars)]
else: else:
with closing(Pool(processes=processes, initializer=initializer, initargs=(globals(), globs,))) as pl: with closing(mp.get_context("fork").Pool(processes=processes, initializer=initializer, initargs=(globals(),
globs))) as pl:
results = pl.map_async(mp_fun, rng, chunksize=1) results = pl.map_async(mp_fun, rng, chunksize=1)
if not fo.disable_progressbars: if not fo.disable_progressbars:
bar = ProgressBar(prefix='\tprogress:') bar = ProgressBar(prefix='\tprogress:')
...@@ -1354,25 +1363,11 @@ def __oe__(ii): ...@@ -1354,25 +1363,11 @@ def __oe__(ii):
__retrieval_noise__[i1, i2, :, :] = np.full((__state_dim__, __state_dim__), np.nan) __retrieval_noise__[i1, i2, :, :] = np.full((__state_dim__, __state_dim__), np.nan)
__smoothing_error__[i1, i2, :, :] = np.full((__state_dim__, __state_dim__), np.nan) __smoothing_error__[i1, i2, :, :] = np.full((__state_dim__, __state_dim__), np.nan)
else: else:
se = __calc_se__(xx=__xa__[i1, i2, :], se = __calc_se__(xx=__xa__[i1, i2, :], dt=__data__[i1, i2, :], pt=__pt__[i1, i2, :], sb=__sb__,
dt=__data__[i1, i2, :], num_bd=len(__fit_wvl__), snr=__snr__, unknowns=__unknowns__)
pt=__pt__[i1, i2, :],
sb=__sb__, res = __inv_func__(yy=__data__[i1, i2, :], fparam=__pt__[i1, i2, :], ll=__ll__, ul=__ul__, xa=__xa__[i1, i2, :],
num_bd=len(__fit_wvl__), sa=__sa__, se=se, gnform=__gnform__, full=__full__, maxiter=__maxiter__, eps=__eps__)
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) 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): ...@@ -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 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. Multidimensional LUT interpolation.
:param lut1: Path radiance LUT :param lut1: LUT containing path radiance
:param lut2: Solar irradiance, spherical albedo and transmittance LUT :param lut2: LUT containing downwelling direct and diffuse surface solar irradiance, spherical albedo and
:param xnodes: Gridpoints for each LUT dimension transmittance
:param nm_nodes: Overall number of LUT gridpoints (2**ndim) :param xnodes: gridpoints for each LUT dimension
:param ndim: Dimension of LUT :param nm_nodes: overall number of LUT gridpoints (2**ndim)
:param x_cell: Interpolation grid (nm_nodes x n_dim) :param ndim: dimension of LUT
:param vtest: Atmospheric state vector (interpolation point); contains VZA, SZA, HSF, AOT, RAA, CWV :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 :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) 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): ...@@ -302,19 +304,21 @@ def interpol_lut_c(lut1, lut2, xnodes, nm_nodes, ndim, x_cell, vtest, intp_wvl):
return f_int 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. LUT interpolation based in a reduced grid.
:param lut1: Path radiance LUT :param lut1: LUT containing path radiance
:param lut2: Solar irradiance, spherical albedo and transmittance LUT :param lut2: LUT containing downwelling direct and diffuse surface solar irradiance, spherical albedo and
:param xnodes: Gridpoints for each LUT dimension transmittance
:param nm_nodes: Overall number of LUT gridpoints (2**ndim) :param xnodes: gridpoints for each LUT dimension
:param ndim: Dimension of LUT :param nm_nodes: overall number of LUT gridpoints (2**ndim)
:param x_cell: Interpolation grid (nm_nodes x n_dim) :param ndim: dimension of LUT
:param vtest: Atmospheric state vector (interpolation point); only surface elevation and CWV :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 :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) 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): ...@@ -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. 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 lut1: LUT containing path radiance
:param lut2: Solar irradiance, spherical albedo and transmittance LUT :param lut2: LUT containing downwelling direct and diffuse surface solar irradiance, spherical albedo and
:param xnodes: Gridpoints for each LUT dimension transmittance
:param nm_nodes: Overall number of LUT gridpoints (2**ndim) :param xnodes: gridpoints for each LUT dimension
:param ndim: Dimension of LUT :param nm_nodes: overall number of LUT gridpoints (2**ndim)
:param x_cell: Interpolation grid (nm_nodes x n_dim) :param ndim: dimension of LUT
:param gp: State vector of fixed scene parameters (VZA, SZA, AOT, RAA) :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 :param intp_wvl: wavelengths for which interpolation should be conducted
:return: 2D LUT including path radiance, solar irradiance, spherical albedo and transmittance interpolated :return: 2D LUT including path radiance, downwelling direct and diffuse surface solar irradiance, spherical
at gp for varying surface elevation and CWV; reduced gridpoints for each LUT dimension; reduced albedo and transmittance interpolated at gp for varying surface elevation and CWV; reduced
overall number of LUT gridpoints; reduced dimension of LUT; reduced interpolation grid 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)) 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): ...@@ -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]): for jj, cwv in enumerate(xnodes[:, 5]):
vtest = np.array([gp[0], gp[1], hsf, gp[2], gp[3], cwv]) 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, 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) vtest=vtest, intp_wvl=intp_wvl)
red_lut[ii, jj, :, 0] = pp[0, :] red_lut[ii, jj, :, 0] = pp[0, :]
red_lut[ii, jj, :, 1] = pp[1, :] red_lut[ii, jj, :, 1] = pp[1, :]
......
...@@ -25,23 +25,21 @@ ...@@ -25,23 +25,21 @@
import numpy as np import numpy as np
import dill
from tqdm import tqdm 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.conversion import generate_filter
from sicor.Tools.EnMAP.metadata import varsol from sicor.Tools.EnMAP.metadata import varsol
def wv_band_ratio(data, water_msk, land_only, 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): disable=False):
"""Band ratio water vapor retrieval. """Band ratio water vapor retrieval.
:param data: image dataset :param data: image dataset
:param water_msk: water mask :param water_msk: water mask
:param land_only: if True, CWV first guess is calculated for land surfaces only;
if False, all image pixels (land + water) are processed
:param fn_table: path to radiative transfer LUT </