1

我的数据结构如下:

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

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

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)

该代码似乎有效。现在,我想使用over()sp 包中的函数从插值网格的指定单元格中提取预测。我尝试使用以下代码来完成此操作:

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_outputis aSpatialPointsDataFrame而不是 a SpatialPixelsDataFrameor SpatialGridDataFrame

有人对如何解决这个问题有任何建议吗?它可以像转换kr.grid$krige_output为 aSpatialPixelsDataFrame或 a一样简单SpatialGridDataFrame。不幸的是,我不知道该怎么做。

4

1 回答 1

2

通过在几何坐标中创建序列xy坐标并随后进行投影,您的网格不再均匀分布。先尝试投影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
于 2014-03-30T03:22:35.707 回答