# 17 Spatial distribution - lecture

## 17.1 The problem

Let’s say we want to calculate the average rainfall on a watershed, and we have data available for 7 stations, as shown in the figure below [Dingman, figure 4.26]:

There are a number of methods for calculating the average precipitation.

## 17.2 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):

## 17.3 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.

## 17.4 Isohyetal method

Brutsaert, Figure 3.12

The same equation of the Thiessen method can be used:

\[ \langle P \rangle = \frac{\sum_i A_i P_i}{\sum_i A_i} \]

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