Compute residuals to the timeseries inversion, save to raster outputs

[scottstanie]

To see whether the time series inversion ran smoothly (which encompasses detecting unwrapping errors, as well as several other upstream errors which may have snuck in), we'd like to see the residuals of

$$ \mathbf{r} := | A\mathbf{x} - \mathbf{b}| $$

where $A$ is the incidence matrix used during the inversion, $\mathbf{b}$ is the vector of unwrapped phases at a single pixel, and $\mathbf{x}$ is the vector of solved phases, one per secondary date.

Since there could be one single date (or a couple sparse dates) which had problems, we the same number of residual rasters as timeseries rasters by doing

$$ |A|^{T} \mathbf{r} $$

where $|A|^{T}$ is the elementwise absolute value of A.

Small test over tests on the walkthrough data:

A histogram of all the residuals for this area (with vertical lines marking multiple of pi)

reader = io.RasterStackReader.from_file_list( sorted(Path("work-walkthrough/timeseries/").glob("residu*")))
residuals = reader[:, :, :]
plt.hist(residuals.ravel(), log=True, bins=30)
for i in range(1, 8):
    plt.vlines(i * np.pi, 0, 1e6, color='k', label=f'{i}pi')