Open cmarmstrong opened 5 years ago
Because vgmArea loops over source/target features, it's easy to parallelize. I've used this function for the past year or so with no issues, so I thought it's time to share it. It saves me bundles of time.
Basic usage is:
## form := formula object ## src/tgt := Spatial* objects ## vgm := gstat::vgm object ## nnodes := passed to parallel::makeCluster function(form, src, tgt, vgm, nnodes, ...) { v <- gstat::variogram(form, data=src) vgm <- gstat::fit.variogram(v, vgm) krg <- gstat::krige0(form, src, tgt, gstat::parVgmArea, vgm=vgm, nnodes=nnodes, ...) }
and a recreation of the example from demo/a2p.R:
Rprof() # import NC SIDS data: library(sp) library(maptools) fname = system.file("shapes/sids.shp", package="maptools")[1] nc = readShapePoly(fname, proj4string = CRS("+proj=longlat +datum=NAD27 +ellps=clrk66")) # reproject to UTM17, so we can use Euclidian distances: library(rgdal) nc = spTransform(nc, CRS("+proj=utm +zone=17 +datum=WGS84 +ellps=WGS84")) # create a target (newdata) grid, and plot: grd = spsample(nc, "regular", n = 1000) class(grd) plot(nc, axes = TRUE) points(grd, pch = 3) library(gstat) # replace this with devtools::load_all('/path/to/gstat/branch') # area-to-point kriging: kr = krige0(SID74 ~ 1, nc, grd, parVgmArea, ndiscr = 9, vgm = vgm(1, "Exp", 1e5, 0), # point variogram, nnodes=4) out = SpatialPixelsDataFrame(grd, data.frame(pred = kr)) pl0 = spplot(nc["SID74"], main = "areas") pl1 = spplot(out, sp.layout = list("sp.polygons", nc, first=F,col='grey'), main = "points on a grid") print(pl0, split = c(1,1,1,2), more = TRUE) print(pl1, split = c(1,2,1,2), more = FALSE)
Because vgmArea loops over source/target features, it's easy to parallelize. I've used this function for the past year or so with no issues, so I thought it's time to share it. It saves me bundles of time.
Basic usage is:
and a recreation of the example from demo/a2p.R: