Spatial Overlay for Kriging Predictions in R

Question

I have data that are structured as follows:

winter.dat<-structure(list(ELON = c(-98.02325, -96.66909, -99.33808, -98.70974, 
-98.29216, -97.08568, -99.90308, -100.53012, -99.05847, -95.86621, 
-97.25452, -102.49713, -96.63121, -97.69394, -96.35404, -94.6244, 
-99.64101, -96.81046, -97.26918, -99.27059, -97.0033, -99.34652, 
-97.28585, -96.33309, -96.80407, -98.36274, -99.7279, -97.91446, 
-95.32596, -95.2487, -95.64617, -94.84896, -95.88553, -96.32027, 
-98.03654, -99.80344, -95.65707, -98.49766, -96.71779, -96.42777, 
-99.14234, -98.46607, -101.6013, -98.743583, -97.47978, -95.64047, 
-96.0024, -98.48151, -99.05283, -96.35595, -99.83331, -101.22547, 
-95.54011, -94.8803, -95.45067, -94.78287, -102.8782, -97.76484, 
-97.95442, -98.11139, -95.99716, -96.94394, -99.42398, -97.21271, 
-99.01109, -95.78096, -97.74577, -98.56936, -94.84437, -97.95553, 
-97.60685, -94.82275, -96.91035, -97.20142, -97.95202, -95.60795, 
-97.46488, -96.49749, -97.46414, -97.5107, -97.58245, -96.26265, 
-95.91473, -97.22924, -96.76986, -97.04831, -95.55976, -95.27138, 
-98.96038, -97.15306, -99.36001, -97.58812, -94.79805, -99.0403, 
-96.94822, -96.03706, -100.26192, -97.34146, -95.18116, -97.09527, 
-96.06982, -96.95048, -94.98671, -95.01152, -99.13755, -96.67895, 
-95.22094, -97.52109, -98.52615, -97.98815, -98.77509, -95.1233, 
-94.64496, -95.34805, -94.68778, -99.41682, -96.34222), NLAT = c(34.80833, 
34.79851, 34.58722, 36.70823, 34.91418, 34.19258, 36.07204, 36.80253, 
35.40185, 35.96305, 36.75443, 36.69256, 35.17156, 36.41201, 35.7805, 
34.0433, 36.83129, 36.63459, 33.89376, 35.5915, 34.8497, 36.02866, 
36.1473, 34.60896, 35.65282, 36.74813, 35.54615, 35.03236, 34.65657, 
34.22321, 36.32112, 35.68001, 36.90987, 33.92075, 35.54848, 35.20494, 
35.30324, 36.26353, 34.55205, 36.84053, 36.72562, 35.14887, 36.60183, 
34.239444, 35.84891, 35.74798, 35.84162, 35.48439, 34.98971, 
35.07073, 34.6855, 36.85518, 34.03084, 33.83013, 36.14246, 36.4821, 
36.82937, 34.52887, 35.85431, 36.38435, 34.30876, 34.03579, 34.83592, 
36.06434, 36.98707, 34.88231, 36.79242, 34.72921, 36.88832, 35.27225, 
36.11685, 34.31072, 36.8981, 34.2281, 34.96774, 36.74374, 35.23611, 
36.03126, 35.47226, 35.5556, 35.47112, 35.43172, 35.58211, 34.7155, 
36.36114, 35.99865, 35.8257, 36.36914, 35.89904, 36.3559, 35.12275, 
34.19365, 35.43815, 36.19033, 35.36492, 36.4153, 36.59749, 35.54208, 
35.26527, 36.12093, 34.87642, 34.5661, 35.97235, 34.7107, 34.43972, 
34.33262, 36.77536, 34.98224, 35.84185, 34.16775, 35.5083, 35.489, 
36.011, 34.90092, 34.98426, 36.42329, 36.51806), RAIN_WINTER14 = c(0.7, 
1.8, 1.63, 1.14, 1.43, 4.2, 0.76, 1.42, 0.65, 2.42, 0.95, 0.24, 
2.82, 1.33, 1.37, 7.5, 1.22, 1.65, 4.3, 0.54, 2.99, 0.78, 1.17, 
5.57, 1.85, 0.99, 0.42, 0.69, 5, 4.23, 1.17, 5.82, 1.28, 4.42, 
1.22, 0.58, 4.28, 0.85, 2.12, 1.72, 1.41, 0.93, 0.47, 2.28, 1.43, 
3.86, 2.69, 1.17, 1.17, 1.6, 1.12, 0.85, 5.27, 7.11, 1.96, 3.12, 
0.25, 1.52, 1.19, 1.07, 3.53, 4.95, 0.87, 1.32, 1.26, 4.53, 0.97, 
0.47, 2.35, 1.56, 1.22, 7.55, 1.23, 2.98, 0.53, 1.41, 0.61, 1.74, 
1.46, 1.62, 1.71, 2.18, 2.43, 1.72, 1.05, 1.39, 2.52, 1.26, 0.61, 
1.4, 1.01, 2.13, 4.95, 0.9, 1.34, 1.69, 1.29, 1.56, 4.4, 1.13, 
4.82, 2.44, 5.29, 5.68, 1.69, 5.38, 2, 2.54, 1.17, 2.21, 0.67, 
4.38, 5.86, 3.79, 6.14, 1.05, 1.03)), .Names = c("ELON", "NLAT", 
"RAIN_WINTER14"), row.names = c(NA, -117L), class = "data.frame")

