# Geostatistics and rainfall with R

With the rainy winter season behind us, it’s time to summarize how much precipitation we got this year. The southern Negev desert enjoyed above average annual rainfall (in stark contrast to the rest of the country, where we saw only 60% of the multi-year average)

I wrote about creating isohyetal contours some years ago. This time I’ll try to improve my interpolation by using the powerful, open source R statistics program.

We begin with the total annual precipitation measurements at discrete locations – rain gauges. Using this point data we interpolate the rainfall everywhere in the region. The interpolation method I choose is kriging, which works with a variance distribution function, the variogram, that we prepare as part of the procedure. Furthermore, the geostatistics procedure called co-kriging allows us to apply a second, independant variable to the analysis to predict even better the interpolated values. Co-kriging is applicable only when we know that there is some correlation between the main, dependant variable and the auxiliary variable. Intuitively, we are aware that in many areas rainfall is stronger in the mountains, and less at lower elevations. So, I will examine elevation of the rain gauges as the independant auxiliary variable.

We begin by starting R, loading some required libraries, and pulling in the rain gauge data.

<br />
library(sp)<br />
library(rgdal)<br />
library(gstat)<br />
rain.data &lt;- &quot;raindata_ttl_2014.csv&quot;<br />

Now I load a border line, just for spatial reference, and make a simple bubble plot of the data.

<br />
# Put data into a Spatial Data Frame<br />
gauges.xy &lt;- gauges[[c(&quot;x&quot;,&quot;y&quot;)]]<br />
coordinates(gauges) &lt;- gauges.xy<br />
# Set the correct CRS for the gauges data<br />
itm_proj4 &lt;- proj4string(border_il)<br />
proj4string(gauges) &lt;- itm_proj4<br />
class(gauges)<br />
summary(gauges)<br />
# Have a look<br />
plot(gauges, pch=16, col='blue', cex=precip/20)<br />

Here’s what we get:

Now we need to determine the variogram. We start by plotting some candidate variogram models. The parameters that we need to set are “cutoff” and “width”. By default, the cutoff (how far from each prediction point to go to look for real points) is chosen as one third the diagonal of the whole analysis region. And the width (the variance between real points within the cutoff) is 1/15 of the cutoff. We can get an idea how these parameters influence the model variograms by plotting some possible combinations. Once we have a good looking variogram model, we use those in the gstat variogram() function to create the model variogram:

<br />
# Create a variogram<br />
# Display some possible models<br />
plot(gstat::variogram(precip ~ 1 , gauges, cutoff=20000, width=2000))<br />
plot(gstat::variogram(precip ~ 1 , gauges, cutoff=30000, width=4000))<br />
plot(gstat::variogram(precip ~ 1 , gauges, cutoff=70000, width=5000))<br />
# Once we have a visually good variogram, use those parameters:<br />
vg &lt;- gstat::variogram(precip ~ 1 , gauges, cutoff=70000)<br />

Now we have a model, and we want to let R to actually fit a variogram to the point pairs that are included in the variogram model to create a fitted curve. The following code prints the final fitted variogram values and gives us the result below:

<br />
vg.fit &lt;- gstat::fit.variogram(vg, vgm(700, &quot;Exp&quot;, 20000, 300))<br />
# Here are the fitted variogram parameters:<br />
print(vg.fit)<br />
plot(vg, vg.fit)<br />
model    psill    range<br />
1   Nug 355.0346     0.00<br />
2   Exp 421.1879 12672.45<br />

With that we are now prepared to run the kriging and produce a raster of predicted rainfall throughout the region.

<br />
# create a grid onto which we will interpolate:<br />
# first get the range in data<br />
x.range &lt;- as.integer(range(gauges@coords[,1]))<br />
y.range &lt;- as.integer(range(gauges@coords[,2]))<br />
# now expand to a grid with 100 meter spacing:<br />
grd &lt;- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=100), y=seq(from=y.range[1], to=y.range[2], by=100) )<br />
# convert to SpatialPixel class<br />
coordinates(grd) &lt;- ~ x+y<br />
gridded(grd) &lt;- TRUE<br />
proj4string(grd) &lt;- itm_proj4</p>
<p># perform ordinary kriging prediction:<br />
# make gstat object to hold the krige result:<br />
g &lt;- gstat(id=&quot;Precipitation&quot;, formula=precip ~ 1, data=gauges, model=vg.fit)<br />
precip.krige &lt;- predict(g, model=vg.fit, newdata=grd)<br />

The predict() function chooses ordinary kriging, based on the basic parameters we passed to the function. Now we can go ahead and plot the result. I load the RColorBrewer library to take advantage of the full set of color palettes available, and I use the spplot function from the sp package for plotting, since it easily adds contours (isohyetal lines) to the plot.

