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

Implemented dimensionality reduction of LUT grid to increase interpolation speed.


Signed-off-by: Niklas Bohn's avatarnbohn <nbohn@gfz-potsdam.de>
parent d66bcbfb
......@@ -41,7 +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, interpol_lut_c, download_LUT
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.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
......@@ -681,6 +682,23 @@ class Fo(object):
self.lut2_vnir = lut2_all_res_vnir
self.lut2_swir = lut2_all_res_swir
# reduce grid dimensionality of LUTs to increase interpolation speed
self.logger.info("Reducing grid dimensionality of LUT to increase interpolation speed...")
self.gp = [enmap_l1b.meta.geom_view_zenith, enmap_l1b.meta.geom_sun_zenith,
options["retrieval"]["default_aot_value"], self.raa.mean()]
self.lut_red_fit, self.xnodes_red, self.nm_nodes_red, self.ndim_red, self.x_cell_red = reduce_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, gp=self.gp, intp_wvl=self.wvl_sel)
self.lut_red_vnir, _, _, _, _ = reduce_lut(lut1=self.lut1_vnir, lut2=self.lut2_vnir, xnodes=self.xnodes,
nm_nodes=self.nm_nodes, ndim=self.ndim, x_cell=self.x_cell,
gp=self.gp, intp_wvl=self.wvl_vnir)
self.lut_red_swir, _, _, _, _ = reduce_lut(lut1=self.lut1_swir, lut2=self.lut2_swir, xnodes=self.xnodes,
nm_nodes=self.nm_nodes, ndim=self.ndim, x_cell=self.x_cell,
gp=self.gp, intp_wvl=self.wvl_swir)
# calculate absorption coefficients of liquid water and ice
warnings.filterwarnings("ignore")
self.logger.info("Calculating absorption coefficients of liquid water and ice...")
......@@ -864,15 +882,11 @@ class Fo(object):
:return: modeled TOA radiance
"""
# 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)
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_l0 = f_int[0, :] * 1.e+3
f_int_edir = f_int[1, :] * 1.e+3
......@@ -904,35 +918,20 @@ class Fo(object):
:return: modeled surface reflectance
"""
# LUT interpolation
vtest = np.asarray([pt[0], pt[1], pt[2], pt[3], pt[4], xx])
vtest = np.asarray([pt[2], xx])
if mode == "vnir":
f_int = interpol_lut_c(lut1=self.lut1_vnir,
lut2=self.lut2_vnir,
xnodes=self.xnodes,
nm_nodes=self.nm_nodes,
ndim=self.ndim,
x_cell=self.x_cell,
vtest=vtest,
intp_wvl=self.wvl_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)
elif mode == "swir":
f_int = interpol_lut_c(lut1=self.lut1_swir,
lut2=self.lut2_swir,
xnodes=self.xnodes,
nm_nodes=self.nm_nodes,
ndim=self.ndim,
x_cell=self.x_cell,
vtest=vtest,
intp_wvl=self.wvl_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)
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_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_l0 = f_int[0, :] * 1.e+3 * self.fac
f_int_edir = f_int[1, :] * 1.e+3 * self.fac
......
......@@ -32,12 +32,12 @@ from os.path import isfile, join, dirname
import inspect
import sys
import numpy as np
from numba import jit
import urllib.request
def get_data_file(module_name, file_basename):
"""Load data file.
"""
Load data file.
:param module_name: name of python module where data file is located
:param file_basename: name of data file
......@@ -80,12 +80,12 @@ def download_LUT(path_LUT_default):
def read_lut_enmap_formatted(file_lut):
"""Read MODTRAN LUT.
"""
Read MODTRAN LUT.
:param file_lut: path to LUT
:return: LUT of atmospheric functions, x and y axes grid points, LUT wavelengths
"""
ndim = 6
with open(file_lut, mode='rb') as fd:
......@@ -238,9 +238,9 @@ def read_lut_enmap_formatted(file_lut):
return luts, axes_x, axes_y, wvl, lut1, lut2, xnodes, nm_nodes, ndim, x_cell
@jit(nopython=True)
def interpol_lut_c(lut1, lut2, xnodes, nm_nodes, ndim, x_cell, vtest, intp_wvl):
"""Multidimensional LUT interpolation.
"""
Multidimensional LUT interpolation based on numba.
:param lut1: Path radiance LUT
:param lut2: Solar irradiance, spherical albedo and transmittance LUT
......@@ -248,14 +248,12 @@ def interpol_lut_c(lut1, lut2, xnodes, nm_nodes, ndim, x_cell, vtest, intp_wvl):
: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)
: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
"""
lim = np.zeros((2, ndim), np.int8)
# TODO compute that on the full cube
for ii in range(ndim):
if vtest[ii] >= xnodes[:, ii].max():
vtest[ii] = 0.99 * xnodes[-1, ii]
......@@ -297,9 +295,105 @@ def interpol_lut_c(lut1, lut2, xnodes, nm_nodes, ndim, x_cell, vtest, intp_wvl):
for i in range(ndim):
vtest[i] = (vtest[i] - xnodes[lim[0, i], i]) / (xnodes[lim[1, i], i] - xnodes[lim[0, i], i])
# weights = np.abs(np.prod(vtest - x_cell[::-1, :], axis=1)).reshape(1, nm_nodes, 1) # unsupported by numba
diffs = vtest - x_cell[::-1, :]
weights = np.abs(np.array([np.prod(diffs[i]) for i in range(nm_nodes)])).reshape(1, nm_nodes, 1)
f_int = np.sum(weights * lut_cell, axis=1)
return f_int
def interpol_lut_c_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
:param intp_wvl: wavelengths for which interpolation should be conducted
:return: Path radiance, solar irradiance, spherical albedo and transmittance interpolated at vtest
"""
lim = np.zeros((2, ndim), np.int8)
for ii in range(ndim):
if vtest[ii] >= xnodes[:, ii].max():
vtest[ii] = 0.99 * xnodes[-1, ii]
wh = np.where(vtest[ii] < xnodes[:, ii])[0]
lim[0, ii] = wh[0] - 1
lim[1, ii] = wh[0]
lut_cell = np.zeros((5, nm_nodes, len(intp_wvl)))
cont = 0
for i in range(2):
iv = lim[i, 0]
for j in range(2):
jv = lim[j, 1]
lut_cell[0, cont, :] = \
lut1[iv, jv, :]
for ind in range(1, 5):
lut_cell[ind, cont, :] = \
lut2[iv, jv, :, ind - 1]
cont += 1
for i in range(ndim):
vtest[i] = (vtest[i] - xnodes[lim[0, i], i]) / (xnodes[lim[1, i], i] - xnodes[lim[0, i], i])
diffs = vtest - x_cell[::-1, :]
weights = np.abs(np.array([np.prod(diffs[i]) for i in range(nm_nodes)])).reshape(1, nm_nodes, 1)
f_int = np.sum(weights * lut_cell, axis=1)
return f_int
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 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
"""
red_lut = np.zeros((len(xnodes[:4, 2]), len(xnodes[:, 5]), len(intp_wvl), lut2.shape[-1] + 1))
for ii, hsf in enumerate(xnodes[:4, 2]):
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)
red_lut[ii, jj, :, 0] = pp[0, :]
red_lut[ii, jj, :, 1] = pp[1, :]
red_lut[ii, jj, :, 2] = pp[2, :]
red_lut[ii, jj, :, 3] = pp[3, :]
red_lut[ii, jj, :, 4] = pp[4, :]
new_xnodes = np.zeros((xnodes.shape[0], 2))
new_xnodes[:, 0] = xnodes[:, 2]
new_xnodes[:, 1] = xnodes[:, 5]
new_nm_nodes = new_xnodes.shape[1] ** 2
new_ndim = new_xnodes.shape[1]
new_x_cell = np.zeros((new_nm_nodes, new_ndim))
for ii in range(new_nm_nodes):
new_x_cell[ii, :] = x_cell[ii][-2:]
return red_lut, new_xnodes, new_nm_nodes, new_ndim, new_x_cell
......@@ -85,7 +85,8 @@ def wv_band_ratio(data, fn_table, vza, sza, dem, aot, raa, intp_wvl, intp_fwhm,
for jj, cwv in enumerate(cwvs):
vtest = np.asarray([vza, sza, hsf, aot, raa, cwv])
f_int = interpol_lut_c(lut1_res, lut2_res, xnodes, nm_nodes, ndim, x_cell, vtest, intp_wvl)
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_l0 = f_int[0, :] * 1.e+3
f_int_edir = f_int[1, :] * 1.e+3
......
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