sensor_points<-structure(list(LON = c(-95.91, -97.51, -98.42, -97.51, -97.34, 
-97.86, -95.95, -96.09, -96.43, -96.26, -97.11, -98.61, -95.61, 
-95.91, -95.91, -94.45, -94.6, -97.51, -95.78, -96.46, -99.5, 
-95.78, -98.42, -95.91, -95.94, -94.97, -95.38, -97.86, -97.37, 
-97.51, -94.6, -97.51, -97.47, -97.34, -95.78, -97.06, -95.91, 
-97.59, -95.66, -97.76, -96.51, -94.66, -99.5, -95.49, -97.22, 
-98.24, -95.52, -95.38, -95.26, -96.14), LAT = c(36.12, 35.46, 
34.6, 35.46, 35.22, 36.4, 35.63, 35.99, 33.93, 34.13, 34.19, 
34.62, 36.31, 36.12, 36.12, 35.33, 35.4, 35.46, 36.03, 36.29, 
34.87, 36.03, 34.6, 36.12, 36.73, 35.91, 35.95, 36.4, 35.01, 
35.46, 34.89, 35.46, 35.65, 35.22, 36.03, 36.72, 36.12, 35.24, 
35.96, 35.51, 34, 35.13, 34.87, 36.28, 34.72, 35.06, 35.47, 35.95, 
36.43, 34.01)), .Names = c("LON", "LAT"), row.names = c(NA, 50L
), class = "data.frame")

I am using the autoKrige() function from the automap package in R to generate interpolated values for a grid that I have defined using the following code:

ok_state<-map("state","ok",plot=F) 
sp1<-seq(min(ok_state$x,na.rm=T),max(ok_state$x,na.rm=T),length=100) 
sp2<-seq(min(ok_state$y,na.rm=T),max(ok_state$y,na.rm=T),length=100) 
sp<-expand.grid(sp1,sp2)
S<-SpatialPoints(sp)
gridded(S)<-TRUE
proj4string(S)<-"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
S<-spTransform(S,CRS("+proj=merc +zone=18s +ellps=WGS84 +datum=WGS84"))

To interpolate, I am using the following code:

coordinates(winter.dat)=~ELON+NLAT
proj4string(winter.dat)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
project_winter.dat=spTransform(winter.dat, CRS("+proj=merc +zone=18s +ellps=WGS84 +datum=WGS84"))
kr.grid<-autoKrige(RAIN_WINTER14~1,project_winter.dat,S)

That code seems to work. Now, I would like to extract predictions from specified cells of the interpolated grid using the over() function in the sp package. I tried to accomplish this using this code:

coordinates(sensor_points)=~LON+LAT
proj4string(sensor_points)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
project_sensor_points=spTransform(sensor_points, CRS("+proj=merc +zone=18s +ellps=WGS84 +datum=WGS84"))
proj4string(kr.grid$krige_output)==proj4string(project_sensor_points)
over(project_sensor_points,kr.grid$krige_output)

But this code yielded a matrix full of NA values. This may have happened because kr.grid$krige_output is a SpatialPointsDataFrame rather than a SpatialPixelsDataFrame or SpatialGridDataFrame.

