\chapter{Spatial Models}
\label{chap:spatialmodels}
\SweaveOpts{keep.source=TRUE, pdf=FALSE, prefix.string=Chap09, grdevice=tikz.Swd}
<<options and clear workspace, echo=FALSE,eval=TRUE>>=
options(digits=3, show.signif.stars=FALSE, width=53)
rm(list=ls())
require(tikzDevice)
source("../SweaveTikZ.R")
@
%data files: load("best.use.RData"); read.table("sstJuly2005.txt", T); read.table("FayRain.txt", T)
%packages: gstat, maptools, sp, rgdal, maps, colorRamps, spdep, spgwr
%source code: NONE
%third party: NONE

\begin{quote} ``Design is a question of substance, not just form.''
\end{quote}
\indent ---Adriano Olivetti\\

In Chapter~\ref{chap:frequencymodels} annual counts are used to create rate models and in Chapter~\ref{chap:intensitymodels} lifetime maximum winds are used to create intensity models.  Here we show you how to use cyclone track data together with climate field data to create spatial models.

Spatial models make use of location information in data.  Geographic coordinates locate the cyclone's center on the surface of the Earth and wind speed provides an attribute.  Spatial models make use of cyclone location separate from cyclone attributes.  Given a common spatial support, these models accommodate climate data including indices (e.g., North Atlantic Oscillation) and fields (e.g., sea-surface temperature).

\section{Track Hexagons}

Here we show you how to create a spatial framework for combining cyclone data with spatially continuous climate data.  The method tessellates the basin with hexagons and populates them with local cyclone and climate information \citep{ElsnerHodgesJagger2012}.

\subsection{Spatial points data frame}

In Chapter~\ref{chap:graphsandmaps} you learned how to create a spatial data frame using functions from the {\bf sp} package \citep{BivandEtAl2008}.  Let's review.  Here you are interested in wind speeds along the entire track for all tropical storms and hurricanes during the 2005 North Atlantic season.  You begin by creating a data frame from the {\it best.use.RData} file where you subset on year and wind speed and convert the speed to meters per second.
<<load data, echo=TRUE, eval=TRUE>>=
load("best.use.RData")
W.df = subset(best.use, Yr==2005 & WmaxS >= 34
  & Type=="*")
W.df$WmaxS = W.df$WmaxS * .5144
@
The asterisk on for \verb@Type@ indicates tropical cyclone as opposed to a tropical wave or extratropical cyclone.  The number of rows in your data frame is the total number of cyclone hours (\Sexpr{nrow(W.df)}) and you save this by typing
<<number of cyclones, echo=TRUE, eval=TRUE>>=
ch = nrow(W.df)
@

Next, assign the \verb@lon@ and \verb@lat@ columns of the data frame as spatial coordinates using the \verb@coordinates@ function ({\bf sp}). Make a copy of your data frame keeping only the spatial coordinates and the wind speed columns.
<<make spatial points data frame, echo=TRUE, eval=TRUE>>=
require(sp)
W.sdf = W.df[c("lon", "lat", "WmaxS")]
coordinates(W.sdf) = c("lon", "lat")
str(W.sdf, max.level=2, strict.width="cut")
@
The result is a spatial points data frame with five slots.  The data slot is a data frame and contains the attributes (here only wind speed).  The coordinate (\verb@coord@) slot contains the longitude and latitude columns from the original data frame and the coordinate numbers (here two spatial dimensions) are given in the \verb@coords.nrs@ slot.  Note that these two columns are not in the data slot.  The bounding box (\verb@bbox@) slot is a matrix containing the maximal extent of the hourly positions as defined by the lower left and upper right longitude/latitude coordinates.

A summary of the information in your spatial points data frame is obtained by typing
<<summary of spatial points data frame, echo=TRUE, eval=TRUE>>=
summary(W.sdf)
@
There are \Sexpr{dim(W.sdf)[1]} hourly cyclone hours.  The projection (\verb@proj4string@) slot contains an \verb@NA@ character indicating that it has not yet been specified.

You give your spatial data frame a coordinate reference system using the \verb@CRS@ function.   Here you specify a geographic reference (intersections of parallels and meridians) as a character string in an object called \verb@ll_crs@, then use the \verb@proj4string@ function to generate a CRS object and assign it to your spatial data frame.
<<projection, echo=TRUE, eval=TRUE>>=
ll = "+proj=longlat +datum=WGS84"
proj4string(W.sdf) = CRS(ll)
@
Check this slot by typing
<<check slot, echo=TRUE, eval=TRUE>>=
slot(W.sdf, "proj4string")
@

Next you transform the geographic CRS into a Lambert conformal conic (LCC) planar projection using the parallels 30 and 60$^\circ$N and a center longitude of 60$^\circ$W.  First save the reference system as a CRS object then use the \verb@spTransform@ function from the {\bf rgdal} package.
<<new projection, echo=TRUE, eval=TRUE>>=
lcc = "+proj=lcc +lat_1=60 +lat_2=30 +lon_0=-60"
require(rgdal)
W.sdf = spTransform(W.sdf, CRS(lcc))
bbox(W.sdf)
@
The coordinates are now planar rather than geographical although the coordinate names remain the same (\verb@lon@ and \verb@lat@).  Transformation to a plane makes it easy to perform geometric calculations.

\subsection{Hexagon tessellation}

Next you construct a hexagon tessellation of the cyclone tracks.  This is done by situating the points representing the center of each hexagon and then drawing the hexagon boundaries.  First you create a set of equally-spaced points in a hexagon grid configuration within the bounding box defined by your spatial points data frame using the \verb@spsample@ function.  Here you specify the number of hexagons, but the actual number will vary depending on the bounding box.  You expand the bounding box by multiplying each coordinate value by 20\% and use an offset vector to fix the position of the hexagons over the tracks.
<<hexagon points, echo=TRUE, eval=TRUE>>=
hpt = spsample(W.sdf, type="hexagonal", n=250,
  bb=bbox(W.sdf) * 1.2, offset=c(1, -1))
@
The expansion and offset values require a bit of trail-and-error.  The methods used in spatial sampling function assume that the geometry has planar coordinates, so your spatial data object should not have geographic coordinates.

Next you call the function to convert the points to hexagons.  The more points the smaller the area of each hexagon.
<<points to polygons, echo=TRUE, eval=TRUE>>=
hpg = HexPoints2SpatialPolygons(hpt)
@
The number of polygons generated and the area of each polygon is obtained by typing
<<number and area of polygons, echo=TRUE, eval=TRUE>>=
np = length(hpg@polygons)
area = hpg@polygons[[1]]@area
np; area
@
This results in \Sexpr{np} equal-area hexagons.  The length unit of the LCC projection is meters.  To convert the area from square meters to square kilometers, multiply by 10$^{-6}$.  Thus the area of each hexagon is approximately \Sexpr{round(area*1.0e-6,0)}~km$^2$.

\subsection{Overlays}

With your hexagons and cyclone locations having the same projection, next you obtain the maximum intensity and number of observations per grid.  The function \verb@over@ combines points (or grids) and polygons by performing point-in-polygon operations on all point-polygon combinations.  First you obtain a vector containing the hexagon identification number for each hourly cyclone observation by typing
<<hexagon ids for each cyclone hour, echo=TRUE, eval=TRUE>>=
hexid = over(x=W.sdf, y=hpg)
@
The length of the vector is the number of hourly observations.  Not all hexagons have cyclones, so you subset your spatial polygons object keeping only those that do.
<<subset spatial polygons, echo=TRUE, eval=TRUE>>=
hpg = hpg[unique(hexid)]
@
Next you create a data frame with a single column listing the maximum wind speed over the set of all cyclone observations in each hexagon.
<<highest per grid intensity, echo=TRUE, eval=TRUE>>=
int = over(x=hpg, y=W.sdf, fn=max)
colnames(int) = c("WmaxS")
head(int)
@
The result is a data frame with a single column representing the maximum wind speed value over all cyclone observations in each hexagon.  The row names are the hexagon numbers prefixed with {\tt ID}.  Finally, you combine the spatial polygons with the maximum per hexagon wind speeds into a spatial polygons data frame.
<<hspdf creation, echo=TRUE, eval=TRUE>>=
hspdf = SpatialPolygonsDataFrame(hpg, int,
  match.ID = TRUE)
@