<br />
# Some parameters for plotting<br />
par(font.main=2, cex.main=1.5,cex.lab=0.4, cex.sub=1)<br />
# Use the ColorBrewer library for color ramps<br />
library(RColorBrewer)<br />
precip.pal &lt;- colorRampPalette(brewer.pal(7, name=&quot;Blues&quot;))<br />
# plot the krige interpolation<br />
spplot(precip.krige, zcol='Precipitation.pred', col.regions=precip.pal, contour=TRUE, col='black', pretty=TRUE, main=&quot;Interpolated Precipitation - 2014&quot;, sub=&quot;Ordinary Kriging&quot;, labels=TRUE)<br />

The plot of krige interpolated precipitation looks like this:

Now, as promised, let’s continue on to elevation data. We import a DEM and use the R over() function to attach an elevation value to each gauge. Then we check correlation between the annual precipitation and elevation. These are the steps:

<br />
class(elev)<br />
summary(elev)<br />
# Add elevation values to each gauge<br />
gauges.elev &lt;- over(gauges, elev)<br />
precip.elev  &lt;-  cbind(gauges, gauges.elev)<br />
coordinates(precip.elev) &lt;- gauges.xy<br />
proj4string(precip.elev) &lt;- itm_proj4<br />
# Check that we have the elevations in a spatial data frame<br />
names(precip.elev)<br />
# the column &quot;band1&quot; contains the elevation values<br />
# Now check for correlation between precipitation and elevation<br />
cor.test(precip.elev\$band1, precip.elev\$precip)<br />

Checking the output of the cor.test() function we see:

<br />
Pearson's product-moment correlation<br />
data:  precip.elev\$band1 and precip.elev\$precip<br />
t = 0.331, df = 96, p-value = 0.7414<br />
alternative hypothesis: true correlation is not equal to 0<br />
95 percent confidence interval:<br />
-0.1657725  0.2306346<br />
sample estimates:<br />
cor<br />
0.03375868<br />

Well, a correlation of 0.033 means essentially no correlation. Values of correlation can range from 1 (prefect correlation) to -1 (perfect inverse correlation). Values near to 0 indicate nearly no correlation. So for this last season the spread of the total annual rainfall had no connection to elevation across the Negev. We will settle with the above kriging map as a good estimate of the distribution of rainfall for 2014. End of story.

But, we don’t give up so easily. Let’s go back a few years to the winter of 2010. Anyone dealing with rainfall and flooding in the Negev desert remembers that season and the storms of January. We had once in 50 year floods in several wadis throughout the south due to very heavy rains in the Negev mountains. Let’s search thru the archives to pull out the annual rainfall data for that year, and try to correlate with elevation. I rerun all the above steps for the gauge data for 2010. Of course the variogram will be different, and I obtain the ordinary kriging result as below:

Now I load the same elevation dataset and as above I test correlation between the precipitation data and elevation values:

<br />
cor.test(precip.2010.elev\$band1, precip.2010.elev\$precip)<br />
Pearson's product-moment correlation</p>
<p>data:  precip.2010.elev\$band1 and precip.2010.elev\$precip<br />
t = 6.6778, df = 71, p-value = 4.552e-09<br />
alternative hypothesis: true correlation is not equal to 0<br />
95 percent confidence interval:<br />
0.4562349 0.7447525<br />
sample estimates:<br />
cor<br />
0.6211078<br />

Well, a correlation of 0.62 is nothing to write home about, but is does indicate some reasonable level of correlation. Additionally, the p-value is very low, meaning that there is a high level of confidence that the resulting correlation is accurate. So let’s forge ahead and use the elevation values to (hopefully) improve our precipitation estimation.

We first need to redo the gstat object to hold both variables:

<br />
# recreate g, the gstat object with a second variable<br />
rm(g)<br />
g &lt;- gstat(id=&quot;Precip&quot;, formula=precip ~ 1, data=precip.2010.elev, model=vgm(psill=2000, model=&quot;Exp&quot;, range=100000, nugget=50))<br />
g &lt;- gstat(g, id=&quot;Elevation&quot;, formula=band1 ~ 1, data=precip.2010.elev, model=vgm(psill=100000, model=&quot;Exp&quot;, range=50000, nugget=10))<br />
vg &lt;- gstat::variogram(g)<br />
# Make the multivariable variogram (Linear Model of Coregionalization)<br />
vg.fit &lt;- fit.lmc(vg, g, vgm(psill=2000, model=&quot;Exp&quot;, range=100000, nugget=50))<br />
# Graph the model and experimental variograms<br />
plot(vg, vg.fit)<br />

Note that we us the fit.lmc() function to create a fitted variogram with multiple variables, each variable with its own set of sill and range parameters. And we must add a non-zero nugget value to insure that the Linear Model of Coregionalization will be chosen.
The plot() shows us three variograms: one for elevation, one for precipitation, and a third for the covariance of the two variables together.

What’s left is to rerun the predict() function, and note that it chooses co-kriging:

<br />
# Now predict() should create a CoKriging interpolation<br />
precip.2010.cokrige &lt;- predict.gstat(vg.fit, newdata=grd)<br />
Linear Model of Coregionalization found. Good.<br />
[using ordinary cokriging]<br />

and here’s our new estimate of distribution of precipitation for 2010. Compared to the ordinary kriging map for 2010 above, we can identify an obvious change of rainfall distribution in the western mountain regions: