Пространственное наложение для прогнозов Кригинга в R
У меня есть данные, которые структурированы следующим образом:
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() функция из пакета automap в R для генерации интерполированных значений для сетки, которую я определил, используя следующий код:
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 это SpatialPointsDataFrame а не SpatialPixelsDataFrame или же SpatialGridDataFrame,
У кого-нибудь есть предложения о том, как решить эту проблему? Это может быть так же просто, как преобразование kr.grid$krige_output к SpatialPixelsDataFrame или SpatialGridDataFrame, К сожалению, я не могу понять, как это сделать.
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