Retrieval 1: Static parameters
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:
a set of constant RT1 parameters from an incidence-angle dependent \(\sigma^0\) dataset.
Show code cell source
%matplotlib widget
from rt1_model import RT1, surface, volume, set_loglevel
from scipy.optimize import least_squares
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 = True # calculate values in decibel
sig0 = True # calculate sigma0 values rather than intensities
noise_sigma = 0.5 if dB is True else 1e-3 # Noise-level (sigma of gaussian noise)
inc = rand.normal(45, 10, (1000,)).clip(20, 70) # Incidence angles
sim_params = dict(tau=0.3, omega=0.4, N=0.1, t_s=0.4) # Simulation parameter values
Set start values and boundaries for the fit
start_vals = dict(tau=0.1, omega=0.2, N=0.3, t_s=0.1)
bnd_vals = dict(tau=(0.01, 0.5), omega=(0.01, 0.5), N=(0.01, 0.5), t_s=(0.01, 0.5))
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))
R.update_params(**sim_params)
tot, surf, vol, inter = R.calc()
tot += rand.normal(0, noise_sigma, tot.size) # Add some random noise
Show code cell output
08:30:18.610 INFO: Evaluating coefficients for interaction-term...
08:30:18.757 INFO: Coefficients extracted, it took 0.02121 sec.
Setup scipy optimize to fit RT1 model to the data
param_names = list(sim_params)
def parse_params(x):
"""Map 1D parameter array to dict {parameter_name: value(s)}."""
return dict(zip(param_names, x))
def fun(x):
"""Calculate residuals."""
R.update_params(**parse_params(x))
res = (R.calc()[0] - tot).ravel()
return res
def jac(x):
"""Calculate jacobian."""
R.update_params(**parse_params(x))
jac = R.jacobian(param_list=list(param_names), format="scipy_least_squares")
return jac
# Unpack start-values and boundaries as required by scipy optimize
x0 = [start_vals[key] for key in param_names]
bounds = list(zip(*[bnd_vals[key] for key in param_names]))
res = least_squares(
fun=fun,
x0=x0,
bounds=bounds,
jac=jac,
ftol=1e-8,
gtol=1e-8,
xtol=1e-3,
verbose=2,
)
found_params = dict(zip(param_names, res.x))
Show code cell output
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 7.5142e+03 1.45e+04
1 2 1.4744e+03 6.04e+03 1.45e-01 4.40e+03
2 3 5.1737e+02 9.57e+02 9.30e-02 1.55e+03
3 4 2.4675e+02 2.71e+02 1.07e-01 7.36e+02
4 5 1.4910e+02 9.76e+01 6.22e-02 2.92e+02
5 6 1.3132e+02 1.78e+01 5.69e-02 1.33e+02
6 8 1.2941e+02 1.91e+00 1.84e-02 4.03e+01
7 9 1.2753e+02 1.88e+00 3.18e-02 1.10e+01
8 11 1.2701e+02 5.23e-01 1.36e-02 1.48e+01
9 12 1.2606e+02 9.44e-01 2.84e-02 3.15e+01
10 13 1.2509e+02 9.71e-01 5.78e-02 3.14e+01
11 15 1.2505e+02 4.34e-02 4.51e-03 2.18e+01
12 16 1.2502e+02 3.23e-02 3.89e-03 3.05e+00
13 18 1.2502e+02 8.79e-04 1.60e-03 4.90e+00
14 19 1.2501e+02 3.79e-03 3.81e-04 7.63e-01
`xtol` termination condition is satisfied.
Function evaluations 19, initial cost 7.5142e+03, final cost 1.2501e+02, first-order optimality 7.63e-01.
| Parameter | Target value | Start value | Retrieved value | (Target - Retrieved) |
|---|---|---|---|---|
| tau | 0.300 | 0.100 | 0.341 | 0.041 |
| omega | 0.400 | 0.200 | 0.380 | -0.020 |
| N | 0.100 | 0.300 | 0.106 | 0.006 |
| t_s | 0.400 | 0.100 | 0.412 | 0.012 |
Initialize analyzer widget and overlay results
Show code cell source
analyze_params = {key: (0.01, 0.5, found_params[key]) for key in param_names}
ana = R.analyze(**analyze_params)
# Plot fit-data on top
ana.ax.scatter(inc, tot, s=10, c="k")
ana.ax.plot(
np.rad2deg(R.t_0).squeeze(),
R.calc(**sim_params)[0].squeeze(),
c="r",
ls="--",
lw=0.5,
zorder=0,
)
ana.ax.plot(
np.rad2deg(R.t_0).squeeze(),
R.calc(**found_params)[0].squeeze(),
c="C0",
ls="--",
lw=0.5,
zorder=0,
)
# Set limits to fit-data range
ana.ax.set_xlim(inc.min() - 2, inc.max() + 2)
ana.ax.set_ylim(tot.min() - 2, tot.max() + 2)
# Indicate fit-results in slider-axes
for key, s in ana.sliders.items():
if key in sim_params:
s.ax.plot(sim_params[key], np.mean(s.ax.get_ylim()), marker="o")
# Add text
t = ana.f.text(
0.6,
0.95,
"\n".join(
[
f"{key:>8} = {found:.3f} ({sim_params[key]:.2f}) "
rf"| $\Delta$ ={found - sim_params[key]: .3f}"
for (key, found) in found_params.items()
]
),
va="top",
fontdict=dict(family="monospace", size=8),
)
08:30:19.265 INFO: Evaluating coefficients for interaction-term...
08:30:19.335 INFO: Coefficients extracted, it took 0.02138 sec.