You plot the hourly storm locations together with an overlay of your hexagons by typing
<<plot polygons, echo=TRUE, eval=FALSE>>=
plot(hspdf)
plot(W.sdf, pch=20, cex=.3, add=TRUE)
@
<<alternative specification of hexagon grid, echo=FALSE, eval=FALSE>>=
W.sdf = W.df[c("lon", "lat", "WmaxS")]
coordinates(W.sdf) = c("lon", "lat")
ll = "+proj=longlat +datum=WGS84"
proj4string(W.sdf) = CRS(ll)
bbox(W.sdf)
lcc = "+proj=lcc +lat_1=60 +lat_2=30 +lon_0=-60"
W.sdf = spTransform(W.sdf, CRS(lcc))
bb = matrix(c(-100, -5, -100, 45, -10, 45, -10, -5, -100, -5), byrow=TRUE, ncol=2)
bb.sp = SpatialPolygons(list(Polygons(list(Polygon(bb)), "1")), proj4string=CRS(ll))
bb.sp = spTransform(bb.sp, CRS(lcc))
hpt = spsample(bb.sp, type="hexagonal", n=250, offset=c(.5, -.5))
hpg = HexPoints2SpatialPolygons(hpt)
hpgk = hpg[sort(unique(over(x=W.sdf, y=hpg)))]
int = over(x=hpgk, y=W.sdf, fn=max)
colnames(int) = c("WmaxS")
hspdf = SpatialPolygonsDataFrame(hpg, int, match.ID = TRUE)
plot(hspdf)
plot(W.sdf, pch=20, cex=.3, add=TRUE)
@
Some hexagons contain many cyclone observations while others contain only one or two.  This difference might be important in modeling intensity so you add total cyclone hours as a second attribute.  First, replace the data in the spatial points data frame with an index of ones.
<<replace attribute data, echo=TRUE, eval=TRUE>>=
W.sdf@data = data.frame(num=rep(1, ch))
@
Then perform an overlay of the hexagons on the cyclone locations with the function argument (\verb@fn@) set to \verb@sum@).  Check that the sum of the counts over all grids equals the total number of cyclone hours.
<<over and sum, echo=TRUE, eval=TRUE>>=
co = over(x=hspdf, y=W.sdf, fn=sum)
sum(co) == ch
@
Finally, add the counts to the spatial polygons data frame.  Here the rows of \verb@co@ correspond to those in the data slot of \verb@hspdf@ so there is no need to match the polygons IDs.
<<add the counts, echo=TRUE, eval=TRUE>>=
hspdf$count = co[, 1]
head(slot(hspdf, "data"))
@
Here you list the first six rows of the data frame that correspond to the first six hexagon grids.

\subsection{Maps}

You now have a spatial polygons data frame with each polygon as an equal-area hexagon on a LCC projection and two attributes (maximum wind speed and cyclone count) in the data slot.  Choropleth maps of cyclone attributes are created using the \verb@spplot@ method introduced in Chapter~\ref{chap:graphsandmaps}.

Before mapping you create a couple of lists that are used by the \verb@sp.layout@ argument.  One specifies the set of hourly locations from the spatial points data frame and assigns values to plot arguments.
<<points list, echo=TRUE, eval=TRUE>>=
l1 = list("sp.points", W.sdf, pch=20, col="gray",
  cex=.3)
@
Another specifies the coastlines as a spatial lines object.  Some work is needed up front.

First require the {\bf maps} and {\bf maptools} packages.  The first package contains the \verb@map@ function to generate an object of country borders and the second contains the \verb@map2SpatialLines@ conversion function.  The conversion is made to geographic coordinates, which are then transformed to the LCC projection of your above spatial objects.  You set geographic coordinate limits on the map domain to limit the amount of conversion and projection calculations.
<<get map borders, echo=TRUE, eval=TRUE>>=
require(maps)
require(maptools)
cl = map("world", xlim=c(-120, 20),
  ylim=c(-10, 70), plot=FALSE)
clp = map2SpatialLines(cl, proj4string=CRS(ll))
clp = spTransform(clp, CRS(lcc))
l2 = list("sp.lines", clp, col="gray")
@

Depicting attribute levels on a map is done using a color ramp.  A color ramp creates a character vector of color hex codes.  The number of colors is specified by the argument in the color ramp function.  Numerous color ramp functions are available in the {\bf colorRamps} package \citep{Keitt2009}.  Acquire the package and assign 20 colors to the vector \verb@cr@ using the blue to yellow color ramp.
<<get color ramps, echo=TRUE, eval=TRUE>>=
require(colorRamps)
cr = blue2yellow(20)
@
If the number of colors is less than the number of levels the colors get recycled.  A map of the cyclone frequency is made by typing
<<frequency map, fig=FALSE, echo=TRUE, eval=FALSE>>=
spplot(hspdf, "count", col="white",
   col.regions=blue2yellow(20),
   sp.layout=list(l1, l2),
   colorkey=list(space="bottom"),
   sub="Cyclone Hours")
@
Similarly a map of the highest cyclone intensity is made by typing
<<intensity map, fig=FALSE, echo=TRUE,eval=FALSE>>=
spplot(hspdf, "WmaxS", col="white",
   col.regions=blue2red(20),
   sp.layout=list(l2),
   colorkey=list(space="bottom"),
   sub="Highest Intensity (m/s)")
@
The results are shown in Fig.~\ref{fig:freqint}.
\begin{figure}
\centering
<<freqintensitymaps, fig=TRUE, echo=FALSE, eval=TRUE, width=4.6, height=6>>=
require(grid)
p1 = spplot(hspdf, "count", col="white", col.regions=blue2yellow(20),
  sp.layout=list(l1, l2), at=seq(0, 140, 20),
  colorkey=list(space="bottom", labels=list(labels=paste(seq(0, 140, 20)), cex=.8)),
  sub=list("Cyclone hours", cex=.8, font=1))
p2 = spplot(hspdf, "WmaxS", col="white", col.regions=blue2red(20),
  sp.layout=list(l2), at=seq(10, 90, 10),
  colorkey=list(space="bottom", labels=list(labels=paste(seq(10, 90, 10)), cex=.8)),
  sub=list("Highest intensity [m s$^{-1}$]", cex=.8, font=1))
p1 = update(p1, main=textGrob("a", x=unit(.05, "npc")))
p2 = update(p2, main=textGrob("b", x=unit(.05, "npc")))
print(p1, split=c(1, 1, 1, 2), more=TRUE)
print(p2, split=c(1, 2, 1, 2), more=FALSE)
@
\vspace{0cm}
\caption{Cyclone frequency and intensity.  (a) Hours and (b) highest intensity.}
\label{fig:freqint}
\end{figure}
A
reas along the southeastern coastline of the United States have the greatest number of cyclone hours while the region from the western Caribbean into ^the Gulf of Mexico has the highest intensities.

\section{SST Data}

^^^^W^e^ retur^^^^^^^^^n to these data here.  Input ^t^hem by typing
<<read SST data, echo=TRUE, eval=TRUE>>=
sst = read.table("JulySST2005.txt", header=TRUE)
head(sst)
@
The data are listed in the column labeled {\tt SST}.  You treat these values as point data because they have longitudes and latitudes although they are regional averages.

At 25$^\circ$N latitude the 2$^\circ$ spacing of the SST values covers^ an area of approximately the s^^ame size as the hexagon grids used in the previous ^^^^^^^^patial data frame.  Recall you need to first assign a projection string to the {\tt proj4string} slot.
<<spatial point data frame and coordinate transform, echo=TRUE, eval=TRUE>>=
coordinates(sst) = c("lon", "lat")
proj4string(sst) = CRS(ll)
sst = spTransform(sst, CRS(lcc))
@
To examine the SST data as spatial points, type
<<plot SST, echo=TRUE, eval=FALSE>>=
spplot(sst, "SST", col.regions=rev(heat.colors(20)))
@
This produces a plot of the spatial points data frame, where the SST attribute is specified as a character string in the second argument.

Your interest is the average SST within each hexagon.  Again, use \verb@over@  to combine your SST values with the track polygons setting the argument \verb@fn@ as \verb@mean@.

<<overlay sst, echo=TRUE, eval=TRUE>>=
ssta = over(x=hspdf, y=sst, fn=mean)
@
The result is a data frame with a single column representing the average over all SST values in each hexagon.  The row names are the hexagon numbers prefixed with {\tt ID}, however, hexagons completely over land have missing SST values.

Next, you add the average SSTs as an attribute to spatial polygon data frame and then remove hexagons with missing SST values.
<<add sst to spatial polygons, echo=TRUE, eval=TRUE>>=]
hspdf$sst = ssta$SST
hspdf = hspdf[!is.na(hspdf$sst),]
str(slot(hspdf, "data"))
@

Finally you create the July SST map corresponding to where the cyclones occurred during the season by typing
<<SST map, fig=FALSE, echo=TRUE, eval=FALSE>>=
spplot(hspdf, "sst", col="white",
   col.regions=blue2red(20),
   sp.layout=list(l1, l2),
   colorkey=list(space="bottom"),
   sub="Sea Surface Temperature (C)")
