Skip to content
Snippets Groups Projects
Commit 4769e9c8 authored by Niklas Bohn's avatar Niklas Bohn
Browse files

Changed some instrument specifications.

parent 530b8740
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
sicor/AC/RtFo_3_phases.py
+ 35
− 74
View file @ 4769e9c8
......@@ -69,8 +69,6 @@ class FoGen(object):
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("Getting observation metadata...")
......@@ -88,19 +86,19 @@ class FoGen(object):
self.hsf = dem
self.saa = np.full(self.data.shape[:2], options["metadata"]["saa"])
# for ii in range(self.saa.shape[0]):
# for jj in range(self.saa.shape[1]):
# if np.full(self.data.shape[:2], options["metadata"]["saa"])[ii, jj] < 0:
# self.saa[ii, jj] = 360 + np.full(self.data.shape[:2],
# options["metadata"]["saa"])[ii, jj]
# else:
# self.saa[ii, jj] = np.full(self.data.shape[:2],
# options["metadata"]["saa"])[ii, jj]
for ii in range(self.saa.shape[0]):
for jj in range(self.saa.shape[1]):
if np.full(self.data.shape[:2], options["metadata"]["saa"])[ii, jj] < 0:
self.saa[ii, jj] = 360 + np.full(self.data.shape[:2],
options["metadata"]["saa"])[ii, jj]
else:
self.saa[ii, jj] = np.full(self.data.shape[:2],
options["metadata"]["saa"])[ii, jj]
self.raa = np.abs(np.full(self.data.shape[:2], 360 - np.abs((options["metadata"]["vza"] - self.saa))))
# for ii in range(self.raa.shape[0]):
# for jj in range(self.raa.shape[1]):
# if self.raa[ii, jj] > 180:
# self.raa[ii, jj] = 360 - self.raa[ii, jj]
for ii in range(self.raa.shape[0]):
for jj in range(self.raa.shape[1]):
if self.raa[ii, jj] > 180:
self.raa[ii, jj] = 360 - self.raa[ii, jj]
# check if observation metadata values are within LUT value ranges
self.logger.info("Checking if observation metadata values are within LUT value ranges...")
......@@ -155,11 +153,6 @@ class FoGen(object):
self.wvl = np.array(options["sensor"]["fit"]["wvl_center"])
self.fwhm = np.array(options["sensor"]["fit"]["fwhm"])
# vnir wvl
if self.instrument == "EnMAP":
self.wvl_vnir = self.wvl[:88]
self.fwhm_vnir = self.fwhm[:88]
# fit wvl
self.fit_wvl = np.array(options["sensor"]["fit"]["idx"])
self.wvl_sel = self.wvl[self.fit_wvl]
......@@ -184,24 +177,21 @@ class FoGen(object):
self.s1 = float(solar_1) * (10e6 / self.wvl_sel[-1]) ** 2
# 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:
# self.path_sicorlib = os.path.dirname(pkgutil.get_loader("sicor").path)
# if os.path.isfile(os.path.join(self.path_sicorlib, "AC", "data", "EnMAP_LUT_MOD5_formatted_1nm")):
# self.fn_table = os.path.join(self.path_sicorlib, "AC", "data", "EnMAP_LUT_MOD5_formatted_1nm")
# else:
# raise FileNotFoundError("LUT file was not found. Make sure to indicate path in options file or to "
# "store the LUT at /sicor/AC/data/ directory, otherwise, the AC will not work.")
# luts, axes_x, axes_y, wvl, lut1, lut2, xnodes, nm_nodes, ndim, x_cell = read_lut_enmap_formatted(
# file_lut=self.fn_table)
self.logger.info("Loading RT LUT...")
if os.path.isfile(options["retrieval"]["fn_LUT"]):
self.fn_table = options["retrieval"]["fn_LUT"]
else:
self.path_sicorlib = os.path.dirname(pkgutil.get_loader("sicor").path)
if os.path.isfile(os.path.join(self.path_sicorlib, "AC", "data", "EnMAP_LUT_MOD5_formatted_1nm")):
self.fn_table = os.path.join(self.path_sicorlib, "AC", "data", "EnMAP_LUT_MOD5_formatted_1nm")
else:
raise FileNotFoundError("LUT file was not found. Make sure to indicate path in options file or to "
"store the LUT at /sicor/AC/data/ directory, otherwise, the AC will not work.")
with open("/Users/bohn/Desktop/hyspex_mjolnir_MOD5_LUT.dill", "rb") as fl:
wvl, xnodes, nm_nodes, ndim, x_cell, lut1, lut2 = dill.load(fl)
luts, axes_x, axes_y, wvl, lut1, lut2, xnodes, nm_nodes, ndim, x_cell = read_lut_enmap_formatted(
file_lut=self.fn_table)
self.wvl_lut = wvl * 1000
self.wvl_lut = wvl
# resample LUT to instrument wavelengths
self.logger.info("Resampling LUT to instrument wavelengths...")
......@@ -214,8 +204,8 @@ class FoGen(object):
self.ndim = ndim
self.x_cell = x_cell
lut2_all_res_fit = np.zeros((2, 2, 2, 2, 1, 7, len(self.wvl_sel), 4))
lut2_all_res_full = np.zeros((2, 2, 2, 2, 1, 7, len(self.wvl), 4))
lut2_all_res_fit = np.zeros((5, 6, 4, 6, 1, 7, len(self.wvl_sel), 4))
lut2_all_res_full = np.zeros((5, 6, 4, 6, 1, 7, len(self.wvl), 4))
lut1_res_fit = lut1[:, :, :, :, :, :, :, 0] @ self.s_norm_fit
lut1_res_full = lut1[:, :, :, :, :, :, :, 0] @ self.s_norm_full
......@@ -929,43 +919,14 @@ def __minimize__(fo, opt_func, unknowns=False, logger=None):
cwv_fg = np.full(fo.data.shape[:2], fo.prior_mean[0])
else:
logger.info("Performing first guess water vapor retrieval based on a common band ratio...")
cwv_fg_vnir = wv_band_ratio(data=fo.data_fg, fn_table=fo.fn_table, vza=fo.pt[:, :, 0].mean(),
cwv_fg = wv_band_ratio(data=fo.data, fn_table=fo.fn_table, vza=fo.pt[:, :, 0].mean(),
sza=fo.pt[:, :, 1].mean(), dem=fo.pt[:, :, 2].mean(), aot=fo.pt[:, :, 3].mean(),
raa=fo.pt[:, :, 4].mean(), intp_wvl=fo.wvl_vnir, intp_fwhm=fo.fwhm_vnir,
jday=fo.jday, month=fo.month, idx=[73, 76, 81], disable=fo.disable_progressbars)
cwv_fg_vnir[np.isnan(cwv_fg_vnir)] = fo.prior_mean[0]
# transform water vapor first guess to SWIR sensor geometry
logger.info("Transforming CWV first guess to SWIR sensor geometry...")
cwv_fg = fo.VS_SGT.transform_sensorgeo_VNIR_to_SWIR(cwv_fg_vnir)
raa=fo.pt[:, :, 4].mean(), intp_wvl=fo.wvl, intp_fwhm=fo.fwhm,
jday=fo.jday, month=fo.month, idx=[942, 991, 1072], disable=fo.disable_progressbars)
cwv_fg[np.isnan(cwv_fg)] = fo.prior_mean[0]
# perform first guess liquid water retrieval based on the NDWI
if fo.use_prior_mean[3]:
liq_fg = np.full(fo.data.shape[:2], fo.prior_mean[3])
else:
logger.info("Performing first guess liquid water retrieval based on the NDWI...")
ndwi = np.zeros(fo.data.shape[:2])
for ii in tqdm(range(fo.data.shape[0]), disable=fo.disable_progressbars):
for jj in range(fo.data.shape[1]):
ndwi[ii, jj] = (fo.data_fg_trans[ii, jj, 72] - fo.data[ii, jj, 94]) / \
(fo.data_fg_trans[ii, jj, 72] + fo.data[ii, jj, 94])
liq_fg = 0.49585423 * ndwi**2 - 0.0396154 * ndwi - 0.0174602
# perform first guess snow retrieval based on the NDSI
if fo.use_prior_mean[4]:
ice_fg = np.full(fo.data.shape[:2], fo.prior_mean[4])
else:
logger.info("Performing first guess snow retrieval based on the NDSI...")
ndsi = np.zeros(fo.data.shape[:2])
ice_fg = np.zeros(fo.data.shape[:2])
for ii in tqdm(range(fo.data.shape[0]), disable=fo.disable_progressbars):
for jj in range(fo.data.shape[1]):
ndsi[ii, jj] = (fo.data_fg_trans[ii, jj, 29] - fo.data[ii, jj, 57]) / \
(fo.data_fg_trans[ii, jj, 29] + fo.data[ii, jj, 57])
if ndsi[ii, jj] > 0.9:
ice_fg[ii, jj] = 0.25
else:
ice_fg[ii, jj] = 0
# perform calculation of first guess for intercept and slope of absorption feature continuum
logger.info("Calculating first guess for intercept and slope of absorption feature continuum...")
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment