Spatial Overlay for Kriging Predictions in R -


i have data structured 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 using autokrige() function automap package in r generate interpolated values grid have defined using 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, using 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 work. now, extract predictions specified cells of interpolated grid using over() function in sp package. tried accomplish using 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 code yielded matrix full of na values. may have happened because kr.grid$krige_output spatialpointsdataframe rather spatialpixelsdataframe or spatialgriddataframe.

does have suggestions how address problem? easy converting kr.grid$krige_output spatialpixelsdataframe or spatialgriddataframe. unfortunately, can't figure out how that.

by creating sequences of x , y coordinates in geometric coordinates , subsequently projecting, grid no longer evenly spaced. try projecting ok_state extent first, , create sequence of evenly-spaced points in 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 

Comments

Popular posts from this blog

android - Get AccessToken using signpost OAuth without opening a browser (Two legged Oauth) -

org.mockito.exceptions.misusing.InvalidUseOfMatchersException: mockito -

google shop client API returns 400 bad request error while adding an item -