Commit 5e9d23f4 authored by Niklas Bohn's avatar Niklas Bohn

Adjusted generic forward operator to version 0.15.0.

parent e2ea7036
......@@ -66,13 +66,14 @@ class FoGen(object):
"""
self.logger = logger or logging.getLogger(__name__)
self.cpu = options["retrieval"]["cpu"]
self.disable_progressbars = options["retrieval"]["disable_progressbars"]
self.instrument = options["sensor"]["name"]
self.data = data
if self.instrument == "EnMAP":
self.data_fg = data[:, :, :88]
# get observation metadata
self.logger.info("get observation metadata...")
self.logger.info("Getting observation metadata...")
self.jday = options["metadata"]["jday"]
self.month = options["metadata"]["month"]
self.dsol = varsol(jday=self.jday, month=self.month)
......@@ -103,7 +104,7 @@ class FoGen(object):
self.raa[ii, jj] = 360 - self.raa[ii, jj]
# check if observation metadata values are within LUT value ranges
self.logger.info("check if observation metadata values are within LUT value ranges...")
self.logger.info("Checking if observation metadata values are within LUT value ranges...")
try:
assert np.min(self.vza) >= 0 and np.max(self.vza) <= 40
except AssertionError:
......@@ -131,8 +132,7 @@ class FoGen(object):
self.pt[:, :, 3] = np.full(self.data.shape[:2], options["metadata"]["aot"])
self.pt[:, :, 4] = self.raa
# self.par_opt = ["cwv", "cwc", "ice", "lai", "dasf", "p_max", "k", "b"]
self.par_opt = ["cwv", "ewt", "ice", "dasf", "p"]
self.par_opt = ["water vapor", "intercept", "slope", "liquid water", "ice"]
self.state_dim = len(self.par_opt)
self.prior_mean, self.use_prior_mean, self.ll, self.ul, self.prior_sigma, self.unknowns = [], [], [], [], [], []
......@@ -168,7 +168,7 @@ class FoGen(object):
self.snr_fit = np.array(options["sensor"]["fit"]["snr"])
# get solar irradiances for absorption feature shoulders wavelengths
self.logger.info("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_table(new_kur_fn, delim_whitespace=True)
freq_0 = 10e6 / self.wvl_sel[0]
......@@ -185,7 +185,7 @@ class FoGen(object):
self.s1 = float(solar_1) * (10e6 / self.wvl_sel[-1]) ** 2
# load RT LUT
self.logger.info("load RT LUT...")
self.logger.info("Loading RT LUT...")
if os.path.isfile(options["retrieval"]["fn_LUT"]):
self.fn_table = options["retrieval"]["fn_LUT"]
else:
......@@ -201,7 +201,7 @@ class FoGen(object):
self.wvl_lut = wvl
# resample LUT to instrument wavelengths
self.logger.info("resample LUT to instrument wavelengths...")
self.logger.info("Resampling LUT to instrument wavelengths...")
self.s_norm_fit = generate_filter(wvl_m=self.wvl_lut, wvl=self.wvl_sel, wl_resol=self.fwhm_sel)
self.s_norm_full = generate_filter(wvl_m=self.wvl_lut, wvl=self.wvl, wl_resol=self.fwhm)
......@@ -226,7 +226,8 @@ class FoGen(object):
self.lut2_full = lut2_all_res_full
# calculate absorption coefficients of liquid water and ice
self.logger.info("calculate absorption coefficients of liquid water and ice...")
self.logger.info("Calculating absorption coefficients of liquid water and ice...")
warnings.filterwarnings("ignore")
path_k = get_data_file(module_name="sicor", file_basename="k_liquid_water_ice.xlsx")
k_wi = pd.read_excel(io=path_k)
......@@ -242,32 +243,22 @@ class FoGen(object):
self.abs_co_i = 4 * np.pi * self.ki / self.wvl_sel
# load mean solar exoatmospheric irradiances
self.logger.info("load 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
def surface_model(self, xx):
"""Nonlinear surface reflectance model using leaf and snow single scattering albedo as well as photon
re-collision probability for the retrieval of liquid water and ice path lengths. This model accounts for
absorption caused by biogeochemical leaf and snow constituents as well as for multiple scattering determined
by structural and geometric parameters.
"""Nonlinear surface reflectance model using the Beer-Lambert attenuation law for the retrieval of liquid water
and ice path lengths.
:param xx: state vector, must be in the order [cwv, cwc, ice, lai, dasf, p_max, k, b]
:param xx: state vector, must be in the order [vapor, intercept, slope, liquid, ice]
:return: modeled surface reflectance
"""
# combined leaf and snow single scattering albedo (Marshak & Knyazikhin 2017)
w0 = np.exp(-xx[1] * 1e7 * self.abs_co_w - xx[2] * 1e7 * self.abs_co_i)
# canopy photon re-collision probability (Smolander & Stenberg 2005)
# p = xx[5] * (1 - np.exp(-xx[6] * xx[3] ** xx[7]))
# canopy single scattering albedo (Lewis & Disney 2007)
W = w0 * (1 - xx[4]) / (1 - xx[4] * w0)
# canopy bidirectional reflectance factor (Knyazikhin et al. 2013)
rho = xx[3] * W
# modeling of surface reflectance
attenuation = np.exp(-xx[3] * 1e7 * self.abs_co_w - xx[4] * 1e7 * self.abs_co_i)
rho = (xx[1] + xx[2] * self.wvl_sel) * attenuation
return rho
......@@ -349,15 +340,15 @@ class FoGen(object):
kb_rt = []
eps = 1e-5
perturb = (1.0 + eps)
unknowns = ["Skyview", "H2O_ABSCO"]
unknowns = ["skyview", "water_vapor_absorption_coefficients"]
for unknown in unknowns:
if unknown == "Skyview":
if unknown == "skyview":
# perturb the sky view
rdne = self.toa_rad(xx=xx, pt=pt, perturb=perturb)
dx = (rdne - rdn) / eps
kb_rt.append(dx)
elif unknown == "H2O_ABSCO":
elif unknown == "water_vapor_absorption_coefficients":
xx_perturb = xx.copy()
xx_perturb[0] = xx[0] * perturb
rdne = self.toa_rad(xx=xx_perturb, pt=pt)
......@@ -379,7 +370,7 @@ class FoGen(object):
kb_surface = []
eps = 1e-5
perturb = (1.0 * eps)
unknowns = ["Liquid_ABSCO", "Ice_ABSCO"]
unknowns = ["liquid_water_absorption_coefficients", "ice_absorption_coefficients"]
for ii, unknown in enumerate(unknowns):
xx_perturb = xx.copy()
......
......@@ -24,41 +24,42 @@
"retrieval": {
"fn_LUT": "",
"cpu":12,
"disable_progressbars": false,
"state_vector": {
"cwv": {
"water_vapor": {
"prior_mean": 2.5,
"use_prior_mean": false,
"ll": 0.0001,
"ul": 4.9999,
"prior_sigma": 0.5
},
"ewt": {
"prior_mean": 0.0111,
"use_prior_mean": true,
"ll": 0.0001,
"ul": 0.0999,
"prior_sigma": 0.006979349
"ll": 0.001,
"ul": 4.999,
"prior_sigma": 1e3
},
"ice": {
"prior_mean": 0,
"use_prior_mean": true,
"ll": 0.0001,
"ul": 0.4999,
"prior_sigma": 1e-100
"intercept": {
"prior_mean": 0.3,
"use_prior_mean": false,
"ll": 0.001,
"ul": 0.999,
"prior_sigma": 1e3
},
"slope": {
"prior_mean": 0.0002,
"use_prior_mean": false,
"ll": -0.0004,
"ul": 0.0004,
"prior_sigma": 1e3
},
"dasf": {
"prior_mean": 0.15,
"liquid_water": {
"prior_mean": 0.02,
"use_prior_mean": true,
"ll": 0.0001,
"ul": 0.9999,
"prior_sigma": 0.000001
"ll": 0.001,
"ul": 0.499,
"prior_sigma": 1e3
},
"p": {
"prior_mean": 0.7954979526021773,
"ice": {
"prior_mean": 0.02,
"use_prior_mean": true,
"ll": 0.4001,
"ul": 0.9999,
"prior_sigma": 0.13466832630044967
"ll": 0.001,
"ul": 0.499,
"prior_sigma": 1e3
}
},
"unknowns": {
......
......@@ -27,6 +27,7 @@
import logging
from tqdm import tqdm
import numpy as np
import warnings
from sicor.AC.RtFo_3_phases import FoGen, FoFunc, __minimize__
......@@ -50,21 +51,22 @@ def make_ac_generic(data_l1b, fo, xx, logger=None):
data_l2a = np.full(data_l1b.shape, np.NaN, dtype=np.float)
ncols, nrows = (data_l1b.shape[1], data_l1b.shape[0])
logger.info("Perform columnwise ac...")
logger.info("Performing columnwise ac...")
warnings.filterwarnings("ignore")
for icol in tqdm(range(ncols)):
for irow in range(nrows):
data_l2a[irow, icol, :] = fo.surf_ref(dt=data_l1b[irow, icol, :], xx=[xx["cwv_model"][irow, icol]],
pt=fo.pt[irow, icol, :], mode="full")
swir_feature = fo.surface_model(xx=[xx["cwv_model"][irow, icol], xx["ewt_model"][irow, icol],
xx["ice_model"][irow, icol], xx["dasf_model"][irow, icol],
xx["p_model"][irow, icol]])
swir_feature = fo.surface_model(xx=[xx["cwv_model"][irow, icol], xx["intercept_model"][irow, icol],
xx["slope_model"][irow, icol], xx["liq_model"][irow, icol],
xx["ice_model"][irow, icol]])
data_l2a[irow, icol, fo.fit_wvl] = swir_feature
return data_l2a
def sicor_ac_generic(data_l1b, options, dem=None, unknowns=False, logger=None):
"""Atmospheric correction for sensor-independent hyperspectral L1 products, including a three phases of water
"""Atmospheric correction for sensor-independent imaging spectroscopy L1 products, including a three phases of water
retrieval.
:arg data_l1b:
......@@ -83,14 +85,14 @@ def sicor_ac_generic(data_l1b, options, dem=None, unknowns=False, logger=None):
"""
logger = logger or logging.getLogger(__name__)
logger.info("set up forward operator...")
logger.info("Setting up forward operator...")
fo_gen = FoGen(data=data_l1b, options=options, dem=dem, logger=logger)
fo_func_gen = FoFunc(fo=fo_gen)
logger.info("start 3 phases of water retrieval...")
logger.info("Starting 3 phases of water retrieval...")
res = __minimize__(fo=fo_gen, opt_func=fo_func_gen, unknowns=unknowns, logger=logger)
logger.info("start surface reflectance retrieval...")
logger.info("Starting surface reflectance retrieval...")
data_l2a = make_ac_generic(data_l1b=data_l1b, fo=fo_gen, xx=res, logger=logger)
logger.info("%s atmospheric correction successfully finished!" % options["sensor"]["name"])
......
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