Retrieval 5: Timeseries with auxiliary data

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 both static and dynamic parameters from a series of incidence-angle dependent \(\sigma^0\) measurements.
use auxiliary timeseries of (incidence-angle dependent) parameter values as input

Specify simulation and fit parameters

Set parameter values that are used to simulate the data

dB, sig0 = True, 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.25, 0.25, (num, 1)).clip(0.01, 0.5)  # NormBRDF values

tau = 0.2 + 1.5 * np.sin(np.linspace(0, 2.*np.pi, num))**2
tau = tau[:,np.newaxis] # broadcast the tau-values among all incidence-angles

bsf = np.cos(np.linspace(0, 2.*np.pi, num))**2
# broadcast the tau-values among all incidence-angles
bsf = np.tile(bsf[:,np.newaxis], incs)* np.linspace(0.3, 1, incs)

sim_params = dict(omega=0.1, N=N)  # Simulation parameter values
const_params = dict(t_s=0.4, tau=tau, bsf=bsf)  # Constant parameters (assumed to be known)

Visualize used auxiliary datasets

Set start values and boundaries for the fit

start_vals = dict(omega=0.2, N=[0.3] * num)
bnd_vals = dict(omega=(0.01, 0.5), N=[(0.01, 0.5)] * num)

Setup RT1 and create a simulated dataset

V = volume.Isotropic()
SRF = surface.HG_nadirnorm(t="t_s", ncoefs=10)

R = RT1(V=V, SRF=SRF, int_Q=True, dB=dB, sig0=True)
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, **const_params)

tot = R.calc()[0]
tot += rand.normal(0, noise_sigma, tot.shape)  # Add some random noise

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(omega=x[0], N=x[1:][:, np.newaxis])

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["omega"], *start_vals["N"]]
bounds = list(zip(*[bnd_vals["omega"], *bnd_vals["N"]]))

res = least_squares(
    fun=fun,
    x0=x0,
    bounds=bounds,
    jac=jac,
    ftol=1e-8,
    gtol=1e-8,
    xtol=1e-3,
    verbose=2,
)

# Unpack found parameters
found_params = parse_params(res.x)
# Calcuate total backscatter based on found parameters
found_tot = R.calc(**found_params)[0]

Retrieved Parameters
Parameter	Target value	Start value	Retrieved value	(Target - Retrieved)
omega	0.100	0.200	0.100	-0.000
N (mean)	0.251	0.300	0.251	0.000

Visualize Results

Plot timeseries

Initialize analyzer widget and overlay results

08:30:42.296 INFO: Evaluating coefficients for interaction-term...

08:30:42.338 INFO: Coefficients extracted, it took 0.01001 sec.