Something like this
Load necessary packages
kpacks <- c('sp','rgdal', 'gstat', 'raster')
new.packs <- kpacks[!(kpacks %in% installed.packages()[,"Package"])]
if(length(new.packs)) install.packages(new.packs)
lapply(kpacks, require, character.only=T)
remove(kpacks, new.packs)
data(wrld_simpl)
a projected coordinate system to work with
p.utm33n <- CRS("+init=epsg:32633") # UTM 33N Landsat Images
A country (I particularly like this one)
ago <- wrld_simpl[wrld_simpl@data$NAME == 'Angola',]
Project it to UTM 33S
ago <- spTransform(ago, p.utm33n)
Sample some points within polygons
ago_p <- spsample(ago, type="random", n=25)
plot(ago, col = 'grey' , axes = T)
plot(ago_p, add = T)
Some imaginary temperature data for 3 days
tdata <- data.frame(x=rep(coordinates(ago_p)[,1], 3),
y=rep(coordinates(ago_p)[,2], 3),
temp=runif(75, 12,35),
day = rep(1:3, each = 25))
Manage to get it as spatialPointDataFrame object
coordinates(tdata) <- ~x+y
proj4string(tdata) <- CRS(proj4string(ago))
Since I don't know your base map, I'll use the country I pick up above The base layer must be a SpatialPixelDataBase. Ill play with a rasterLayer
rago <- raster(extent(ago))
res(rago) <- c(10000,10000)
rago[] <- 1
proj4string(rago) <- CRS(proj4string(ago))
r_ago <- mask(rago, ago)
#plot(r_ago)
grid_ago <- as(r_ago, 'SpatialPointsDataFrame')
grid_ago <- grid_ago[!is.na(grid_ago@data$layer), ]
gridded(grid_ago) <- TRUE
I can now run idw()
from gstat. I'll run with data from day == 1
idw_ago <- idw(temp ~ 1, tdata[tdata$day == 1, ], grid_ago, idp = 2.5)
And finally plot it
spplot(idw_ago, "var1.pred")
Now with your data, that I was missing from your question. The same approach
library(latticeExtra)
p.dutch <- CRS("+init=epsg:28991") # Dutch National Grid EPSG:28991
load(url('http://gadm.org/data/rda/DEU_adm0.RData'))
ger <- gadm
ger <- spChFIDs(ger, paste("ger", rownames(ger), sep = "_"))
ger <- spTransform(ger, p.dutch)
ger_p <- spsample(ger, type="random", n=25)
plot(ger, col = 'yellow', border = NA, axes = T, cex.axis = 0.6)
plot(ger_p, add = T, pch = 20)
tdata <- data.frame(x=rep(coordinates(ger_p)[,1], 3),
y=rep(coordinates(ger_p)[,2], 3),
temp=runif(75, 12,35),
day = rep(1:3, each = 25))
coordinates(tdata) <- ~x+y
proj4string(tdata) <- CRS(proj4string(ger))
rger <- raster(extent(ger))
res(rger) <- c(10000,10000)
rger[] <- 1
proj4string(rger) <- CRS(proj4string(ger))
r_ger <- mask(rger, ger)
plot(r_ger)
grid_ger <- as(r_ger, 'SpatialPointsDataFrame')
grid_ger <- grid_ger[!is.na(grid_ger@data$layer), ]
gridded(grid_ger) <- TRUE
idw_ger <- idw(temp ~ 1, tdata[tdata$day == 1, ], grid_ger, idp = 2.5)
spplot(idw_ger, "var1.pred") +
latticeExtra::layer(sp.polygons(ger, fill = NA, col = 'blue')) +
latticeExtra::layer(sp.points(tdata[tdata$day == 1, ],
fill = NA, col = 'red'))
Hope it helps