Retrieval 6: Multi-temporal parameter timeseries
About this retrieval example
This example shows how to use the rt1 python package together with scipy optimize to setup a retrieval procedure to
obtain multiple dynamic parameters with different temporal frequencies from a series of incidence-angle dependent \(\sigma^0\) measurements.
IMPORTANT Make sure to checkout “Retrieval 3: Multi-parameter timeseries” first!
This example will use the same start conditions, buttauvalues are retrieved at a different temporal frequency compared toN.
Show code cell source
%matplotlib widget
from rt1_model import RT1, surface, volume, set_loglevel
from scipy.optimize import least_squares
import matplotlib.pyplot as plt
import numpy as np
rand = np.random.RandomState(123456) # initialize a reproducible random state
set_loglevel("info")
Specify simulation and fit parameters
Set parameter values that are used to simulate the data
dB, sig0 = False, True
num = 100 # Number of measurements
incs = 30 # Available incidence angles per measurement
noise_sigma = 0.5 if dB is True else 1e-3 # Noise-level (sigma of gaussian noise)
inc = rand.normal(45, 10, (num, incs)).clip(20, 70) # Incidence angles
N = rand.normal(0.1, 0.1, (num, 1)).clip(0.01, 0.25) # NormBRDF values
tau = (0.1 + 0.5 * np.sin(np.linspace(0, 2*np.pi, num))**2)[:,np.newaxis] # Optical Depth values
sim_params = dict(tau=tau, omega=0.2, N=N) # Simulation parameter values
const_params = dict(t_s=0.4) # Constant parameters (assumed to be known)
Show code cell source
f, ax = plt.subplots(figsize=(10, 2))
f.canvas.header_visible = False
f.suptitle(r"Input timeseries for 'N' and 'tau' parameters")
ax.plot(N, marker=".", lw=0.5, label="N")
ax.plot(tau, marker=".", lw=0.5, label="tau")
ax.set_ylabel("N / tau")
ax.set_xlabel("# measurement")
ax.legend(loc="upper left")
f.tight_layout()
Set start values and boundaries for the fit
Note
For the parameter
omega, a single value is optimized for the entire timeseries.For the parameter
N, a unique value is optimized for each timestamp.For the parameter
taua unique value is optimized only everytau_freqtimestamps!
(assuming thattauremains constant during the period)
# set the temporal frequency at which "tau" values should be retrieved
tau_freq = 5
assert num%tau_freq==0, f"The length of the timeseries (num={num}) must be divisible by the temporal frequency of 'tau' (tau_freq={tau_freq})!"
n_tau = num // tau_freq # number of unique tau-values with respect to the chosen frequency
start_vals = dict(omega=[0.2], tau=[0.3] * n_tau, N=[0.1] * num)
bnd_vals = dict(omega=[(0.01, 0.5)], tau=[(0.01, 1.)] * n_tau, N=[(0.01, 0.5)] * num)
Setup RT1 and create a simulated dataset
V = volume.Rayleigh()
SRF = surface.HG_nadirnorm(t="t_s", ncoefs=10)
R = RT1(V=V, SRF=SRF, int_Q=True, dB=dB, sig0=sig0)
R.set_monostatic(p_0=0)
R.NormBRDF = "N" # Use a synonym for NormBRDF parameter
R.set_geometry(t_0=np.deg2rad(inc))
tot = R.calc(**sim_params, **const_params)[0]
tot += rand.normal(0, noise_sigma, tot.shape) # Add some random noise
Show code cell output
08:30:44.970 INFO: Evaluating coefficients for interaction-term...
08:30:45.109 INFO: Coefficients extracted, it took 0.02145 sec.
Setup scipy optimize to fit RT1 model to the data
Since we obtain tau values only every tau_freq timestamps, we need to distribute the obtained values to the full timeseries!
def parse_params(x):
"""Map 1D parameter array to dict {parameter_name: value(s)}."""
return dict(
omega=x[0],
tau=np.repeat(x[1:n_tau + 1], tau_freq)[:, np.newaxis], # make sure to re-distribute tau to the data-index
N=x[n_tau+1:][:, np.newaxis]
)
def fun(x):
"""Calculate residuals."""
R.update_params(**parse_params(x), **const_params)
res = (R.calc()[0] - tot).ravel() # Ravel result because scipy requires 1D arrays
return res
Note
In order to obtain a representation suitable to retrieve tau at a reduced frequency, we have to re-shape the jacobian accordingly!
from scipy.sparse import block_diag, hstack
def jac(x):
"""Calculate jacobian."""
R.update_params(**parse_params(x), **const_params)
# obtain block-diagonal sparse jacobians for "omega" and "N".
jac_omega = R.jacobian(param_list=["omega"], format="scipy_least_squares")
jac_N = R.jacobian(param_list=["N"], format="scipy_least_squares")
# re-shape jacobian for "tau" to desired retrieval frequency
jac_tau = R.jacobian(param_list=["tau"])[0].reshape(n_tau, -1)
jac_tau = block_diag(jac_tau.tolist(), "csr").T
# stack jacobians
return hstack((jac_omega, jac_tau, jac_N), "csr")
# Unpack start-values and boundaries as required by scipy optimize
x0 = [*start_vals["omega"], *start_vals["tau"], *start_vals["N"]]
bounds = list(zip(*[*bnd_vals["omega"], *bnd_vals["tau"], *bnd_vals["N"]]))
res = least_squares(
fun=fun,
x0=x0,
bounds=bounds,
jac=jac,
ftol=1e-5,
gtol=1e-5,
xtol=1e-5,
verbose=2,
)
# Unpack found parameters
found_params = parse_params(res.x)
# Calcuate total backscatter based on found parameters
found_tot = R.calc(**found_params, **const_params)[0]
Show code cell output
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 3.2617e+00 1.28e+00
1 2 6.4669e-01 2.61e+00 8.01e-01 2.92e-01
2 3 1.6216e-01 4.85e-01 5.50e-01 7.38e-01
3 4 4.6368e-02 1.16e-01 7.45e-01 2.95e-01
4 5 1.3026e-02 3.33e-02 2.38e-01 1.64e-01
5 7 1.0593e-02 2.43e-03 9.33e-02 5.29e-02
6 8 8.5685e-03 2.02e-03 1.85e-01 1.68e-01
7 9 7.4982e-03 1.07e-03 1.59e-01 1.52e-01
8 10 6.9020e-03 5.96e-04 8.01e-02 9.42e-02
9 11 6.7323e-03 1.70e-04 8.00e-02 8.43e-02
10 12 6.5787e-03 1.54e-04 6.66e-02 7.28e-02
11 13 6.4829e-03 9.58e-05 5.75e-02 6.87e-02
12 14 6.4130e-03 7.00e-05 5.04e-02 6.34e-02
13 15 6.3663e-03 4.66e-05 4.35e-02 5.73e-02
14 16 6.3383e-03 2.81e-05 3.50e-02 5.01e-02
15 17 6.3198e-03 1.85e-05 2.80e-02 4.15e-02
16 18 6.3062e-03 1.36e-05 1.94e-02 3.12e-02
17 19 6.2992e-03 7.06e-06 1.20e-02 2.27e-02
18 20 6.2877e-03 1.14e-05 2.80e-04 1.02e-03
19 25 6.2877e-03 4.36e-08 4.50e-05 3.58e-04
`ftol` termination condition is satisfied.
Function evaluations 25, initial cost 3.2617e+00, final cost 6.2877e-03, first-order optimality 3.58e-04.
| Parameter | Target value | Start value | Retrieved value | (Target - Retrieved) |
|---|---|---|---|---|
| omega | 0.200 | 0.200 | 0.207 | 0.007 |
| tau (mean) | 0.347 | 0.300 | 0.311 | -0.037 |
| N (mean) | 0.108 | 0.100 | 0.101 | -0.007 |
Visualize Results
Plot timeseries
Show code cell source
f, (ax, ax2) = plt.subplots(2, figsize=(10, 4), sharex=True)
f.canvas.header_visible = False
# Plot retrieved parameter timeseries
ax.set_ylabel("N / tau")
ax.plot(sim_params["N"], marker=".", lw=0.25, label="target N", c="C0")
ax.plot(found_params["N"], marker="o", lw=0.25, markerfacecolor="none", label="retrieved N", c="C0")
ax.plot(sim_params["tau"], marker=".", lw=0.25, label="target tau", c="C1")
ax.plot(found_params["tau"], marker="o", lw=0.25, markerfacecolor="none", label="retrieved tau", c="C1")
# Plot backscatter timeseries
ax2.set_ylabel(r"$\sigma_0$ [dB]")
ax2.set_xlabel("# measurement")
ax2.plot(tot, lw=0, marker=".", c="C0", ms=3)
ax2.plot(found_tot, lw=0, marker="o", markerfacecolor="none", c="C1", ms=3)
ax.legend(loc="upper center", ncols=3, bbox_to_anchor=(0.5, 1.5))
f.tight_layout()
Initialize analyzer widget and overlay results
Show code cell source
analyze_params = {key: (*np.mean(np.atleast_2d(bnd_vals[key]), axis=0), found_params[key].mean()) for key in found_params}
ana = R.analyze(**analyze_params)
# Plot fit-data on top
ana.ax.scatter(inc, tot, c="k", s=3, zorder=0)
ana.ax.scatter(inc, found_tot, c="C0", s=1, zorder=0)
# Indicate fit-results in slider-axes
for key, s in ana.sliders.items():
if key in ["omega"]:
s.ax.plot(sim_params[key], np.mean(s.ax.get_ylim()), marker="o")
# Add text for static parameters
t = ana.f.text(
0.6,
0.95,
"\n".join(
[
f"{key:>8} = {found_params[key]:.3f} ({sim_params[key]:.2f}) "
rf"| $\Delta$ = {found_params[key] - sim_params[key]: .3f}"
for key in ["omega"]
]
),
va="top",
fontdict=dict(family="monospace", size=8),
)
08:30:46.673 INFO: Evaluating coefficients for interaction-term...
08:30:46.745 INFO: Coefficients extracted, it took 0.02140 sec.