11-interpolation.qmd

# Simple interpolation

```{julia}
#| echo: false
#| output: false
import Pkg
Pkg.activate(".")
using GeoStats
import CairoMakie as Mke
```

A very common task in geospatial data science is **geospatial interpolation**,
i.e., predicting variables on geometries that lie between two or more geometries that
have measurements. In this chapter, we will exploit **geospatial correlation**
to make good predictions of continuous variables over an entire domain based on
sparse measurements, which are usually stored on a `GeometrySet`.

The basic idea behind most geostatistical interpolation methods is weighted
combination of values from neighboring geometries. Given a geometry $u$, we want
to estimate the value of the variable at this geometry $z(u)$ using weighted
combinations of measurements from neighboring geometries $u_i,\ i=1,2,\ldots,n$:

$$
z(u) = \lambda_1 z(u_1) + \lambda_2 z(u_2) + \cdots + \lambda_n z(u_n)
$$

The methods differ in the way they compute the weights $\lambda_i,\ i=1,2,\ldots,n$,
and we will cover two basic methods from classical literature: IDW and Kriging.

## IDW

In Inverse Distance Weighting (IDW), the weights are computed in terms of
distances $d(u, u_i)$ to the neighboring geometries:

$$
\lambda_i = \frac{1}{{d(u, u_i)}^\beta}
$$

This basic idea was proposed by @Shepard1968, who also studied the effect of the
exponent $\beta$ in the interpolation results. Here, we will visualize the results
using synthetic data:

```{julia}
data = georef((z=[1.0, 0.0, 1.0],), [(25, 25), (50, 75), (75, 50)])
```

```{julia}
viewer(data, pointsize = 10)
```

First, we need to the define the domain of interpolation, i.e., the geometries
where we want to estimate the variable `z`. In this case, we will perform
interpolation on a 2D `CartesianGrid`:

```{julia}
grid = CartesianGrid(100, 100)
```

With the measurements of the variable `z` in the geotable, and the domain of
interpolation, we can use the `Interpolate` transform with the `IDW` model:

```{julia}
interp = data |> Interpolate(grid, IDW())
```

```{julia}
interp |> viewer
```

To visualize the effect of the exponent, let's extract the interpolation results
along the line segment between two of the measurements, and visualize it:

```{julia}
seg = Segment((25, 25), (50, 75))
```

```{julia}
z = interp[seg, "z"]
```

```{julia}
Mke.lines(z)
```

We observe that the exponent $\beta=1$ leads to a gradual transition
from the value $z=1$ to the value $z=0$. Let's repeat the process with
increasing values of the exponent:

```{julia}
fig = Mke.Figure()
Mke.Axis(fig[1,1])
for β in [1,2,3,4,5]
  interp = data |> Interpolate(grid, IDW(β))
  Mke.lines!(interp[seg, "z"], label = "β=$β")
end
Mke.axislegend(position = :lb)
Mke.current_figure()
```

The larger is the exponent, the more abrupt is the transition of values
between the two locations. In addition, the IDW solution will converge
to the nearest neighbor solution as $\beta \to \infty$:

```{julia}
data |> Interpolate(grid, IDW(100)) |> viewer
```

