#! /usr/bin/env python3
#
# Copyright 2018 California Institute of Technology
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# ISOFIT: Imaging Spectrometer Optimal FITting
# Author: David R Thompson, david.r.thompson@jpl.nasa.gov
#
from __future__ import annotations
import numpy as np
from scipy.linalg import block_diag
from isofit.core.common import eps, svd_inv_sqrt
from isofit.surface.surface_multicomp import MultiComponentSurface
[docs]
class GlintModelSurface(MultiComponentSurface):
"""A model of the surface based on a collection of multivariate
Gaussians, extended with two surface glint terms (sun + sky glint)."""
def __init__(self, full_config: Config):
super().__init__(full_config)
config = full_config.forward_model.surface
# TODO: Enforce this attribute in the config, not here (this is hidden)
self.statevec_names.extend(["SKY_GLINT", "SUN_GLINT"])
[docs]
self.glint_ind = len(self.statevec_names) - 2
[docs]
self.sky_glint_ind = self.glint_ind
[docs]
self.sun_glint_ind = self.glint_ind + 1
self.scale.extend([1.0, 1.0])
# Numbers from Marcel Koenig; used for prior mean
self.init.extend([1 / np.pi, 0.02])
# Gege (2021), WASI user manual
self.bounds.extend([[0, 10], [0, 10]])
[docs]
self.n_state = self.n_state + 2
# Useful indexes to track
self.glint_ind = len(self.statevec_names) - 2
self.sky_glint_ind = self.glint_ind
self.sun_glint_ind = self.glint_ind + 1
[docs]
self.idx_surface = np.arange(len(self.statevec_names))
# Change this if you don't want to analytical solve for all the full statevector elements.
[docs]
self.analytical_iv_idx = np.arange(len(self.statevec_names))
[docs]
self.sun_glint_sigma = (
config.sun_glint_prior_sigma * np.array(self.scale[self.sun_glint_ind])
) ** 2
[docs]
self.sky_glint_sigma = (
config.sky_glint_prior_sigma * np.array(self.scale[self.sky_glint_ind])
) ** 2
# Compute and and cache normalized Sa inversions for glint case
Cov = np.array([[self.sky_glint_sigma, 0], [0, self.sun_glint_sigma]])
self.Sa_inv_glint, self.Sa_inv_sqrt_glint = svd_inv_sqrt(
Cov / np.mean(np.diag(Cov))
)
[docs]
def xa(self, x_surface, geom):
"""Mean of prior distribution, calculated at state x."""
mu = MultiComponentSurface.xa(self, x_surface, geom)
mu[self.glint_ind :] = self.init[self.glint_ind :]
return mu
[docs]
def Sa(self, x_surface, geom):
"""Covariance of prior distribution, calculated at state x. We find
the covariance in a normalized space (normalizing by z) and then un-
normalize the result for the calling function."""
Sa_unnormalized, Sa_inv_normalized, Sa_inv_sqrt_normalized = (
MultiComponentSurface.Sa(self, x_surface, geom)
)
Sa_unnormalized[self.sun_glint_ind, self.sun_glint_ind] = self.sun_glint_sigma
Sa_unnormalized[self.sky_glint_ind, self.sky_glint_ind] = self.sky_glint_sigma
# Append normalized Sa inv and sqrt from glint model
Sa_inv_normalized = block_diag(Sa_inv_normalized, self.Sa_inv_glint)
Sa_inv_sqrt_normalized = block_diag(
Sa_inv_sqrt_normalized, self.Sa_inv_sqrt_glint
)
return Sa_unnormalized, Sa_inv_normalized, Sa_inv_sqrt_normalized
[docs]
def fit_params(self, rfl_meas, geom, *args):
"""Given a reflectance estimate and one or more emissive parameters,
fit a state vector.
"""
# Make sure the glint estimation is not
# driving short wavelengths into the negative
# causes issues in inversion
blue_band_1 = np.argmin(np.abs(450 - self.wl))
blue_band_2 = np.argmin(np.abs(500 - self.wl))
# Estimate additive glint. Uses multiple options to handle different sensors.
# Try SWIR, then NIR
blue_band_1 = np.argmin(np.abs(450 - self.wl))
blue_band_2 = np.argmin(np.abs(500 - self.wl))
if np.max(self.wl) >= 2300:
glint_band_1 = np.argmin(np.abs(2200 - self.wl))
glint_band_2 = np.argmin(np.abs(2300 - self.wl))
else:
glint_band_1 = np.argmin(np.abs(1000 - self.wl))
glint_band_2 = np.argmin(np.abs(1020 - self.wl))
glint_est = np.min(
[
np.median(rfl_meas[blue_band_1:blue_band_2]),
np.median(rfl_meas[glint_band_1:glint_band_2]),
]
)
# hard-coded glint bounds from experience #TODO - get from config
bounds_glint_est = [
0,
1.0,
]
glint_est = max(
bounds_glint_est[0] + eps,
min(bounds_glint_est[1] - eps, glint_est),
)
lamb_est = rfl_meas - glint_est
x = MultiComponentSurface.fit_params(self, lamb_est, geom) # Bounds reflectance
# Get estimate for g_dd and g_dsf parameters
# Set sky glint to a static number; use prior mean.
g_dsf_est = self.init[self.sky_glint_ind]
g_dd_est = glint_est / self.fresnel_rf(geom.observer_zenith)
# Updating self.init will set the prior mean (xa) to this value
self.init[self.sun_glint_ind] = g_dd_est
# SKY_GLINT g_dsf
x[self.sky_glint_ind] = g_dsf_est
# SUN_GLINT g_dd
x[self.sun_glint_ind] = g_dd_est
return x
[docs]
def calc_rfl(self, x_surface, geom):
"""Direct and diffuse Reflectance (includes sun and sky glint).
Inputs:
x_surface : np.ndarray
Surface portion of the statevector element
geom : Geometry
Isofit geometry object
Outputs:
rho_dir_dir : np.ndarray
Reflectance quantity for downward direct photon paths
rho_dif_dir : np.ndarray
Reflectance quantity for downward diffuse photon paths
NOTE:
Here, we treat direct and diffuse photon path reflectance
differently. The sun and sky glint magnitudes are statevector
elements that interact with the two reflectance quantities
independently.
"""
# fresnel reflectance factor (approx. 0.02 for nadir view)
rho_ls = self.fresnel_rf(geom.observer_zenith)
# Enforce bounds, also turns -9999. fill into 0.
# Note if null_value > bounds[1], will return bounds[1]
sun_glint = (
np.max(
[
self.bounds[self.sun_glint_ind][0],
np.min(
[
self.bounds[self.sun_glint_ind][1],
x_surface[self.sun_glint_ind],
]
),
]
)
* rho_ls
)
sky_glint = (
np.max(
[
self.bounds[self.sky_glint_ind][0],
np.min(
[
self.bounds[self.sky_glint_ind][1],
x_surface[self.sky_glint_ind],
]
),
]
)
* rho_ls
)
rho_dir_dir = self.calc_lamb(x_surface, geom) + sun_glint
rho_dif_dir = self.calc_lamb(x_surface, geom) + sky_glint
return rho_dir_dir, rho_dif_dir
[docs]
def drfl_dsurface(self, x_surface, geom):
"""Partial derivative of reflectance with respect to state vector,
calculated at x_surface.
"""
drfl = self.dlamb_dsurface(x_surface, geom)
rho_ls = self.fresnel_rf(geom.observer_zenith)
# TODO make the indexing better for the surface state elements
drfl[:, self.sun_glint_ind] = rho_ls
drfl[:, self.sky_glint_ind] = rho_ls
return drfl
[docs]
def drdn_dglint(self, L_tot, L_down_dir, s_alb, rho_dif_dir):
"""Derivative of radiance with respect to
the direct and diffuse glint terms"""
drdn_dgdd = L_down_dir
drdn_dgdsf = (L_tot / ((1.0 - (s_alb * rho_dif_dir)) ** 2)) - drdn_dgdd
return drdn_dgdd, drdn_dgdsf
[docs]
def drdn_dsurface(
self,
rho_dif_dir,
drfl_dsurface,
dLs_dsurface,
s_alb,
t_total_up,
L_tot,
L_down_dir,
):
"""Derivative of radiance with respect to
full surface vector"""
# Construct the output matrix
# Dimensions should be (len(RT.wl), len(x_surface))
# which is correctly handled by the instrument resampling
drdn_dsurface = np.zeros(drfl_dsurface.shape)
drdn_drfl = self.drdn_drfl(L_tot, s_alb, rho_dif_dir)
drdn_dsurface[:, : self.n_wl] = np.multiply(
drdn_drfl[:, np.newaxis], drfl_dsurface[:, : self.n_wl]
)
# Glint derivatives
drdn_dgdd, drdn_dgdsf = self.drdn_dglint(L_tot, L_down_dir, s_alb, rho_dif_dir)
# Store the glint derivatives as last two rows in drdn_drfl
drdn_dsurface[:, self.sun_glint_ind] = (
drdn_dgdd * drfl_dsurface[:, self.sun_glint_ind]
)
drdn_dsurface[:, self.sky_glint_ind] = (
drdn_dgdsf * drfl_dsurface[:, self.sky_glint_ind]
)
# Get the derivative w.r.t. surface emission
drdn_dLs = np.multiply(self.drdn_dLs(t_total_up)[:, np.newaxis], dLs_dsurface)
return np.add(drdn_dsurface, drdn_dLs)
[docs]
def analytical_model(
self,
background,
L_down_dir,
L_down_dif,
L_tot,
geom,
L_dir_dir=None,
L_dir_dif=None,
L_dif_dir=None,
L_dif_dif=None,
):
"""
Linearization of the glint terms to use in AOE inner loop.
Function will fetch the linearization of the rho terms and
add the matrix components for the direct glint term.
Currently we set the diffuse glint scaling term to constant
value, which makes the AOE inner loop inversion possible.
"""
rho_ls = self.fresnel_rf(geom.observer_zenith)
# Construct the H matrix from:
# theta (rho portion)
# gam (sun glint portion)
# ep (sky glint portion)
H = super().analytical_model(
background,
L_down_dir,
L_down_dif,
L_tot,
geom,
L_dir_dir,
L_dir_dif,
L_dif_dir,
L_dif_dif,
)
# NOTE: The order of ep and gam respectively is important
# It must match the alphabeitcal order of the glint terms
# Diffuse portion
ep = (
(L_dif_dir + L_dif_dif) + ((L_tot * background) / (1 - background))
) * rho_ls
# If you ignore multi-scattering
# ep = (L_dif_dir + L_dif_dif) * rho_ls
ep = np.reshape(ep, (len(ep), 1))
H = np.append(H, ep, axis=1)
# Direct portion
gam = (L_dir_dir + L_dir_dif) * rho_ls
gam = np.reshape(gam, (len(gam), 1))
H = np.append(H, gam, axis=1)
return H
[docs]
def summarize(self, x_surface, geom):
"""Summary of state vector."""
return MultiComponentSurface.summarize(
self, x_surface, geom
) + " Sun Glint: %5.3f, Sky Glint: %5.3f" % (
x_surface[self.sun_glint_ind],
x_surface[self.sky_glint_ind],
)
@staticmethod
[docs]
def fresnel_rf(vza):
"""Calculates reflectance factor of sky radiance based on the
Fresnel equation for unpolarized light as a function of view zenith angle (vza).
"""
if vza > 0.0:
n_w = 1.33 # refractive index of water
theta = np.deg2rad(vza)
# calculate angle of refraction using Snell′s law
theta_i = np.arcsin(np.sin(theta) / n_w)
# reflectance factor of sky radiance based on the Fresnel equation for unpolarized light
rho_ls = 0.5 * np.abs(
((np.sin(theta - theta_i) ** 2) / (np.sin(theta + theta_i) ** 2))
+ ((np.tan(theta - theta_i) ** 2) / (np.tan(theta + theta_i) ** 2))
)
else:
rho_ls = 0.02 # the reflectance factor converges to 0.02 for view angles equal to 0.0°
return rho_ls