Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- library(gstat)
- library(sp)
- library(raster)
- # create a regular grid
- nx=100 # number of columns
- ny=100 # number of rows
- srgr <- expand.grid(1:ny, nx:1)
- names(srgr) <- c('x','y')
- gridded(srgr)<-~x+y
- # generate a spatial process (unconditional simulation)
- g<-gstat(formula=z~x+y, locations=~x+y, dummy=T, beta=15, model=vgm(psill=3, range=10, nugget=0,model='Exp'), nmax=20)
- sim <- predict(g, newdata=srgr, nsim=1)
- r<-raster(sim)
- # generate sample data (Poisson process)
- library(spatstat)
- int<-0.02
- rpp<-rpoispp(int,win=owin(c(0,nx),c(0,ny)))
- df<-as.data.frame(rpp)
- coordinates(df)<-~x+y
- # assign raster values to sample data
- dfpp <-raster::extract(r,df,df=TRUE)
- smp<-cbind(coordinates(df),dfpp)
- smp<-smp[complete.cases(smp), ]
- coordinates(smp)<-~x+y
- # fit variogram to sample data
- vs <- variogram(sim1~1, data=smp)
- m <- fit.variogram(vs, vgm("Exp"))
- plot(vs, model = m)
- # generate 2 conditional simulations with one core processor
- one <- krige(formula = sim1~1, locations = smp, newdata = srgr, model = m,nmax=12,nsim=2)
- # plot simulation 1 and 2: statistics (min, max) are ok, simulations are also ok.
- spplot(one["sim1"], main = "conditional simulation")
- spplot(one["sim2"], main = "conditional simulation")
- # generate 2 conditional with parallel processing
- library("parallel")
- no_cores<-detectCores()
- cl<-makeCluster(no_cores)
- parts <- split(x = 1:length(srgr), f = 1:no_cores)
- clusterExport(cl = cl, varlist = c("smp", "srgr", "parts","m"), envir = .GlobalEnv)
- clusterEvalQ(cl = cl, expr = c(library('sp'), library('gstat')))
- par <- parLapply(cl = cl, X = 1:no_cores, fun = function(x) krige(formula=sim1~1, locations=smp, model=m, newdata=srgr[parts[[x]],], nmax=12, nsim=2))
- stopCluster(cl)
- # merge all parts
- mergep <- maptools::spRbind(par[[1]], par[[2]])
- mergep <- maptools::spRbind(mergep, par[[3]])
- mergep <- maptools::spRbind(mergep, par[[4]])
- # create SpatialPixelsDataFrame from mergep
- mergep <- SpatialPixelsDataFrame(points = mergep, data = mergep@data)
- # plot mergep: statistics (min, max) are ok, but simulated maps show "vertical lines". i don't understand why.
- spplot(mergep[1], main = "conditional simulation")
- spplot(mergep[2], main = "conditional simulation")
Add Comment
Please, Sign In to add comment