Commit 48411995 authored by Niklas Bohn's avatar Niklas Bohn
Browse files

Added "if" statement to run water vapor first guess only for land pixels.

parent 4e8ad032
......@@ -9,6 +9,7 @@ History
New features:
* Dimensionality reduction of LUT grid to increase interpolation speed.
* Updated final log message of SICOR AC for EnMAP.
* First guess water vapor retrieval is only applied to land pixels if land_only is set to true.
Bugfixes:
* Fixed bug in empirical line function which produced unrealistic peaks in water reflectance spectra.
......
......@@ -482,10 +482,10 @@ class Fo(object):
self.data_swir = enmap_l1b.swir.data
self.resamp_alg = options["sensor"]["resamp_alg"]
self.land_only = options["retrieval"]["land_only"]
self.water_mask_vnir = enmap_l1b.vnir.mask_landwater[:, :]
if self.land_only:
logger.info("SICOR is applied to land pixels only. This may result in edge effects, e.g., at coastlines...")
self.water_mask_vnir = enmap_l1b.vnir.mask_landwater[:, :]
self.water_mask_swir = enmap_l1b.transform_vnir_to_swir_raster(array_vnirsensorgeo=self.water_mask_vnir,
resamp_alg=self.resamp_alg,
respect_keystone=False)
......@@ -731,6 +731,8 @@ class Fo(object):
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,
land_only=self.land_only,
fn_table=self.fn_table,
vza=self.pt[:, :, 0].mean(),
sza=self.pt[:, :, 1].mean(),
......
......@@ -33,10 +33,14 @@ from sicor.Tools.EnMAP.conversion import generate_filter
from sicor.Tools.EnMAP.metadata import varsol
def wv_band_ratio(data, fn_table, vza, sza, dem, aot, raa, intp_wvl, intp_fwhm, jday, month, idx, disable=False):
def wv_band_ratio(data, water_msk, land_only, fn_table, 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 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
:param vza: viewing zenitg angle
:param sza: sun zenith angle
......@@ -60,6 +64,8 @@ def wv_band_ratio(data, fn_table, vza, sza, dem, aot, raa, intp_wvl, intp_fwhm,
toa_sub[:, 1] = np.ndarray.flatten(data[:, :, idx[2]])
cnt_land = len(toa_sub[:, 0])
water_msk_flat = np.ndarray.flatten(water_msk)
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
......@@ -144,32 +150,35 @@ def wv_band_ratio(data, fn_table, vza, sza, dem, aot, raa, intp_wvl, intp_fwhm,
wv_arr[:] = np.nan
for ind in tqdm(range(0, cnt_land), disable=disable):
dem_fac_pix = dem_fac[ind]
if ll_cal[ind] <= rhos[0]:
ll_cal[ind] = rhos[0] + 0.01 * rhos[0]
ll_cal_low = ll_cal[ind] >= rhos
ll_cal_high = ll_cal[ind] <= rhos
idx_low = np.where(ll_cal_low)
idx_high = np.where(ll_cal_high)
rfl_fac = (ll_cal[ind] - rhos[idx_low[0][-1]]) / (rhos[idx_high[0][0]] - rhos[idx_low[0][-1]])
rfl_fac_pix = rfl_fac
rfl_sl = np.ndarray.flatten(rfl_sl_img)[ind]
rfl_sl_pix = (rfl_sl - rfl_sl_gr[0]) / (rfl_sl_gr[1] - rfl_sl_gr[0])
cf_int = (1 - dem_fac_pix) * (1 - rfl_fac_pix) * (1 - rfl_sl_pix) * cf_arr[:, 0, idx_low[0][-1], 0] + (
1 - dem_fac_pix) * (1 - rfl_fac_pix) * rfl_sl_pix * cf_arr[:, 0, idx_low[0][-1], 1] + (
1 - dem_fac_pix) * rfl_fac_pix * (1 - rfl_sl_pix) * cf_arr[:, 0, idx_high[0][0], 0] + (
1 - dem_fac_pix) * rfl_fac_pix * rfl_sl_pix * cf_arr[:, 0, idx_high[0][0], 1] + dem_fac_pix * (
1 - rfl_fac_pix) * (1 - rfl_sl_pix) * cf_arr[:, 1, idx_low[0][-1], 0] + dem_fac_pix * (
1 - rfl_fac_pix) * rfl_sl_pix * cf_arr[:, 1, idx_low[0][-1], 1] + dem_fac_pix * rfl_fac_pix * (
1 - rfl_sl_pix) * cf_arr[:, 1, idx_high[0][0],
0] + dem_fac_pix * rfl_fac_pix * rfl_sl_pix * cf_arr[:, 1, idx_high[0][0], 1]
wv = cf_int[2] + alog_rat_wv_img[ind] * cf_int[1] + alog_rat_wv_img[ind] * alog_rat_wv_img[ind] * cf_int[0]
wv_arr[ind] = wv
if land_only and water_msk_flat[ind] != 1:
pass
else:
dem_fac_pix = dem_fac[ind]
if ll_cal[ind] <= rhos[0]:
ll_cal[ind] = rhos[0] + 0.01 * rhos[0]
ll_cal_low = ll_cal[ind] >= rhos
ll_cal_high = ll_cal[ind] <= rhos
idx_low = np.where(ll_cal_low)
idx_high = np.where(ll_cal_high)
rfl_fac = (ll_cal[ind] - rhos[idx_low[0][-1]]) / (rhos[idx_high[0][0]] - rhos[idx_low[0][-1]])
rfl_fac_pix = rfl_fac
rfl_sl = np.ndarray.flatten(rfl_sl_img)[ind]
rfl_sl_pix = (rfl_sl - rfl_sl_gr[0]) / (rfl_sl_gr[1] - rfl_sl_gr[0])
cf_int = (1 - dem_fac_pix) * (1 - rfl_fac_pix) * (1 - rfl_sl_pix) * cf_arr[:, 0, idx_low[0][-1], 0] + (
1 - dem_fac_pix) * (1 - rfl_fac_pix) * rfl_sl_pix * cf_arr[:, 0, idx_low[0][-1], 1] + (
1 - dem_fac_pix) * rfl_fac_pix * (1 - rfl_sl_pix) * cf_arr[:, 0, idx_high[0][0], 0] + (
1 - dem_fac_pix) * rfl_fac_pix * rfl_sl_pix * cf_arr[:, 0, idx_high[0][0], 1] + dem_fac_pix * (
1 - rfl_fac_pix) * (1 - rfl_sl_pix) * cf_arr[:, 1, idx_low[0][-1], 0] + dem_fac_pix * (
1 - rfl_fac_pix) * rfl_sl_pix * cf_arr[:, 1, idx_low[0][-1], 1] + dem_fac_pix * rfl_fac_pix * (
1 - rfl_sl_pix) * cf_arr[:, 1, idx_high[0][0],
0] + dem_fac_pix * rfl_fac_pix * rfl_sl_pix * cf_arr[:, 1, idx_high[0][0], 1]
wv = cf_int[2] + alog_rat_wv_img[ind] * cf_int[1] + alog_rat_wv_img[ind] * alog_rat_wv_img[ind] * cf_int[0]
wv_arr[ind] = wv
wv_arr_img = np.reshape(wv_arr, (data.shape[:2]))
......
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