Skip to content

adds parallelized vgmArea function #49

New issue

Have a question about this project? # for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “#”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? # to your account

Open
wants to merge 1 commit into
base: main
Choose a base branch
from

Conversation

cmarmstrong
Copy link

@cmarmstrong cmarmstrong commented Jun 28, 2019

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)

# for free to join this conversation on GitHub. Already have an account? # to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

1 participant