R中Kriging预测的空间叠加

Question

我的数据结构如下：

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")

我使用R中的automap包中的autoKrige()函数为我使用以下代码定义的网格生成插值值：

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"))

要插值，我使用以下代码：

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)

这段代码似乎有效。现在，我想使用sp包中的over()函数从插值网格的指定单元中提取预测。我试图使用以下代码来完成这一任务：

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)

但是这段代码产生了一个充满NA值的矩阵。这可能是因为kr.grid$krige_output是SpatialPointsDataFrame而不是SpatialPixelsDataFrame或SpatialGridDataFrame。

有没有人对如何解决这个问题有任何建议？它可以像将kr.grid$krige_output转换为SpatialPixelsDataFrame或SpatialGridDataFrame一样容易。不幸的是，我不知道该怎么做。

Answer 1

通过在几何坐标中创建x和y坐标的序列并随后投影，您的网格不再是均匀的。首先尝试投影ok_state范围，然后在新的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)

这产生了：

   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