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

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)

そのコードはうまくいくようです。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_outputために発生した可能性があります。SpatialPointsDataFrameSpatialPixelsDataFrameSpatialGridDataFrame

この問題に対処する方法について何か提案はありますか? または に変換kr.grid$krige_outputするのと同じくらい簡単です。残念ながら、その方法がわかりません。SpatialPixelsDataFrameSpatialGridDataFrame

4

1 に答える 1

2

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
于 2014-03-30T03:22:35.707 に答える