@
Results are shown in Fig.~\ref{fig:sst}.  Ocean temperatures exceed 26$^\circ$C over a wide swath of the North Atlantic extending into the Caribbean Sea and Gulf of Mexico.  Coldest waters are noted off the coast of Nova Scotia and Newfoundland.
\begin{figure}
\centering
<<SSTmapfigure, fig=TRUE, echo=FALSE, eval=TRUE, width=4.6, height=3>>=
print(spplot(hspdf, "sst", col="white", col.regions=blue2red(20),
  sp.layout=list(l1, l2), at=seq(14, 30, 2),
  colorkey=list(space="bottom", labels=list(cex=.8)),
  sub=list("Sea surface temperature [$^\\circ$C]", cex=.8, font=1)))
@
\vspace{0cm}
\caption{Sea-surface temperature in the cyclone hexagons.}
\label{fig:sst}
\end{figure}

\section{SST and Intensity}

Analyzing and modeling your SST and cyclone data is easy with your spatial polygons data frame.   For instance, the average maximum wind speed in regions where the SST exceeds 28$^\circ$C is obtained by typing
<<conditional intensity, echo=TRUE, eval=TRUE>>=
mean(hspdf$WmaxS[hspdf$sst > 28])
@
Here you treat your spatial data frame as you would a regular data frame.  Continuing you create side-by-side box plots of wind speed conditional on the collocated values of SST (see Chapter~\ref{chap:graphsandmaps}) by typing
<<box plots, echo=TRUE, eval=FALSE>>=
boxplot(hspdf$WmaxS~hspdf$sst > 28)
@

Importantly, the spatial information allows you to analyze relationships on a map.  Figure~\ref{fig:grouped} shows the hexagons colored by groups defined by a two-way table of cyclone intensity and SST using above and below median values.   The median SST and storm intensity values calculated from your data are \Sexpr{round(median(hspdf$sst),1)}$^\circ$C and \Sexpr{round(median(hspdf$WmaxS),1)}~m~s$^{-1}$, respectively.
\begin{figure}
\centering
<<groupedmapfigure, fig=TRUE, echo=FALSE, eval=TRUE, width=4.1, height=2.8>>=
mT = median(hspdf$sst); mI = median(hspdf$WmaxS)
hspdf$group = ifelse(hspdf$sst >= mT & hspdf$WmaxS >= mI, 4,
   ifelse(hspdf$sst >= mT & hspdf$WmaxS < mI, 3,
   ifelse(hspdf$sst <= mT & hspdf$WmaxS >= mI, 2,
   ifelse(hspdf$sst <= mT & hspdf$WmaxS < mI, 1, 0))))
mycolors = c("blue", "cyan", "magenta", "red")
txt = list(c("Low SST,   Low intensity", "Low SST,  High intensity", "High SST,  Low intensity", "High SST, High intensity"))
print(spplot(hspdf, "group", col="white", col.regions=mycolors,
  sp.layout=list(l1, l2), colorkey=F, key=list(space="bottom",
  rect=list(col=mycolors), text=txt, cex=.7)))
@
\vspace{0cm}
\caption{SST and cyclone intensity.  Groups are based on median values.}
\label{fig:grouped}
\end{figure}
Red hexagons indicate regions of high intensity and relatively high ocean temperature and blue hexagons indicate regions of low intensity and relatively low ocean temperature.  More interesting are regions of mismatch.  Magenta hexagons show low intensity coupled with relatively high ocean temperature indicating `under-performing' cyclones (cyclones weaker than the thermodynamic potential of their environment).  In contrast, cyan hexagons show high intensity coupled with relatively low temperature indicating `over-performing' cyclones (cyclones stronger than the thermodynamic potential of their environment).

Maps like these provide insight into the relationship between hurricanes and climate not accessible with basin-level analyzes.  Next we show how you how to model these spatial data.  We begin with a look at spatial correlation and then use a spatial regression model to quantify the regional change in intensity as a function of SST.

\section{Spatial Autocorrelation}

Values in neighboring hexagons will tend to be more similar than values in hexagons farther away.  Spatial autocorrelation quantifies the degree of similarity across geographic space.  It is a bit more complex than the familiar temporal autocorrelation because space is multi-dimensional and multi-directional.

