Commit 8fcb6516 authored by Niklas Bohn's avatar Niklas Bohn
Browse files

Merge branch 'enhancement/Issue#88' into 'master'

Implemented dimensionality reduction of LUT grid to increase interpolation speed.

See merge request !90
parents d66bcbfb f89f484e
...@@ -41,7 +41,8 @@ import platform ...@@ -41,7 +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, 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.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
...@@ -681,6 +682,23 @@ class Fo(object): ...@@ -681,6 +682,23 @@ class Fo(object):
self.lut2_vnir = lut2_all_res_vnir self.lut2_vnir = lut2_all_res_vnir
self.lut2_swir = lut2_all_res_swir 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 # calculate absorption coefficients of liquid water and ice
warnings.filterwarnings("ignore") warnings.filterwarnings("ignore")
self.logger.info("Calculating absorption coefficients of liquid water and ice...") self.logger.info("Calculating absorption coefficients of liquid water and ice...")
...@@ -864,15 +882,11 @@ class Fo(object): ...@@ -864,15 +882,11 @@ class Fo(object):
:return: modeled TOA radiance :return: modeled TOA radiance
""" """
# LUT interpolation # LUT interpolation
vtest = np.asarray([pt[0], pt[1], pt[2], pt[3], pt[4], xx[0]]) vtest = np.asarray([pt[2], xx[0]])
f_int = interpol_lut_c(lut1=self.lut1_fit,
lut2=self.lut2_fit, f_int = interpol_lut_c_red(lut1=self.lut_red_fit[:, :, :, 0], lut2=self.lut_red_fit[:, :, :, 1:],
xnodes=self.xnodes, xnodes=self.xnodes_red, nm_nodes=self.nm_nodes_red, ndim=self.ndim_red,
nm_nodes=self.nm_nodes, x_cell=self.x_cell_red, vtest=vtest, intp_wvl=self.wvl_sel)
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
...@@ -904,35 +918,20 @@ class Fo(object): ...@@ -904,35 +918,20 @@ class Fo(object):
:return: modeled surface reflectance :return: modeled surface reflectance
""" """
# LUT interpolation # 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": if mode == "vnir":
f_int = interpol_lut_c(lut1=self.lut1_vnir, f_int = interpol_lut_c_red(lut1=self.lut_red_vnir[:, :, :, 0], lut2=self.lut_red_vnir[:, :, :, 1:],
lut2=self.lut2_vnir, xnodes=self.xnodes_red, nm_nodes=self.nm_nodes_red, ndim=self.ndim_red,
xnodes=self.xnodes, x_cell=self.x_cell_red, vtest=vtest, intp_wvl=self.wvl_vnir)
nm_nodes=self.nm_nodes,
ndim=self.ndim,
x_cell=self.x_cell,
vtest=vtest,
intp_wvl=self.wvl_vnir)
elif mode == "swir": elif mode == "swir":
f_int = interpol_lut_c(lut1=self.lut1_swir, f_int = interpol_lut_c_red(lut1=self.lut_red_swir[:, :, :, 0], lut2=self.lut_red_swir[:, :, :, 1:],
lut2=self.lut2_swir, xnodes=self.xnodes_red, nm_nodes=self.nm_nodes_red, ndim=self.ndim_red,
xnodes=self.xnodes, x_cell=self.x_cell_red, vtest=vtest, intp_wvl=self.wvl_swir)
nm_nodes=self.nm_nodes,
ndim=self.ndim,
x_cell=self.x_cell,
vtest=vtest,
intp_wvl=self.wvl_swir)
else: else:
f_int = interpol_lut_c(lut1=self.lut1_fit, f_int = interpol_lut_c_red(lut1=self.lut_red_fit[:, :, :, 0], lut2=self.lut_red_fit[:, :, :, 1:],
lut2=self.lut2_fit, xnodes=self.xnodes_red, nm_nodes=self.nm_nodes_red, ndim=self.ndim_red,
xnodes=self.xnodes, x_cell=self.x_cell_red, vtest=vtest, intp_wvl=self.wvl_sel)
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
......
...@@ -32,12 +32,12 @@ from os.path import isfile, join, dirname ...@@ -32,12 +32,12 @@ from os.path import isfile, join, dirname
import inspect import inspect
import sys import sys
import numpy as np import numpy as np
from numba import jit
import urllib.request import urllib.request
def get_data_file(module_name, file_basename): 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 module_name: name of python module where data file is located
:param file_basename: name of data file :param file_basename: name of data file
...@@ -80,12 +80,12 @@ def download_LUT(path_LUT_default): ...@@ -80,12 +80,12 @@ def download_LUT(path_LUT_default):
def read_lut_enmap_formatted(file_lut): def read_lut_enmap_formatted(file_lut):
"""Read MODTRAN LUT. """
Read MODTRAN LUT.
:param file_lut: path to LUT :param file_lut: path to LUT
:return: LUT of atmospheric functions, x and y axes grid points, LUT wavelengths :return: LUT of atmospheric functions, x and y axes grid points, LUT wavelengths
""" """
ndim = 6 ndim = 6
with open(file_lut, mode='rb') as fd: with open(file_lut, mode='rb') as fd:
...@@ -238,9 +238,9 @@ def read_lut_enmap_formatted(file_lut): ...@@ -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 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): 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 lut1: Path radiance LUT
:param lut2: Solar irradiance, spherical albedo and transmittance 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): ...@@ -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 nm_nodes: Overall number of LUT gridpoints (2**ndim)
:param ndim: Dimension of LUT :param ndim: Dimension of LUT
:param x_cell: Interpolation grid (nm_nodes x n_dim) :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 :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, solar irradiance, spherical albedo and transmittance interpolated at vtest
""" """
lim = np.zeros((2, ndim), np.int8) lim = np.zeros((2, ndim), np.int8)
# TODO compute that on the full cube
for ii in range(ndim): for ii in range(ndim):
if vtest[ii] >= xnodes[:, ii].max(): if vtest[ii] >= xnodes[:, ii].max():
vtest[ii] = 0.99 * xnodes[-1, ii] 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): ...@@ -297,9 +295,105 @@ def interpol_lut_c(lut1, lut2, xnodes, nm_nodes, ndim, x_cell, vtest, intp_wvl):
for i in range(ndim): for i in range(ndim):
vtest[i] = (vtest[i] - xnodes[lim[0, i], i]) / (xnodes[lim[1, i], i] - xnodes[lim[0, i], i]) 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, :] 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) 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) f_int = np.sum(weights * lut_cell, axis=1)
return f_int 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, ...@@ -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): for jj, cwv in enumerate(cwvs):
vtest = np.asarray([vza, sza, hsf, aot, raa, cwv]) 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_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
......
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