Spatial interpolation of geographical data in R

Current task is to create raster map of Estonian air temperature for 12/May/2010.

This example is based on: http://rpubs.com/adam_dennett/10873

Data: www.geo.ut.ee/aasa/LOOM02331/interpolate.rar

Import estonia_air_temperature_2.csv to your RStudio R project (Tools - Import Dataset - From Text File).

We use in current example the IDW (Inverse Distnace Weighting) interpolation method. More about: http://en.wikipedia.org/wiki/Inverse_distance_weighting

Run libraries:

library(ggplot2)
library(gstat)
library(sp)
library(maptools)
## Checking rgeos availability: TRUE

Import data:

estonia_air_temperature_2 <- read.csv(file = "estonia_air_temperature_2.csv", 
    header = TRUE)

Let's duplicate the data table to avoid accidental overwriting. Define latitude and longitude:

estonia_air_temperature_2_test <- estonia_air_temperature_2  # duplicate air temp. data file
estonia_air_temperature_2_test$x <- estonia_air_temperature_2_test$lon  # define x & y as longitude and latitude
estonia_air_temperature_2_test$y <- estonia_air_temperature_2_test$lat

Set spatial coordinates to create a Spatial object:

coordinates(estonia_air_temperature_2_test) = ~x + y

Plot the results:

plot(estonia_air_temperature_2_test)

plot of chunk unnamed-chunk-5

Define the grid extent:

x.range <- as.numeric(c(21.76, 28.21))  # min/max longitude of the interpolation area
y.range <- as.numeric(c(57.45, 59.72))  # min/max latitude of the interpolation area

Create a data frame from all combinations of the supplied vectors or factors. See the description of the return value for precise details of the way this is done. Set spatial coordinates to create a Spatial object. Assign gridded structure:

grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 0.1), y = seq(from = y.range[1], 
    to = y.range[2], by = 0.1))  # expand points to grid
coordinates(grd) <- ~x + y
gridded(grd) <- TRUE

Plot the weather station locations and interpolation grid:

plot(grd, cex = 1.5, col = "grey")
points(estonia_air_temperature_2_test, pch = 1, col = "red", cex = 1)

plot of chunk unnamed-chunk-8

Interpolate surface and fix the output:

idw <- idw(formula = may12 ~ 1, locations = estonia_air_temperature_2_test, 
    newdata = grd)  # apply idw model for the data
## [inverse distance weighted interpolation]
idw.output = as.data.frame(idw)  # output is defined as a data table
names(idw.output)[1:3] <- c("long", "lat", "var1.pred")  # give names to the modelled variables

Plot the results:

ggplot() + geom_tile(data = idw.output, aes(x = long, y = lat, fill = var1.pred)) + 
    geom_point(data = estonia_air_temperature_2, aes(x = lon, y = lat), shape = 21, 
        colour = "red")

plot of chunk unnamed-chunk-10

Not to beautiful? Let's add Estonian contour, change the colour palette and add locations of weather stations:

est_contour <- readShapePoly("C:/ANTO/loengud/GIS Maps and Spatial Analyses for Urban Planning/data/population_in_municipalities_2011_wgs84.shp")
est_contour <- fortify(est_contour, region = "name")
## Loading required package: rgeos
## rgeos version: 0.3-4, (SVN revision 438)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Polygon checking: TRUE
ggplot() + geom_tile(data = idw.output, alpha = 0.8, aes(x = long, y = lat, 
    fill = round(var1.pred, 0))) + scale_fill_gradient(low = "cyan", high = "orange") + 
    geom_path(data = est_contour, aes(long, lat, group = group), colour = "grey") + 
    geom_point(data = estonia_air_temperature_2, aes(x = lon, y = lat), shape = 21, 
        colour = "red") + labs(fill = "Air temp.", title = "Air temperature in Estonia, 15/May/2010")

plot of chunk unnamed-chunk-11

And maybe you noticed, that the interpolated surface is not continuous any more. This is beacause I rounded all interpolated values: fill=round(var1.pred, 0). More about using ggplot2: http://docs.ggplot2.org/current/

Author: Anto Aasa, anto.aasa@ut.ee

Guide is compiled in RStudio with R Markdown.