\subsection{Moran's I}

A measure of spatial autocorrelation is Moran's $I$ \citep{Moran1950} defined as
\begin{equation}
I = \frac{m}{s} \frac{y^TWy}{y^Ty}
\label{eq:moransI}
\end{equation}
where $m$ is the number of hexagons, $y$ is the vector of values within each hexagon (e.g., cyclone intensities) where the values are deviations from the overall mean,  $W$ is a weights matrix, $s$ is the sum over all the weights, and the subscript $T$ indicates the transpose operator.

Values of Moran's $I$ range from $-$1 to $+$1 with a value of zero indicating a pattern with no spatial autocorrelation.  Although not widely used in climate studies, \cite{deBeursHenebry2008} use it to identify spatially coherent eco-regions and biomes related to the North Atlantic oscillation.

To compute Moran's $I$ you need a weights matrix.  The weights matrix is square with the number of rows and columns equal to the number of hexagons.  Here the weight in row $i$, column $j$ of the matrix is assigned a non-zero value if hexagon $i$ is contiguous with hexagon $j$, otherwise it is assigned a zero.  The {\bf spdep} package \citep{Bivand2011} has functions for creating weights based on contiguity (and distance) neighbors.  The process requires two steps.

First, you use the \verb@poly2nb@ function on your spatial polygons data frame to create a contiguity-based neighborhood list object.
<<neighbor list, echo=TRUE, eval=TRUE, results=hide>>=
require(spdep)
hexnb = poly2nb(hspdf)
@
A summary of the neighborhood list is obtained by typing
<<summary of neighborhood list, echo=TRUE, eval=TRUE>>=
hexnb
@
The list is ordered by hexagon number starting with the southwest-most hexagon.  It has \Sexpr{length(hexnb[[1]])} neighbors; hexagon numbers \Sexpr{hexnb[[1]][1]} and \Sexpr{hexnb[[1]][2]}.  Hexagon numbers increase to the east and north.  A hexagon has at most six contiguous neighbors.  Hexagons at the borders have fewer neighbors.   A graph of the hexagon connectivity, here defined by the first-order contiguity, is made by typing
<<neighbor plot, echo=TRUE, eval=FALSE>>=
plot(hexnb, coordinates(hspdf))
plot(hspdf, add=TRUE)
@
A summary method applied to the neighborhood list (\verb@summary(hexnb)@) reveals the average number of neighbors and the distribution of connectivity among the hexagons.

You turn the neighborhood list object into a \verb@listw@ object using the \verb@nb2listw@ function that duplicates the neighborhood list and adds the weights.  The \verb@style@ argument determines the weighting scheme.  With the argument value set to \verb@W@, the weights are the inverse of the number of neighbors.  For instance, the six neighbors of a fully connected hexagon each get a weight of 1/6.
<<weight matrix, echo=TRUE, eval=TRUE, results=hide>>=
wts = nb2listw(hexnb, style="W")
summary(wts)
@

Now you are ready to compute the value of Moran's $I$.  This is done using the \verb@moran@ function.  The first argument is the variable of interest followed by the name of the \verb@listw@ object.  Also needed is the number of hexagons and the global sum of the weights, which is obtained using the \verb@Szero@ function.
<<morans I, echo=TRUE, eval=TRUE>>=
m = length(hspdf$WmaxS)
s = Szero(wts)
moran(hspdf$WmaxS, wts, n=m, S0=s)
@
The function returns Moran's $I$ and the sample kurtosis.\footnote{Kurtosis is a measure of the peakedness of the distribution.  A normal distribution has a kurtosis of 3.}  The value of \Sexpr{round(moran(hspdf$WmaxS,wts,n=m,S0=s)$I,2)} indicates fairly high spatial autocorrelation in cyclone intensity.  This is expected since the strength of a cyclone as it moves across the hexagons often does not vary by much and because stronger cyclones tend to occur at lower latitudes.

\subsection{Spatial lag variable}

Insight into the Moran's $I$ statistic is obtained by noting that it is equivalent to the slope coefficient in a regression model of $Wy$ on $y$, where $Wy$ is called the spatial lag variable (see Eq.~\ref{eq:moransI}).  Let \verb@y@ be the set of cyclone intensities in each hexagon, then you create a spatial lag intensity variable using the \verb@lag.listw@ function.
<<spatial lag intensities, echo=TRUE, eval=TRUE>>=
y = hspdf$WmaxS
Wy = lag.listw(wts, y)
@
Thus for each cyclone intensity value in the vector object {\tt y} there is a corresponding value in the vector object {\tt Wy} representing the mean intensity over the neighboring hexagons.  The neighbor average does not include the value in {\tt y}, so for a completely connected hexagon the average is taken over the adjoining six neighboring values.

A scatter plot of the neighborhood average intensity versus the intensity in each hexagon shows the spatial autocorrelation relationship.  The slope of a least-squares regression line through the points is the value of Moran's $I$.  Use the code below to create the scatter plot shown in Fig.~\ref{fig:moranscatterplot}.
<<morans scatter plot, echo=TRUE, eval=FALSE>>=
par(las=1, pty="s")
plot(y, Wy, pch=20, xlab="Intensity (m/s)",
  ylab="Neighborhood Avg Intensity (m/s)")
abline(lm(Wy ~ y), lwd=2, col="red")
@
\begin{figure}
\centering
<<moranscatterplot, fig=TRUE, echo=FALSE, eval=TRUE, width=3.4, height=3.4>>=
par(las=1, pty="s", mgp=c(2, .4, 0), tcl=-.3)
plot(y, Wy, pch=20,
  xlab="Intensity [m s$^{-1}$]",
  ylab="Neighborhood avg intensity [m s$^{-1}$]")
grid()
abline(lm(Wy ~ y), lwd=2, col="red")
points(y, Wy, pch=20)
@
\vspace{-.7cm}
\caption{Moran's scatter plot for cyclone intensity.}
\label{fig:moranscatterplot}
\end{figure}

The scatter plot (Moran's scatter plot) suggests the possibility of a significant amount of spatial autocorrelation.  High intensity hexagons tend to be surrounded, on average, by other hexagons with high intensity and vice versa as evidenced by the upward slope of the regression line.

\subsection{Statistical significance}

The expected value of Moran's $I$ under the hypothesis of no spatial autocorrelation is
\begin{equation}
E(I) = \frac{-1}{n-1}
\end{equation}
where $n$ is the number of hexagons.  This allows you to test the significance of your sample Moran's $I$.  The test is available in the \verb@moran.test@ function.  The first argument is the vector of intensities and the second is the spatial weights matrix in the weights list form.
<<moran test, echo=TRUE, eval=TRUE>>=
moran.test(y, wts)
@
The output shows the standard deviate computed by taking the difference between the estimated $I$ and it's expected value under the null hypothesis of no autocorrelation.  The difference is divided by the square root of the difference between the variance of $I$ and the square of the mean of $I$.  The $p$-value is the chance of observing a standard deviate this large or larger assuming there is no spatial autocorrelation.  The $p$-value is extremely small leading you to conclude there is significant autocorrelation.  The output also includes the value of Moran's $I$ along with its expected value and variance.

By default the variance of $I$ is computed by randomizing the intensities across the hexagon tiling.  If the values are normally distributed, a version of the test that has more statistical power (higher probability that the test will reject the null hypothesis of no autocorrelation when there is, in fact, autocorrelation) by adding the argument \verb@random=FALSE@.  Moran's $I$ and the corresponding significance test are sensitive to the definition of neighbors and to the neighborhood weights, so conclusions should be stated as conditional on your definition of neighborhoods.

\section{Spatial Regression Models}

Models that take advantage of spatial autocorrelation are called spatial regression models.  If significant autocorrelation exists, spatial regression models have parameters that are more stable and statistical tests that are more reliable than non-spatial alternatives.  For instance, confidence intervals on a regression slope will have the proper coverage probabilities and prediction errors will be smaller.  Autocorrelation is included in a regression model by adding a spatial lag dependent variable or by including a spatially correlated error term \citep{AnselinEtAl2006}.

Spatial autocorrelation can also enter a model by allowing the relationship between the response and the explanatory variable to vary across the domain.  This is called geographically-weighted regression \citep{BrunsdonEtAl1998, FotheringhamEtAl2000}.  Geographically-weighted regression (GWR) allows you to see where an explanatory variable contributes strongly to the relationship and where it contributes weakly.  It is similar to a local linear regression.

For example, you compare a standard linear regression model of intensity on SST with a local linear regression of the same relationship by typing
<<local linear regression, echo=TRUE, eval=FALSE>>=
x = hspdf$sst
par(las=1, pty="s")
plot(x, y, pch=20, xlab="SST (C)",
  ylab="Intensity (m/s)")
abline(lm(y ~ x), lwd=2)
lines(loess.smooth(x, y, span=.85), lwd=2, col="red")
@
\begin{figure}
\centering
<<locallinear, fig=TRUE, eval=TRUE, echo=FALSE, width=3.4, height=3.4>>=
x = hspdf$sst
par(las=1, pty="s", mgp=c(2, .4, 0), tcl=-.3)
plot(x, y, pch=20,
   xlab="SST [$^\\circ$C]",
   ylab="Intensity [m s$^{-1}$]")
grid()
abline(lm(y ~ x), lwd=2)
lines(loess.smooth(x, y, span=.85), lwd=2,
  col="red")
points(x, y, pch=20)
legend("topleft", legend=c("Linear regression", "Local linear regression"),
  lwd=2, col=c(1, 2), cex=.8, bg="white")
@
\vspace{-.7cm}
\caption{Linear and local linear regressions of cyclone intensity on SST.}
\label{fig:locallinear}
\end{figure}
With the local regression the relationship between intensity and SST changes for different values of SST (Fig.~\ref{fig:locallinear}).  For lower values of SST the relationship has a gentle slope and with higher values of SST the relationship has a steeper slope.  In contrast, the `global' linear regression results in a single moderate slope across all values of SST.

Fitting is done locally in the domain of SST. That is for a point $s$ along the SST axis, the fit is made using points in a neighborhood of $s$ weighted by their distance to $s$. The size of the neighborhood is controlled by the \verb@span@ argument.  With \verb@span=.85@ the neighborhood includes 85\% of the SST values.

With GWR this localization is done in the domain of geographic space.  For example, a regression of cyclone intensity on SST is performed using the paired intensity and SST values at each of the \Sexpr{length(y)} hexagons.  For each hexagon the weight associated with a paired value in another hexagon is inversely proportional to physical distance between the two hexagons.  In this way the relationship between intensity and SST is localized.

\subsection{Linear regression}

The standard regression model consists of a vector $y$ of response values and a matrix $X$ containing the set of explanatory variables plus a row vector of one's.  The relationship is modeled as
\begin{equation}
y = X\beta+ \varepsilon
\end{equation}
where $\beta$ is a vector of regression coefficients and $\varepsilon \sim N(0, \sigma^2I)$ is a vector of independent and identically distributed residuals with variance $\sigma^2$.  The maximum likelihood estimate of $\beta$ is given by
\begin{equation}
\hat \beta = (X^TX)^{-1}X^T y.
\end{equation}

You begin with a linear regression of cyclone intensity on SST across the set of hexagons.   Although your interest is the relationship between intensity and SST you know that cyclone intensity within the hexagon will likely also depend on the number of cyclone hours.  In general a hexagon with more cyclone hours will have a higher intensity.  Thus your regression model includes SST and cyclone hours as explanatory variables in which case the SST coefficient from the regression is an estimate of the effect of SST on intensity after accounting for cyclone occurrence.

You use the \verb@lm@ function to create a linear regression model (see Chapter~\ref{chap:classicalstatistics}).  The \verb@summary@ method is used to obtain statistical information about the model.
<<global linear model, echo=TRUE, eval=TRUE>>=
model = lm(WmaxS ~ sst + count, data=hspdf)
summary(model)
@
The formula is specified as \verb@WmaxS ~ sst + count@ to indicate the mean \verb@WmaxS@ is related to \verb@sst@ and \verb@count@.  Results show SST and cyclone hours are important in statistically explaining cyclone intensity across the set of hexagons.  The parameter on the SST variable is interpreted to mean that for every 1$^\circ$C increase in sea-surface temperature cyclone intensity goes up by \Sexpr{round(model$coefficients[2],2)} m~s$^{-1}$ after accounting for cyclone hours.  The model explains \Sexpr{round(summary(model)$r.squared*100,1)}\% of the variation in cyclone intensity.  These results represent the average relationship over the basin.

\subsection{Geographically-weighted regression}

The SST parameter is replaced by a vector of parameters in GWR one for each hexagon.  The relationship between the response vector and the explanatory variables is expressed mathematically as
\begin{equation}
y = X\beta(g)+ \varepsilon
\end{equation}
where $g$ is a vector of geographic locations, here the set of hexagons with cyclone intensities and
\begin{equation}
\hat \beta(g) = (X^TWX)^{-1}X^TW y
\end{equation}
where $W$ is a weights matrix given by
\begin{equation}
W =\exp(-D^2/h^2)
\end{equation}
where $D$ is a matrix of pairwise distances between the hexagons and $h$ is the bandwidth.  The elements of the weights matrix, $w_{ij}$, are proportional to the influence hexagons $j$ have on hexagons $i$ in determining the relationship between $X$ and $y$.  Weights are determined by an inverse-distance function (kernel) so that values in nearby hexagons have greater influence on the local relationship compared with values in hexagons farther away.

The bandwidth controls the amount of smoothing.  It is chosen as a trade-off between variance and bias.   A bandwidth too narrow (steep gradients on the kernel) results in large variations in the parameter estimates (large variance).  A bandwidth too wide leads to a large bias as the parameter estimates are influenced by processes that do not represent the local conditions.

Functions for selecting a bandwidth and running GWR are available in the {\bf spgwr} package \citep{BivandYu2011}.  First require the package and select the bandwidth using the \verb@gwr.sel@ function.  The first argument is the model formula as in \verb@lm@ and the second is the data frame.  The data frame can be a spatial points or spatial polygons data frame.
<<select bandwidth, eval=TRUE, echo=TRUE, results=hide>>=
require(spgwr)
bw = gwr.sel(WmaxS ~ sst + count, data=hspdf)
bw * .001
@
The procedure is an iterative optimization with improvements made based on previous values of the bandwidth. Values for the bandwidth and the cross-validation score are printed. After several iterations, no improvement in the score occurs.  The bandwidth has dimensions of length representing a spatial distance.  The units of the LCC projection in the spatial polygons data frame is meters so to convert the bandwidth to kilometers you multiply by 10$^{-3}$.

GWR is performed using the \verb@gwr@ function.  The first two arguments are the same as in the function to select the bandwidth.  The value of the bandwidth is supplied with the bandwidth argument.
<<gwr, echo=TRUE, eval=TRUE>>=
model = gwr(WmaxS ~ sst + count, data=hspdf,
  bandwidth=bw)
@
The output is saved in an object of class {\tt gwr}.  To print a brief summary of the output type
<<gwr summary, echo=TRUE, eval=TRUE>>=
model
@

The output repeats the function call, which includes the form of the model, the kernel function (here Gaussian), and the bandwidth (units of meters).  The output also includes a summary of the regression parameter values across the hexagons.  In general, intensity is positively related to SST, but the minimum parameter value indicates that in at least one hexagon the relationship is negative.  The units on this regression parameter are m~s$^{-1}$/$^\circ$C.  The summary also includes the parameter values from a standard regression under the column heading {\tt Global}.  These values are the same as output above using the \verb@lm@ function.

The object \verb@model$SDF@ inherits the spatial polygons data frame class from the model object along with the corresponding map projection.  A summary method on this object provides additional information about the GWR.
<<summary SDF, echo=TRUE, eval=TRUE>>=
summary(model$SDF)
@
Here the coordinate bounding box is given along with the projection information.  Data attributes are stored in the data slot.  The attributes include the model coefficients (including the intercept term), the sum of the weights, predicted values, prediction errors, and R-squared values.  A histogram of the R-squared values is made by typing
<<histogram of R2, echo=TRUE, eval=FALSE>>=
hist(model$SDF$localR2, main="",
  xlab="Local R-Squared Values")
@
\begin{figure}
\centering
<<histogramR2, fig=TRUE, echo=FALSE, eval=TRUE, width=3.6, height=3.6>>=
par(las=1, mgp=c(2, .4, 0), tcl=-.3)
hist(model$SDF$localR2,
   main="", xlab="R-squared")
@
\vspace{-.5cm}
\caption{R-squared values from a GWR of intensity on SST and cyclone hours.}
\label{fig:histogramR2}
\end{figure}
The R-squared values are centered near a value 0.5, but there are \Sexpr{sum(model$SDF$localR2>.9)} hexagons with a value above 0.9 (Fig.~\ref{fig:histogramR2}).

Additional insight is obtained by mapping the results.  For instance, to make a choropleth map of the SST parameter values, type
<<SST coefficients map, fig=FALSE, echo=TRUE, eval=FALSE>>=
spplot(model$SDF, "sst", col="white",
  col.regions=blue2red(10), at=seq(-25, 25, 5),
  sp.layout=list(l2), colorkey=list(space="bottom"),
   sub="SST Effect on Intensity (m/s per C)")
@
The results are shown in Fig.~\ref{fig:sstcoef}.  The marginal influence of SST on cyclone intensity is shown along with corresponding $t$ values.  The $t$ value is calculated as the ratio of the parameter value to its standard error.  Large values of $|t|$ indicate a statistical significant relationship.  Standard errors are available by specifying the argument \verb@hatmatrix=TRUE@ in the \verb@gwr@ function.
\begin{figure}
\centering
<<SSTcoefficientsmapfigure, fig=TRUE, echo=FALSE, eval=TRUE, width=4.6, height=6>>=
model = gwr(WmaxS ~ sst + count, data=hspdf, bandwidth=bw, hatmatrix=TRUE)
p1 = spplot(model$SDF, "sst", col="white", col.regions=blue2red(10),
  at=seq(-25, 25, 5), sp.layout=list(l2),
  colorkey=list(space="bottom", labels=list(cex=.8)),
  sub=list("SST effect on intensity [m s$^{-1}$/$^\\circ$C]", cex=.8, font=1))
model$SDF$sig = model$SDF$sst / model$SDF$sst_se
al = colorRampPalette(c("blue","white","red"), space="Lab")
p2 = spplot(model$SDF, "sig", col="white", col.regions=al(10),
  at=seq(-5, 5), sp.layout=list(l2),
  colorkey=list(space="bottom", labels=list(cex=.8)),
  sub=list("Significance of effect ($t$-value)", cex=.8, font=1))
plot(update(p1, main=textGrob("a", x=unit(.05, "npc"))),
  split=c(1, 1, 1, 2), more=TRUE)
plot(update(p2, main=textGrob("b", x=unit(.05, "npc"))),
  split=c(1, 2, 1, 2), more=FALSE)
@
\vspace{0cm}
\caption{Effect of SST on cyclone intensity. (a) SST coefficient and (b) $t$-value.}
\label{fig:sstcoef}
\end{figure}

Hexagons are colored according to the SST coefficient value.  The SST coefficient represents a local `trend' of intensity as a function of SST holding cyclone occurrence constant.  Hexagons with positive coefficients, indicating a direct relationship between cyclone strength and ocean warmth in m~s$^{-1}$/$^\circ$C, are displayed using red hues and those with negative coefficients are shown with blue hues.  The divergent color ramp \verb@blue2red@ creates the colors.

Hexagons with the largest positive coefficients (greater than 5~m~s$^{-1}$/$^\circ$C) are found over the Caribbean Sea extending into the southeast Gulf of Mexico and east of the Lesser Antilles.  Coefficients above zero extend over much of the Gulf of Mexico northeastward into the southwestern Atlantic.  A region of coefficients less than zero is noted over the central North Atlantic extending from the middle Atlantic coast of the United States eastward across Bermuda and into the central Atlantic.  This map quantifies the grouped results shown in Fig.~\ref{fig:grouped}.

Local statistical significance is estimated by dividing the coefficient value by its standard error.  The ratio, called the $t$ value, is described by a $t$ distribution under the null hypothesis of a zero coefficient value (see Chapter~\ref{chap:classicalstatistics}).  Regions of high $t$ values (absolute value greater than 2) denote areas of statistical significance and generally correspond with regions of large upward trends including most of the Caribbean sea and the eastern Gulf of Mexico.  Regions with negative trends over the central Atlantic extending to the coastline fail to show significance.  To some degree these results depend on the size of your hexagons (modifiable areal unit problem).  One strategy is to rerun the GWR model with larger and smaller hexagons.  In this way you can check how much the results change and whether the changes influence your conclusions.

\subsection{Model fit}

To assess model fit, you examine the residuals.  Here a residual is the difference between observed and predicted intensity in each hexagon.  Residuals above zero indicate the model under predicts the intensity and residuals below zero indicate the model over predicts intensity.  Residuals are saved in \verb@model$SDF$gwr.e@.  A histogram of the residuals (Fig.~\ref{fig:reshist}) shows the values are centered on zero, but with a long left tail.  The long left tail indicates the model anomalously over predicts intensity in some hexagons.  Use the code below to check yourself.
<<residual histogram, echo=TRUE, eval=FALSE>>=
rsd = model$SDF$gwr.e
hist(rsd, main="", xlab="Residual (m/s)")
hist(rsd[rsd <= -10], add=TRUE, col="brown")
@
\begin{figure}
\centering
<<residualhist, fig=TRUE, echo=FALSE, eval=TRUE, width=3.4, height=3.4>>=
rsd = model$SDF$gwr.e
par(las=1, mgp=c(2, .4, 0), tcl=-.3)
hist(rsd, main="", las=1, xlab="Residual [m s$^{-1}$]")
hist(rsd[rsd <= -10], add=TRUE, col="brown")
@
\vspace{-.5cm}
\caption{Histogram of GWR residuals. Residuals less than $-$10 m~s$^{-1}$ are brown.}
\label{fig:reshist}
\end{figure}

An interesting question concerns the location of these anomalous predictions.  This is answered by coloring the hexagons having the largest negative residuals.  First create a new vector for the spatial polygons data frame consisting of 1's indicating largest negative residuals and 0's indicating otherwise.  You then use the \verb@spplot@ method and restrict the color region argument to \verb@"brown"@ and \verb@"white"@ and turn off the color key.
<<residual map, fig=FALSE, echo=TRUE, eval=FALSE>>=
model$SDF$op = as.integer(rsd <= (-10))
spplot(model$SDF, "op", col="white",
  col.regions=c("white", "brown"),
  sp.layout=list(l2), colorkey=FALSE)
@
\begin{figure}
\centering
<<residualmapfigure, fig=TRUE, echo=FALSE, eval=TRUE, width=4.1, height=2.8>>=
model$SDF$op = as.integer(rsd <= (-10))
print(spplot(model$SDF, "op", col="white", col.regions=c("white", "brown"),
  sp.layout=list(l2), colorkey=FALSE))
@
\vspace{-.5cm}
\caption{Regions of largest negative residuals.}
\label{fig:resmap}
\end{figure}
Results are shown in Fig.~\ref{fig:resmap}.  In general, the largest over predictions occur at low latitudes in hexagons near and over land.  This makes sense; in these regions, although SST values are high, cyclone intensities are limited by influences associated with land like drier air and greater surface friction.  The map suggests a way the model can be improved.  A factor indicating the presence or absence of land or a covariate indicating the proportion of hexagon covered by land would be a reasonable thing to try.

Other types of spatial regression models fit this hexagon framework.  For instance if interest is cyclone counts then a Poisson spatial model \citep{JaggerEtAl2002} can be used.  And if interest is prediction rather than explanation a model that includes spatial autocorrelation through a correlated error term or a lag response variable is possible.  In Chapter~\ref{chap:bayesianmodels} we show you how to use the hexagon framework to construct a space-time model for hurricane occurrence.

\section{Spatial Interpolation}

Spatial data often needs to be interpolated.  Cyclone rainfall is a good example. You know how much it rained at observing sites but not everywhere across the affected region.  Spatial interpolation estimates the rainfall values at arbitrary locations.  Estimates on a regular grid are used to create contour maps.

This can be done in a various ways.  Here we show you how to do it statistically.  Statistical interpolation is preferable because it includes uncertainty estimates.  The procedure is called `kriging' after a mining engineer \citep{Matheron1963}.

A kriged estimate (prediction) of some variable $z$ at a given location is a weighted average of the $z$ values over the entire domain where the weights are proportional to the spatial correlation.  The estimates are optimal in the sense that they minimize the variance between the observed and interpolated values.  In short, kriging involves estimating and modeling the spatial autocorrelation and then using the model together with the observations to interpolate.

Here you work through an example using rainfall from tropical cyclone Fay in 2008.  Fay formed near the Dominican Republic as a tropical wave, passed over the island of Hispaniola and Cuba before making landfall on the Florida Keys.  Fay then crossed the Florida peninsula and moved westward across portions of the Florida panhandle producing heavy rains in parts of the state.

\subsection{Preliminaries}

Here you will use a spatial interpolation model to generate a continuous isohyet surface of rainfall totals from Fay.  The data are in {\it FayRain.txt} and are a compilation of reports from the U.S. NOAA/NWS official weather sites and coop sites.  The coop sites are the Community Collaborative Rain, Hail and Snow Network (CoCoRaHS), a community-based, high density precipitation network made up of volunteers who take measurements of precipitation in their backyards.  The data were obtained from NOAA/NCEP/HPC and from the Florida Climate Center.

You make use of functions in the {\bf gstat} package \citep{Pebesma2004}.  Require the package and read in the data by typing
<<require gstat and read data, echo=TRUE, eval=TRUE>>=
require(gstat)
FR = read.table("FayRain.txt", header=TRUE)
names(FR)
@
The data frame contains \Sexpr{dim(FR)[1]} rainfall sites.  Longitude and latitude coordinates of the sites are given in the first two columns and total rainfall in inches and millimeters are given in the second two columns.

Create a spatial points data frame by specifying columns that contain the spatial coordinates.  Then assign a geographic coordinate system and convert the rainfall from millimeters to centimeters.
<<make a spdf, echo=TRUE, eval=TRUE>>=
coordinates(FR) = c("lon", "lat")
ll2 = "+proj=longlat +datum=NAD83"
proj4string(FR) = CRS(ll2)
FR$tpm = FR$tpm/10
@
A summary of the observed rainfall totals in cm is given by typing
<<data summary, echo=TRUE, eval=TRUE>>=
summary(FR$tpm)
@
The median value is \Sexpr{round(fivenum(FR$tpm)[3],1)}~cm and the maximum is \Sexpr{round(fivenum(FR$tpm),1)[5]}~cm.

Next import the Florida shapefile containing the county polygons as used in Chapter~\ref{chap:graphsandmaps}.  You use the \verb@readShapeSpatial@ function from the {\bf maptools} package to accomplish this task (see Chapter~\ref{chap:graphsandmaps}).
<<read FL map boundaries, echo=TRUE, eval=TRUE>>=
require(maptools)
FLpoly = readShapeSpatial("FL/FL",
  proj4string=CRS(ll2))
@

Next create a character string specifying the tags for a planar projection and transform the geographic coordinates of the site locations and map polygons to the projected coordinates.  Here you use Albers equal-area projection with true scale at latitudes 23 and 30$^\circ$N and you include a tag to specify meters as the unit of distance.
<<Alberts proj4string, echo=TRUE, eval=TRUE>>=
require(rgdal)
aea = "+proj=aea +lat_1=23 +lat_2=30 +lat_0=26,
  +lon_0=-83 +units=m"
FR = spTransform(FR, CRS(aea))
FLpoly = spTransform(FLpoly, CRS(aea))
@

A map of the rainfall sites and storm totals with the state boundaries is made by typing
<<station map, fig=FALSE, echo=TRUE, eval=TRUE>>=
l3 = list("sp.polygons", FLpoly, lwd=.3,
  first=FALSE)
spplot(FR, "tpm", sp.layout=l3)
@
Two areas of extreme rainfall are noted.  One running north-south along the east coast and another one over the north.  Rainfall collection sites are clustered in urban areas.  With the preliminaries out of the way you are ready to interpolate.

\subsection{Empirical variogram}

Rainfall is an example of geostatistical data.  In principle it can be measured anywhere, but typically you have values at a sample of sites.  The pattern of sites is not of much interest as it is a consequence of constraints (convenience, opportunity, economics, etc) unrelated to the phenomenon.  Instead, your interest centers on inference about how much rain fell across the region.  This is done using kriging---the cornerstone of geostatistics \cite{Cressie1993}.

The tradition with kriging is to model the spatial correlation using a variogram rather than a correlation function (correlogram). The sample variogram is given as
\begin{equation}
\hat \gamma(h) = \frac{1}{2N(h)}\sum^{N(h)}_{i,j} (z_i -z_j)^2
\label{eq:samplevariogram}
\end{equation}
where $N(h)$ is the number of distinct pairs of observation sites a lag distance $h$ apart, and $z_i$ and $z_j$ are the rainfall totals at sites $i$ and $j$.  Note that the $h$ is an approximate distance used with a lag tolerance $\delta h$.  The assumption is that the rainfall field is stationary. This means that the relationship between rainfall at two locations depends only on the relative positions of the sites and not on where the sites are located.  This relative position refers to lag distance and orientation between events.

You further assume the variance of rainfall between sites depends only on their lag distance and not their orientation relative to one another (isotropy assumption).  Smaller variances are expected between nearby sites compared with variances between more distant sites.  The variogram is the inverse of the spatial correlation.  You expect larger correlation in rainfall amounts between sites that are nearby and smaller correlation between sites farther apart.  This means for instance that if a site has a large rain total, nearby sites will also tend to have large totals (small difference), but farther away the totals decrease.

By definition $\gamma(h)$ is the semivariogram and $2\gamma(h)$ is the variogram.  However, for conciseness $\gamma(h)$ is often referred to as the variogram.  The variogram plots the semivariance as a function of lag distance.  Since your rainfall values have units of centimeters, the units of the semivariance are cm$^2$.

The sample variogram is computed using the \verb@variogram@ function.  The first argument is the model formula specifying the rainfall column from the data frame and the second argument is the data frame name.  Here \verb@~1@ in the model formula indicates no covariates or trends in the data.  Trends can be included by specifying coordinate names (e.g., \verb@~lon+lat@).  Recall that although your data have a planar projection, the coordinate names are not changed from the original data set.

You compute the empirical variogram for Fay's rainfall and save it by typing
<<sample variogram, echo=TRUE, eval=TRUE>>=
v = variogram(tpm ~ 1, data=FR)
@
To see the variogram values as a function of lag distance use the plot method on the variogram object.
<<plot sample variogram, echo=TRUE, eval=FALSE>>=
plot(v)
@
Values start low (50 cm$^2$) at short lag distance, then increase to over 200~cm$^2$ at lag distance of about 200~km (Fig.~\ref{fig:samplevariogram}).  Data are from \Sexpr{dim(FR@data)[1]} observations over Florida and adjacent regions. The zero-lag semivariance is called the `nugget' and the semivariance level where the variogram values no longer increase is called the `sill.'  The lag distance to the sill is called the `range.'  These three empirical parameters (nugget, sill, and range) are used to model the variogram.
\begin{figure}
\centering
<<samplevariogramfigure, echo=FALSE, eval=TRUE, fig=TRUE, width=3.8, height=3.4>>=
par(las=1, mgp=c(2, .4, 0), tcl=-.3)
plot(v$dist/1000, v$gamma, xlim=c(0, 400), ylim=c(0, 220),
   xlab="Lagged distance (h) [km]",
   ylab="Semivariance ($\\gamma$) [cm$^2$]",
   las=1, pch=19)
grid()
points(v$dist/1000, v$gamma, pch=19)
text(v$dist/1000, v$gamma, pos=1, labels=as.character(v$np), cex=.5)
@
\vspace{-.5cm}
\caption{Empirical variogram of storm-total rainfall from tropical cyclone Fay.}
\label{fig:samplevariogram}
\end{figure}

Semivariances are calculated using all pairs of rainfall observations within a lag distance (plus a lag tolerance). The number of pairs (indicated below each point on the plot) varies with lag distance.  There are many pairs at short range but fewer at longer ranges.

You check on the assumption of isotropy by plotting empirical variograms at specified azimuths.  For example, separate variograms are estimated for four directions (0, 45, 90, and 135$^\circ$) using the argument \verb@alpha@ in the \verb@variogram@ function, where 0$^\circ$ is north-south.  These directional variograms restrict observational pairs to those having similar relative orientation within an angle tolerance (think of a pie slice).  If the directional variograms look similar, then the assumption of isotropy is valid.

\subsection{Variogram model}

Next you fit a model to the empirical variogram.  The variogram model is a mathematical relationship defining the semivariance as a function of lag distance.  There are several choices for model type (functions).  To see the selection of functions as a trellis plot, type \verb@show.vgms()@.  The nugget is shown as an open circle.  The Gaussian function (`Gau') appears reasonable for your empirical variogram  since the variogram increases slowly at short lag distances then more rapidly at larger distances before leveling off at the sill.

You save the function and initial parameter values in a variogram model object by typing
<<save initial model, echo=TRUE, eval=TRUE>>=
vmi = vgm(model="Gau", psill=150, range=200 * 1000,
  nugget=50)
@
The \verb@psill@ argument is the partial sill as the difference between the sill and the nugget.  You get eyeball estimates of the parameter values from your empirical variogram.

You then use the \verb@fit.variogram@ function to fit the model to the empirical variogram.  Specifically, given the Gaussian function together with initial parameter values, the default method of weighted least squares improves the parameter estimates.  Ordinary least squares is not appropriate as the semivariances are correlated across the lag distances and the precision on the estimates varies depending on the number of site pairs for a given lag.
<<fit variogram model, echo=TRUE, eval=TRUE>>=
v.fit = fit.variogram(v, vmi)
v.fit
@
The result is a variogram model with a nugget of \Sexpr{round(v.fit$psill[1],1)}~cm$^2$, a partial sill of \Sexpr{round(v.fit$psill[2],1)}~cm$^2$, and a range of \Sexpr{round(v.fit$range[2]/1000,1)}~km. You plot the variogram model on top of the empirical variogram by typing
<<variogram model, echo=TRUE, eval=FALSE>>=
plot(v, v.fit)
@
Let $r$ be the range, $c$ the partial sill and $c_o$ the nugget, then the equation defining the curve over the set of lag distances $h$ in Fig.~\ref{fig:gaussianvariogrammodel} is
\begin{equation}
\gamma(h)=c\left(1-\exp\left(-\frac{h^2}{r^2}\right)\right)+c_o
\label{eq:gaussianvariogrammodel}
\end{equation}
\begin{figure}
\centering
<<plotmodelvariogramfigure, fig=TRUE, echo=FALSE, eval=TRUE, width=3.8, height=3.4>>=
par(las=1, mgp=c(2, .4, 0), tcl=-.3)
plot(v$dist/1000, v$gamma, xlim=c(0, 400), ylim=c(0, 220),
   xlab="Lag distance (h) [km]",
   ylab="Semivariance ($\\gamma$) [cm$^2$]",
   las=1, pch=19)
grid()
nug = v.fit$psill[1]
ps = v.fit$psill[2]
r = v.fit$range[2] / 1000
h = seq(0, 400, .2)
fit = ps * (1 - exp(-h^2 / (r^2))) + nug
lines(h, fit, lwd=2)
@
\vspace{-.5cm}
\caption{Gaussian variogram model for Fay's rainfall across Florida.}
\label{fig:gaussianvariogrammodel}
\end{figure}

There are a variety of variogram functions.  You can try the spherical function by replacing the \verb@model="Gau"@ with \verb@model="Sph"@ in the \verb@vgm@ function above.
<<spherical model, echo=TRUE, eval=FALSE>>=
vmi2 = vgm(model="Sph", psill=150,
   range=200 * 1000, nugget=50)
v2.fit = fit.variogram(v, vmi2)
plot(v, v2.fit)
@
The {\bf geoR} package contains the \verb@eyefit@ function that can make choosing a function easier however, the interpolated values are not overly sensitive assuming a reasonable choice.

\subsection{Kriging}

The final step is to use the variogram model together with the rainfall values at the observation sites to create an interpolated surface.  The process is called kriging.  As Edzer Pebesma notes, `krige' is to `kriging' as `predict' is to `prediction.'    Here you use ordinary kriging as there are no spatial trends in the rainfall.  Universal kriging is used when trends are present.

Interpolation is done using the \verb@krige@ function.  The first argument is the model specification and the second is the data.  Two other arguments are needed.  One is the variogram model using the argument name \verb@model@ and the other is a set of locations identifying where the interpolations are to be made.  This is specified with the argument name \verb@newdata@.

Here you interpolate first to locations (point kriging) on a regular grid and then to the county polygons (block kriging).  To create a grid of locations within the boundary of Florida type
<<prediction grid, echo=TRUE, eval=TRUE>>=
grd = spsample(FLpoly, n=5000, type="regular")
@
You specify the number of locations using the argument \verb@n@.  The actual number will be slightly different because of the irregular boundary.  To view the locations simply use the plot method on the grid object.  More options for plotting are available if the spatial points object is converted to a spatial pixels object.
<<points2pixels, echo=TRUE, eval=TRUE>>=
grd = as(grd, "SpatialPixels")
@

First use the \verb@krige@ function to interpolate (predict) the observed rainfall to the grid locations.  For a given location, the interpolation is a weighted average of the rainfall across the entire region where the weights are determined by the variogram model.
<<point krige interpolation>>=
ipl = krige(tpm ~ 1, FR, newdata=grd, model=v.fit)
@
The function recognizes the type of kriging being performed.  For instance, if the variogram model is not included then inverse distance weighted interpolation is performed.  Note that the function will not work if different values have the same location.

The object (\verb@ipl@) inherits the spatial pixels object specified in the \verb@newdata@ argument, but extends it to a spatial pixels data frame by adding a data slot.  The data slot is a data frame with two variables.  The first {\tt var1.pred} is the interpolated rainfall and the second {\tt var1.var} is the prediction variance.

You plot the interpolated field using the \verb@spplot@ method.  You specify an inverted topographic color ramp to highlight in blue regions with the highest rain totals.
<<plot interpolations, echo=TRUE, eval=FALSE>>=
spplot(ipl, "var1.pred", col.regions=
  rev(topo.colors(20)), sp.layout=l3)
@
The map (see Fig.~\ref{fig:krigerainfallmaps}) makes it easy to see that parts of east central and north Florida were deluged by Fay.

You use block kriging to estimate rainfall amounts within each county.  The county-wide rainfall average is relevant for water resource managers.  Block kriging produces a smoothed estimate of this area average, which will differ from a simple arithematic average over all sites within the county.

You use the same function to interpolate but specify the spatial polygons rather than the spatial grid as the new data.
<<block krige interpolation, echo=TRUE, eval=TRUE>>=
ipl2 = krige(tpm ~ 1, FR, newdata=FLpoly,
  model=v.fit)
@
You then map the results using the same plot arguments as before.
<<plot county averages, echo=TRUE, eval=FALSE>>=
spplot(ipl2, "var1.pred", col.regions=
  rev(topo.colors(20)))
@
Rainfall maps using point and block kriging are shown in Fig.~\ref{fig:krigerainfallmaps}.
\begin{figure}
\centering
<<krigerainfallmaps, fig=TRUE, echo=FALSE, eval=TRUE, width=4.3, height=3.1>>=
mar = 30000
p1 = spplot(ipl, "var1.pred", col.regions=rev(topo.colors(20)),
  sp.layout=l3, at=seq(0, 40, 5),
  colorkey=list(space="bottom", labels=list(cex=.7)),
  sub=list("Storm total rainfall [cm]", cex=.7, font=1),
  xlim=c(FLpoly@bbox[1] - mar, FLpoly@bbox[3] + mar),
  ylim=c(FLpoly@bbox[2] - mar, FLpoly@bbox[4] + mar))
p2 = spplot(ipl2, "var1.pred", col.regions=rev(topo.colors(20)),
  at=seq(0, 40, 5),
  colorkey=list(space="bottom", labels=list(cex=.7)),
  sub=list("County average rainfall [cm]", cex=.7, font=1),
  xlim=c(FLpoly@bbox[1] - mar,
  FLpoly@bbox[3] + mar), ylim=c(FLpoly@bbox[2] - mar, FLpoly@bbox[4] + mar))
p1 = update(p1, main=textGrob("a", x=unit(.05, "npc")))
p2 = update(p2, main=textGrob("b", x=unit(.05, "npc")))
plot(p1, split=c(1, 1, 2, 1), more=TRUE)
plot(p2, split=c(2, 1, 2, 1), more=FALSE)
@
\vspace{0cm}
\caption{Rainfall from tropical cyclone Fay.  (a) Point and (b) block kriging.}
\label{fig:krigerainfallmaps}
\end{figure}
The overall pattern of rainfall from Fay featuring the largest amounts along the central east coast and over the eastern panhandle are similar in both maps.

You compute the arithematic average of county-wide rainfall again using the \verb@over@ function by typing
<<arithmetic average rainfall, echo=TRUE, eval=TRUE>>=
ipl3 = over(x=FLpoly, y=FR, fn=mean)
@
The function returns a data frame of the average rainfall in each county.  The state-wide mean of the kriged estimates is \Sexpr{round(mean(ipl2$var1.pred),1)}~cm, which compares with a state-wide mean of the arithematic averages of \Sexpr{round(mean(ipl3$tpm),1)}~cm.  The correlation between the two estimates across the 67 counties is \Sexpr{round(cor(ipl3$tpm, ipl2$var1.pred),2)}.  The variogram model reduces the standard deviation of the kriged estimate (\Sexpr{round(sd(ipl2$var1.pred),2)}~cm) relative to the standard deviation of the simple averages (\Sexpr{round(sd(ipl3$tpm),2)}~cm) because of the smoothing.

\subsection{Uncertainty}

An advantage of kriging as a method of spatial interpolation is the accompanying uncertainty estimates.  The prediction variances are listed in a column in the spatial data frame saved from your application of the \verb@krige@ function. Variances are smaller in regions with a greater number of rainfall observations.  Prediction variances are also smaller with block kriging as much of the variability within the county averages out.  To compare the distribution characteristics of the prediction variances for the point and block kriging of the rainfall observations, type
<<prediction variances, echo=TRUE, eval=TRUE>>=
round(summary(ipl$var1.var), 1)
round(summary(ipl2$var1.var), 1)
@
The median prediction variance for your point kriging is \Sexpr{round(fivenum(ipl$var1.var)[3],1)}~cm$^2$, which is close to the value of the nugget.  In contrast the median prediction variance for your block kriging is a much smaller \Sexpr{round(fivenum(ipl2$var1.var)[3],1)}~cm$^2$.

Simulations exploit this uncertainty and provide synthetic data for use in deterministic models.  A rainfall-runoff model, for instance, can be run with simulated rainfall fields providing a realistic representation of the variation in the amount of runoff resulting from the uncertainty in the rainfall field \citep{BivandEtAl2008}.

Conditional simulation, where the simulated field (realization) is generated given the data and the variogram model, is done using the same \verb@krige@ function by adding the argument \verb@nsim@ to specify the number of simulations.  For a large number it may be necessary to limit the number neighbors in the kriging.  This is done using the \verb@nmax@ argument.  For a given location, the weights assigned to observations far away are very small, so it is efficient to limit how many are used in the simulation.

As an example, here you generate four realizations of the county level storm total rainfall for Fay and limit your neighborhood to 50 of the closest observation sites.  Note that it may take a few minutes to finish processing this function.
<<simulations, echo=TRUE, eval=TRUE, results=hide>>=
ipl.sim = krige(tpm ~ 1, FR, newdata=FLpoly,
  model=v.fit, nsim=4, nmax=50)
@
Maps of the four realizations are shown in Fig.~\ref{fig:simulationsplot}.  Simulations are conditional on the observed rainfall and the variogram model using block kriging on the counties.  Note the overall pattern of rainfall remains the same, but there are differences especially in counties with relative few observations and where the rainfall gradients are steep.   The functionality is limited to Gaussian simulation, but the {\bf RandomFields} package \citep{Schlather2011} provides additional simulation algorithms.
\begin{figure}
\centering
<<simulationsplot, echo=FALSE, eval=TRUE, fig=TRUE, width=4.7, height=5.2>>=
mar = 30000
p1 = spplot(ipl.sim, "sim1", col.regions=rev(topo.colors(20)), at=seq(0, 40, 5),
  colorkey=list(space="bottom", labels=list(cex=.7)),
  sub=list("Simulated rainfall [cm]", cex=.7, font=1),
  xlim=c(FLpoly@bbox[1] - mar,
  FLpoly@bbox[3] + mar), ylim=c(FLpoly@bbox[2] - mar, FLpoly@bbox[4] + mar))
p2 = spplot(ipl.sim, "sim2", col.regions=rev(topo.colors(20)), at=seq(0, 40, 5),
  colorkey=list(space="bottom", labels=list(cex=.7)),
  sub=list("Simulated rainfall [cm]", cex=.7, font=1),
  xlim=c(FLpoly@bbox[1] - mar,
  FLpoly@bbox[3] + mar),ylim=c(FLpoly@bbox[2] - mar, FLpoly@bbox[4] + mar))
p3 = spplot(ipl.sim, "sim3", col.regions=rev(topo.colors(20)), at=seq(0, 40, 5),
  colorkey=list(space="bottom", labels=list(cex=.7)),
  sub=list("Simulated rainfall [cm]", cex=.7, font=1),
  xlim=c(FLpoly@bbox[1] - mar,
  FLpoly@bbox[3] + mar), ylim=c(FLpoly@bbox[2] - mar, FLpoly@bbox[4] + mar))
p4 = spplot(ipl.sim, "sim4", col.regions=rev(topo.colors(20)), at=seq(0, 40, 5),
  colorkey=list(space="bottom", labels=list(cex=.7)),
  sub=list("Simulated rainfall [cm]", cex=.7, font=1),
  xlim=c(FLpoly@bbox[1] - mar,
  FLpoly@bbox[3] + mar), ylim=c(FLpoly@bbox[2] - mar, FLpoly@bbox[4] + mar))
plot(p1, position=c(0, .5, .5, 1))
plot(p2, position=c(.5, .5, 1, 1), newpage=FALSE)
plot(p3, position=c(0, 0, .5, .5), newpage=FALSE)
plot(p4, position=c(.5, 0, 1, .5), newpage=FALSE)
@
\vspace{0cm}
\caption{Rainfall simulations from tropical cyclone Fay.}
\label{fig:simulationsplot}
\end{figure}

In summary, kriging is a statistical method for spatial data interpolation involving three steps.  First you estimate an empirical variogram that describes the spatial autocorrelation structure of your observations.  This step includes checking for trends and isotropy.  Second you determine a variogram model using the method of weighted least squares.  This step involves a somewhat subjective choice of the model but the interpolation is not too sensitive to this choice assuming it is reasonable.  Third you use the variogram model together with your observed values to interpolate values at specified locations, on a grid, or over an area.  Finally, simulation methods generate synthetic data for expressing and propagating uncertainty associated with the spatial variability in the data and with the process of specifying a model.  Kriging is a nuanced process, but with practice it can be an important tool in your toolbox.

This chapter demonstrated ways to analyze and model hurricane data spatially.  We began by showing how to combine hurricane track data with climate data. We then showed how to analyze the data including estimating the spatial autocorrelation.  This was followed by a demonstration on how to use geographically-weighted regression (GWR) to model hurricane intensity across space.  Quantifying and mapping the changing relationship between the response variable and covariates can lead to new insights about how hurricanes work.  This area of research is wide open.  Finally we showed how to perform a spatial interpolation using a prediction method from geostatistics.  Next we look at some ways to analyze and model time series data.