Retrieval 3: Multi-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 from a series of incidence-angle dependent \(\sigma^0\) measurements.
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
start_vals = dict(omega=[0.2], tau=[0.3] * num, N=[0.1] * num)
bnd_vals = dict(omega=[(0.01, 0.5)], tau=[(0.01, 1.)] * num, 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:25.711 INFO: Evaluating coefficients for interaction-term...
08:30:25.849 INFO: Coefficients extracted, it took 0.02155 sec.
Setup scipy optimize to fit RT1 model to the data
def parse_params(x):
"""Map 1D parameter array to dict {parameter_name: value(s)}."""
return dict(omega=x[0], tau=x[1:num+1][:, np.newaxis], N=x[num+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
def jac(x):
"""Calculate jacobian."""
R.update_params(**parse_params(x), **const_params)
jac = R.jacobian(param_list=["omega", "tau", "N"], format="scipy_least_squares")
return jac
# 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]
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 3.2617e+00 1.28e+00
1 2 5.8472e-01 2.68e+00 1.12e+00 5.62e-01
2 3 1.8169e-01 4.03e-01 2.19e+00 1.61e+00
3 4 6.5275e-02 1.16e-01 1.27e+00 6.60e-01
4 5 1.6658e-02 4.86e-02 9.11e-01 2.13e-01
5 6 1.0778e-02 5.88e-03 9.29e-01 3.51e-01
6 7 2.7516e-03 8.03e-03 2.30e-01 1.27e-01
7 8 1.6904e-03 1.06e-03 1.87e-01 2.69e-02
8 9 1.6264e-03 6.40e-05 3.55e-02 3.17e-02
9 10 1.5583e-03 6.81e-05 4.71e-02 3.54e-02
10 11 1.4991e-03 5.92e-05 5.76e-02 3.63e-02
11 12 1.4618e-03 3.72e-05 4.99e-02 3.45e-02
12 13 1.4368e-03 2.51e-05 4.27e-02 3.11e-02
13 14 1.4222e-03 1.46e-05 4.11e-02 2.91e-02
14 15 1.4128e-03 9.47e-06 3.62e-02 2.60e-02
15 16 1.4070e-03 5.77e-06 2.66e-02 2.32e-02
16 17 1.4032e-03 3.77e-06 2.54e-02 2.14e-02
17 18 1.4004e-03 2.79e-06 2.52e-02 2.04e-02
18 19 1.3934e-03 7.06e-06 2.09e-04 8.70e-04
19 22 1.3932e-03 1.49e-07 6.14e-03 3.19e-03
20 23 1.3919e-03 1.36e-06 6.24e-03 2.50e-03
21 24 1.3916e-03 2.97e-07 1.25e-02 5.64e-03
22 25 1.3904e-03 1.19e-06 1.25e-02 5.33e-03
23 27 1.3896e-03 8.07e-07 6.16e-03 2.49e-03
24 28 1.3894e-03 1.98e-07 1.25e-02 5.64e-03
25 29 1.3884e-03 9.88e-07 1.25e-02 5.27e-03
26 31 1.3876e-03 7.53e-07 5.40e-03 2.16e-03
27 32 1.3875e-03 1.57e-07 1.09e-02 4.97e-03
28 33 1.3867e-03 7.38e-07 1.10e-02 4.64e-03
29 35 1.3861e-03 5.77e-07 4.73e-03 1.91e-03
30 36 1.3860e-03 1.07e-07 9.53e-03 4.40e-03
31 37 1.3855e-03 5.63e-07 9.62e-03 4.11e-03
32 38 1.3853e-03 1.38e-07 8.41e-03 4.21e-03
33 39 1.3849e-03 3.98e-07 9.38e-03 3.80e-03
34 40 1.3849e-03 3.25e-08 7.00e-03 4.21e-03
35 41 1.3843e-03 6.55e-07 8.12e-04 2.54e-04
36 42 1.3842e-03 9.57e-08 1.69e-03 8.27e-04
37 43 1.3840e-03 1.22e-07 3.46e-03 1.45e-03
38 44 1.3840e-03 6.75e-08 6.95e-03 3.20e-03
39 45 1.3837e-03 2.78e-07 6.94e-03 3.00e-03
40 46 1.3836e-03 1.14e-07 6.78e-03 3.20e-03
41 47 1.3834e-03 1.56e-07 6.89e-03 2.84e-03
42 49 1.3831e-03 3.06e-07 1.73e-03 6.20e-04
43 50 1.3831e-03 5.52e-08 3.52e-03 1.69e-03
44 51 1.3830e-03 9.39e-08 3.51e-03 1.58e-03
45 52 1.3829e-03 7.79e-08 3.44e-03 1.57e-03
46 53 1.3828e-03 6.66e-08 3.46e-03 1.59e-03
47 54 1.3828e-03 6.41e-08 3.44e-03 1.57e-03
48 55 1.3827e-03 5.06e-08 3.44e-03 1.54e-03
49 56 1.3827e-03 3.11e-08 3.44e-03 1.58e-03
50 57 1.3826e-03 8.83e-08 2.95e-04 1.19e-04
51 58 1.3826e-03 2.23e-08 8.51e-04 3.43e-04
52 59 1.3825e-03 2.26e-08 1.74e-03 8.02e-04
53 60 1.3825e-03 1.27e-08 3.49e-03 1.67e-03
54 61 1.3825e-03 5.35e-08 2.99e-04 2.00e-05
55 62 1.3825e-03 1.50e-08 8.33e-04 3.75e-04
56 63 1.3824e-03 1.84e-08 1.74e-03 8.29e-04
57 64 1.3824e-03 6.34e-09 3.48e-03 1.67e-03
58 65 1.3824e-03 4.69e-08 2.77e-04 1.47e-05
59 66 1.3824e-03 1.19e-08 8.34e-04 4.41e-04
`ftol` termination condition is satisfied.
Function evaluations 66, initial cost 3.2617e+00, final cost 1.3824e-03, first-order optimality 4.41e-04.
| Parameter | Target value | Start value | Retrieved value | (Target - Retrieved) |
|---|---|---|---|---|
| omega | 0.200 | 0.200 | 0.200 | 0.000 |
| tau (mean) | 0.347 | 0.300 | 0.347 | -0.000 |
| N (mean) | 0.108 | 0.100 | 0.108 | 0.000 |
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:29.907 INFO: Evaluating coefficients for interaction-term...
08:30:29.978 INFO: Coefficients extracted, it took 0.02132 sec.