3.4.4 Weighted interpolation

Table of Contents

  • 3.4.4 Weighted interpolation

    • 3.4.4.1 Smoothing

    • 3.4.4.2 weighted smoothing

3.4.4.1 Smoothing

There are many approaches to weighted interpolation. One such is to use a convolution smoothing operation.

In convolution, we combine a signal \(y\) with a filter \(f\) to achieve a filtered signal. For example, if we have an noisy signal, we will attempt to reduce the influence of high frequency information in the signal (a ‘low pass’ filter, as we let the low frequency information pass).

In this approach, we define a digital filter (convolution filter) that should be some sort of weighted average function. A typical filter is the Gaussian, defined by the parameter \(\sigma\). The larger the value of \(\sigma\), the ‘wider’ the filter, which means that the weighted average will be takjen over a greater extent.

We do not expect you to be overly concerned with the code in these sections below, as the main effort should be directed at understanding smoothing and interpolation at this point.

Let’s look at this filter:

from IPython.display import HTML
from geog0111.demofilt1 import demofilt1

HTML(demofilt1().to_jshtml())


Once Loop Reflect

A digital (discrete) convolution uses a sampled signal and filter (represented on a grid).

Without going into the maths, the convolution operates by running the filter centred on \(x_c\) over the extent of the signal as illustrated below.

At each value of :math:`x_c`, the filter is multiplied by the signal. Clearly, where the filter is zero, the influence of the signal is zero. So, the filter effectively selects a local window of data points (shown as green crosses below). The result, i.e. the filtered signal, is simply the weighted average of these local samples. This can be seen in the lower panel, where the green dot shows the filtered signal at \(x_c\).

from IPython.display import HTML
from geog0111.demofilt3 import demofilt3

HTML(demofilt3().to_jshtml())


Once Loop Reflect
_images/Chapter3_4a_GDAL_stacking_and_interpolating-convolution_5_1.png

The convolution of a low pass (‘smoothing’) filter with a signal results in a smoothing of the signal:

from IPython.display import HTML
from geog0111.demofilt2 import demofilt2

HTML(demofilt2().to_jshtml())


Once Loop Reflect
_images/Chapter3_4a_GDAL_stacking_and_interpolating-convolution_7_1.png

The degree of smoothing is greater, the larger the value of \(\sigma\) (i.e. the broader the convolution filter). Eventually, the signal would go ‘flat’, i.e. be a constant (weighted mean) value, if we used a very broad filter.

For small values of \(\sigma\) though, the signal is still rather noisy.

There is a trade-off then between noise suppression and the degree of generalisation.

We could define some concept of an ‘optimium’ value of \(\sigma\), for example using cross validation. But often we can decide empirically an appropriate filter width.

We perform the convolution using:

convolve1d(y, filter,mode='wrap')

from the scipy.ndimage.filters library.

The first argument is the signal to be filtered. The second is the filter (the Gaussian here).

We set

mode='wrap'

which defines the boundary conditions to be periodic (‘wrapped around’). This would generally be appropriate for time series e.g. of one year.

3.4.4.2 weighted smoothing

In the convolution examples above, we apply the smoothing equally to all samples of the signal.

What if we knew that some samples were ‘better’ (quality) than others? What if some samples were missing? How could we incorporate this information?

The answer is to weight the convolution.

For each sample point in the signal, we define a weight, where the weight is high if we trust the data point and low if we don’t trust it much. We apply a weight of zero if a data point is missing.

Then, the weighted convolution is simply the result of applying the filter to (weight \(\times\) signal), divided by the result of applying the filter to the weight alone. You can think of the denominator here as a form of ‘re-normalisation’ of the filter.

We illustrate this by removing samples from the example above, and giving a weight of zero for the ‘missing’ observations.

from IPython.display import HTML
from geog0111.demofilt4 import demofilt4

HTML(demofilt4().to_jshtml())


Once Loop Reflect
_images/Chapter3_4a_GDAL_stacking_and_interpolating-convolution_10_1.png

The result is not as good as we got above, but that is hardly surprising: we are interpolating here and it is hard to interpolate over large gaps.

Since the gaps are quite large, we might benefit from using a larger filter extent. This is then liable to degrade the quality somewhat (over smooth) in data rich areas. These are typical trade-offs we must balance.

from IPython.display import HTML
from geog0111.demofilt5 import demofilt5

HTML(demofilt5().to_jshtml())


Once Loop Reflect
_images/Chapter3_4a_GDAL_stacking_and_interpolating-convolution_12_1.png