# Spatial distribution - lecture

- The problem
- Thiessen method [Voronoi diagram]
- Inverse distance method
- Isohyetal method
- How it is actually done

## Thiessen method [Voronoi diagram]

Brutsaert, Figure 3.11

How to compute the areas:

Average areal precipitation is a weighted sum:

$$ \langle P \rangle = \frac{\sum_i A_i P_i}{\sum_i A_i} $$A nice way to understand the Thiessen method is depicted in the gif below (from Wikipedia):

## Inverse distance method

Brutsaert, Figure 3.12

The precipitation for square 17 is

$$ P_{17} = \displaystyle\frac {\displaystyle\sum_\text{$i$ = all stations}\frac{P_i}{d_{i,17}^2}} {\displaystyle\sum_\text{$i$ = all stations}\frac{1}{d_{i,17}^2}} $$The average precipitation for the whole watershed is the weighted average of all squares, where the weight is their area:

$$ \langle P \rangle = \displaystyle\frac {\displaystyle\sum_\text{$j$ = all squares} A_j P_j} {\displaystyle\sum_\text{$j$ = all squares} A_j} $$Brutsaert, page 93:

Dean and Snyder (1977) found that the exponent (for the distance $d^{-b}$) b = 2 yielded the best results in the Piedmont region of the southeastern United States, whereas Simanton and Osborn (1980) concluded from measurements in Arizona that b can range between 1 and 3 without significantly affecting the results.

## How it is actually done

Most often, Geographic Information System (GIS) software is used to analyze spatial data. Two of the most used programs are ArcGIS (proprietary) and QGIS (free).

A good discussion of the different methods can be found on Manuel Gimond's website, *Intro to GIS and Spatial Analysis*.

**Attention**, Don't mix precision with accuracy. There are many ways of interpolating, just because a result seems detailed, it does not imply that it is accurate! See below three interpolation methods.

Below you can find a simple Python code that exemplifies some of the methods, producing the following figure:

```
%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import griddata
from scipy.spatial import Voronoi, voronoi_plot_2d, ConvexHull
fig, ax = plt.subplots(1, 3, figsize=(10,7))
fig.subplots_adjust(left=0.0, right=1.0, top=0.96, bottom=0.05,
hspace=0.02, wspace=0.02)
N = 6
PI = '3141592653589793'
points = np.random.rand(N, 2)
points = np.vstack([points,[0,0], [0,1], [1,0], [1,1]])
values = np.array([int(x) for x in list(PI)])[:(N+4)]
# values = np.array([3, 1, 4, 1, 5, 9, 2, 6, 5, 3])
grid_x, grid_y = np.mgrid[0:1:100j, 0:1:200j]
grid_z_nearest = griddata(points, values, (grid_x, grid_y), method='nearest')
grid_z_cubic = griddata(points, values, (grid_x, grid_y), method='cubic')
ax[0].plot(points[:,0], points[:,1], 'o', ms=3, markerfacecolor="red", markeredgecolor="red")
ax[0].set_aspect('equal', 'box')
ax[0].set(xlim=[0,1], ylim=[0,1])
ax[0].set_title("the stations")
for i, v in enumerate(values):
ax[0].text(points[i,0], points[i,1], str(v))
ax[1].imshow(grid_z_nearest.T, extent=(0,1,0,1), origin='lower')
ax[1].plot(points[:,0], points[:,1], 'o', ms=3, markerfacecolor="red", markeredgecolor="red")
vor = Voronoi(points)
voronoi_plot_2d(vor, show_vertices=False, line_colors='cyan',
line_width=3, line_alpha=1, point_size=0, ax=ax[1])
ax[1].set_title("Thiessen Method")
ax[2].plot(points[:,0], points[:,1], 'o', ms=3, markerfacecolor="red", markeredgecolor="red")
nlines = int((values.max()-values.min()+1)/2)
ax[2].contourf(grid_x, grid_y, grid_z_cubic, nlines)
cont = ax[2].contour(grid_x, grid_y, grid_z_cubic, nlines, colors="black")
ax[2].clabel(cont, inline=1, colors='white', fmt='%.0f')
ax[2].set_title("Isohyetal Method")
for i, a in enumerate(ax):
a.set(xlim=[-0.2,1.2], ylim=[-0.2,1.2])
a.axis('off')
a.set_aspect('equal', 'box')
fig.savefig("spatial-distribution.png", dpi=500)
```