Custom distances from [Distances.jl](https://github.com/JuliaStats/Distances.jl)
may be used in place of the `Euclidean` distance to meet specific application
requirements (e.g. `Haversine` distance on the sphere).

## Kriging

In Kriging [@Matheron1971], the weights are computed using **geospatial correlation**.
More specifically, they are the solution to a linear system of equations produced
with a theoretical variogram model $\gamma$:

$$
\begin{bmatrix}
\mathbf{G} & \mathbf{1} \\
\mathbf{1}^\top & 0
\end{bmatrix}
\begin{bmatrix}
\mathbf{\lambda} \\
\nu
\end{bmatrix} =
\begin{bmatrix}
\mathbf{g} \\
1
\end{bmatrix}
$$

where $\mathbf{G}_{ij} = \gamma(u_i, u_j)$ and $\mathbf{g}_i = \gamma(u, u_i)$ and
$\nu$ is the [Lagrange multiplier](https://en.wikipedia.org/wiki/Lagrange_multiplier)
associated with the constraint $\mathbf{1}^\top \mathbf{\lambda} = 1$. The system
of equations above is known as *Ordinary Kriging*, but many other variants are
supported by the framework.

::: {.callout-note}

The book by @Olea1999 is a good resource to learn the different systems of
of equations associated with Kriging interpolation. Names such as *Simple Kriging*
*Ordinary Kriging*, *Universal Kriging* are quite popular.

:::

Unlike `IDW`, the `Kriging` solution is a function of pairwise evaluations of distances
between geometries with measurements, represented in the matrix $\mathbf{G}$. The
pairwise evaluations account for possible redundancy in the measurements, which leads
to improvements in the estimates:

```{julia}
γ = GaussianVariogram(range=30.0)

data |> Interpolate(grid, Kriging(γ)) |> viewer
```

In the previous chapter, we learned how the **range** of the variogram determines
the average size of the "blobs" in the image. Let's illustrate this concept again
for increasing values of this parameter:

```{julia}
fig = Mke.Figure()
Mke.Axis(fig[1,1])
for r in [10,20,30,40,50]
  γ = GaussianVariogram(range=r)
  interp = data |> Interpolate(grid, Kriging(γ))
  Mke.lines!(interp[seg, "z"], label = "range=$r")
end
Mke.axislegend(position = :lb)
Mke.current_figure()
```

The larger is the range, the less abrupt is the transition of values between the
two locations. Similar visualizations can be produced by varying the **sill**, the
**nugget** and the **model** of the variogram.

## Example

In order to solidify the concepts learned so far, let's look into an example.
We will cover all the steps that a geospatial data scientist has to perform to
extract **geospatial correlation** from samples and to use this information in
**geospatial interpolation**.

Let's consider an image of the [Walker Lake](https://en.wikipedia.org/wiki/Walker_Lake_(Nevada))
by @Mariethoz2014 as groundtruth. To avoid visualization of large images with CairoMakie.jl,
we will consider a subdomain within a `Box`:


```{julia}
using GeoIO

img = GeoIO.load("data/walkerlake.gslib")

img = img[Box((0, 0), (200, 200)), :]
```

```{julia}
img |> viewer
```

Let's assume that we only have access to 10000 samples from the image:

```{julia}
using Random

samples = img |> Sample(10000, replace=false, rng=MersenneTwister(123))

samples |> viewer
```

Our goal is to interpolate the variable `Z` over the original domain.
Let's start by estimating the `EmpiricalVariogram` from the samples.
Because the distribution of values in the Walker Lake is skewed, the
default `:matheron` estimator of the variogram shows a high nugget
effect:

```{julia}
g = EmpiricalVariogram(samples, "Z", maxlag = 100.0)
```

```{julia}
varioplot(g)
```

A better alternative in this case is to use the robust `:cressie`
estimator:

```{julia}
g = EmpiricalVariogram(samples, "Z", maxlag = 100.0, estimator = :cressie)
```

```{julia}
varioplot(g)
```

After estimating the empirical variogram, the next step consists of fitting
a theoretical model. The behavior near the origin resembles a `SphericalVariogram`:

```{julia}
γ = GeoStatsFunctions.fit(SphericalVariogram, g)
```

```{julia}
varioplot(γ, maxlag = 100.0)
```

Now that we extracted the geospatial correlation from the samples, we can
use this information in `Kriging` interpolation. Instead of fitting all the
samples at once like it is done in the `Interpolate` transform, we will
fit the `Kriging` model with a maximum number of neighbors with the
`InterpolateNeighbors` transform:

```{julia}
interp = samples |> InterpolateNeighbors(img.geometry, Kriging(γ))
```

```{julia}
interp |> viewer
```

::: {.callout-note}

## Tip for all users

The `InterpolateNeighbors` is recommended in 3D applications with
hundreds of thousands of measurements and very large grids.

:::

::: {.callout-note}

## Tip for all users

The `InterpolateMissing` transform can be used to interpolate `missing`
values in a geotable using the same algorithm of `InterpolateNeighbors`.
Likewise, the `InterpolateNaN` can be used to interpolate `NaN` values.

:::

## Congratulations!

Congratulations on finishing **Part IV** of the book. The interpolation
models introduced here are simple, yet very useful. Before we start our
journey with real-world applications of the framework, let's review what
we learned:

- **Geospatial interpolation** can be achieved with the `Interpolate`
  and `InterpolateNeighbors` transforms, and geostatistical models such
  as `IDW` and `Kriging`.
- Models such as `Kriging` exploit **geospatial correlation** to improve
  interpolation results. We can extract this information from samples using
  a two-step procedure:
  1. Estimate the `EmpiricalVariogram` from the available samples
  2. Perform `fit` of theoretical model with result from previous step
- `Interpolate` and `InterpolateNeighbors` are examples of geostatistical
  transforms. They can be easily inserted in more advanced pipelines as
  discussed in **Part II**.

In the next chapters, we will use the framework that we learned with real
data to illustrate how advanced geospatial data science can be done with just
a few lines of code. Once a solution is written in terms of the high-level tools
covered in previous chapters, it is trivial to improve computational performance
in pure Julia.

::: {.callout-note}

Feature and performance requests are **very welcome**. We invite all users
of the framework to [submit issues](https://github.com/JuliaEarth/GeoStats.jl/issues)
and contribute with our open source software stack.

:::