Does anyone have any suggestions about how to address this problem? It could be as easy as converting kr.grid$krige_output to a SpatialPixelsDataFrame or a SpatialGridDataFrame. Unfortunately, I can't figure out how to do that.

Answer 1

By creating the sequences of x and y coordinates in geometric coordinates and subsequently projecting, your grid is no longer evenly spaced. Try projecting the ok_state extent first, and then create the sequence of evenly-spaced points in the new CRS.

library(maps)
library(sp)
library(rgdal)
library(automap)

merc18s <- CRS("+proj=merc +zone=18s +ellps=WGS84 +datum=WGS84")
wgs84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

ok_state <- map("state", "ok", plot=F) 
e <- data.frame(x=range(ok_state$x, na.rm=TRUE),
                y=range(ok_state$y, na.rm=TRUE))
coordinates(e) <- ~x+y
proj4string(e)<- wgs84
e.merc <- spTransform(e, merc18s)

S <- SpatialPoints(expand.grid(seq(e.merc$x[1], e.merc$x[2], length=100),
                                seq(e.merc$y[1], e.merc$y[2], length=100)))

gridded(S) <- TRUE

coordinates(winter.dat) <- ~ELON+NLAT
proj4string(winter.dat) <- wgs84
winter.dat.merc <- spTransform(winter.dat, merc18s)
kr.grid <- autoKrige(RAIN_WINTER14~1, winter.dat.merc, S)


coordinates(sensor_points) <- ~LON+LAT
proj4string(sensor_points) <- wgs84
sensor_points.merc <- spTransform(sensor_points, merc18s)

over(sensor_points.merc, kr.grid$krige_output)

This yields:

   var1.pred  var1.var var1.stdev
1  1.8893073 0.3804578  0.6168126
2  1.4171934 0.3588603  0.5990495
3  1.2733813 0.4004269  0.6327930
4  1.4171934 0.3588603  0.5990495
5  1.4940596 0.3722470  0.6101204
6  1.1237834 0.3879326  0.6228424
7  2.7571025 0.3726498  0.6104505
8  1.9355300 0.3784474  0.6151808
9  4.7656947 0.4286066  0.6546806
10 4.5619002 0.4072064  0.6381272
11 3.6205513 0.3698636  0.6081641
12 1.2280841 0.3974779  0.6304585
13 1.7566100 0.3780431  0.6148521
14 1.8893073 0.3804578  0.6168126
15 1.8893073 0.3804578  0.6168126
16 6.2926628 0.5055258  0.7110034
17 5.8679727 0.4378686  0.6617164
18 1.4171934 0.3588603  0.5990495
19 2.2931583 0.3732854  0.6109708
20 1.3982178 0.3871845  0.6222415
21 1.0272094 0.3830222  0.6188879
22 2.2931583 0.3732854  0.6109708
23 1.2733813 0.4004269  0.6327930
24 1.8893073 0.3804578  0.6168126
25 1.3115832 0.3948189  0.6283462
26 4.4892552 0.3843763  0.6199809
27 3.1293649 0.3792566  0.6158381
28 1.1237834 0.3879326  0.6228424
29 1.6571943 0.3777782  0.6146366
30 1.4171934 0.3588603  0.5990495
31 6.3472936 0.4459952  0.6678287
32 1.4171934 0.3588603  0.5990495
33 1.4031739 0.3635883  0.6029828
34 1.4940596 0.3722470  0.6101204
35 2.2931583 0.3732854  0.6109708
36 1.2112401 0.3877834  0.6227226
37 1.8893073 0.3804578  0.6168126
38 1.3631983 0.3672577  0.6060179
39 2.5845513 0.3715311  0.6095335
40 1.3130873 0.3674918  0.6062111
41 4.6494661 0.4154409  0.6445470
42 5.8924573 0.4228850  0.6502961
43 1.0272094 0.3830222  0.6188879
44 2.0579067 0.3776652  0.6145447
45 2.1571974 0.3750187  0.6123877
46 0.9718268 0.3733108  0.6109916
47 3.8462812 0.3836691  0.6194103
48 3.1293649 0.3792566  0.6158381
49 2.3244060 0.3853137  0.6207364
50 4.7678701 0.4246547  0.6516554