Skip to content
Prev 27049 / 29559 Next

Problem with Sequential Gaussian Simulation

I've run into a problem trying to produce a number of SGS realizations 
using gstat::krige. (this was definitely working two weeks ago, and now 
I'm getting a seg fault, and I can't figure out why. Any help would be 
appreciated.


Here are the details:


My input "locations":

-------------------------

 > str(gauges_spdf)

Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
 ? ..@ data?????? :'data.frame':??? 24 obs. of? 1 variable:
 ? .. ..$ precip: num [1:24] 13.21 2.66 12.67 8.97 18.69 ...
 ? ..@ coords.nrs : num(0)
 ? ..@ coords???? : num [1:24, 1:2] -37.5 11.5 -9.5 20.5 -27.5 -3.5 
-12.5 33.5 -4.5 34.5 ...
 ? .. ..- attr(*, "dimnames")=List of 2
 ? .. .. ..$ : NULL
 ? .. .. ..$ : chr [1:2] "x" "y"
 ? ..@ bbox?????? : num [1:2, 1:2] -45.5 -49.5 46.5 49.5
 ? .. ..- attr(*, "dimnames")=List of 2
 ? .. .. ..$ : chr [1:2] "x" "y"
 ? .. .. ..$ : chr [1:2] "min" "max"
 ? ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
 ? .. .. ..@ projargs: chr NA


My target "newdata" grid

-------------------------

 > str(grd)
Formal class 'SpatialPixels' [package "sp"] with 5 slots
 ? ..@ grid?????? :Formal class 'GridTopology' [package "sp"] with 3 slots
 ? .. .. ..@ cellcentre.offset: Named num [1:2] -49.5 -49.5
 ? .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
 ? .. .. ..@ cellsize???????? : Named num [1:2] 1 1
 ? .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
 ? .. .. ..@ cells.dim??????? : Named int [1:2] 100 100
 ? .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
 ? ..@ grid.index : int [1:10000] 1 2 3 4 5 6 7 8 9 10 ...
 ? ..@ coords???? : num [1:10000, 1:2] -49.5 -48.5 -47.5 -46.5 -45.5 
-44.5 -43.5 -42.5 -41.5 -40.5 ...
 ? .. ..- attr(*, "dimnames")=List of 2
 ? .. .. ..$ : chr [1:10000] "1" "2" "3" "4" ...
 ? .. .. ..$ : chr [1:2] "x" "y"
 ? ..@ bbox?????? : num [1:2, 1:2] -50 -50 50 50
 ? .. ..- attr(*, "dimnames")=List of 2
 ? .. .. ..$ : chr [1:2] "x" "y"
 ? .. .. ..$ : chr [1:2] "min" "max"
 ? ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
 ? .. .. ..@ projargs: chr NA


My variogram

-------------------------

 > precip_fit
 ? model????? psill??? range
1?? Nug? 0.6391787? 0.00000
2?? Exp 12.0158808 17.39194



The command that I am running

(any value of nsim causes the crash below. But oridinary kriging 
*without* the nsim param works fine)

-------------------------

precip_SGS = krige(formula = precip~1,

 ???????????????????? locations = gauges_spdf,
 ???????????????????? newdata = grd,
 ???????????????????? model = precip_fit,
 ???????????????????? nsim = num_sims,
 ???????????????????? nmax = 4)



The error output:
-------------------------

drawing 10 GLS realisations of beta...

[using conditional Gaussian simulation]

 ?*** caught segfault ***
address (nil), cause 'memory not mapped'

Traceback:
 ?1: predict.gstat(g, newdata = newdata, block = block, nsim = nsim,???? 
indicators = indicators, na.action = na.action, debug.level = debug.level)
 ?2: predict(g, newdata = newdata, block = block, nsim = nsim, 
indicators = indicators,???? na.action = na.action, debug.level = 
debug.level)
 ?3: .local(formula, locations, ...)
 ?4: krige(formula = precip ~ 1, locations = gauges_spdf, newdata = 
grd,???? model = precip_fit, nsim = num_sims, nmax = 4)
 ?5: krige(formula = precip ~ 1, locations = gauges_spdf, newdata = 
grd,???? model = precip_fit, nsim = num_sims, nmax = 4)
 ?6: CreateRainRealizations(gauges_spdf, num_sims)
 ?7: eval(ei, envir)
 ?8: eval(ei, envir)
 ?9: withVisible(eval(ei, envir))
10: source("~/Studies/Research/synthetic/code/run_synthetic.R")


My seesion Info:
-------------------------

 > sessionInfo()

R version 3.5.2 (2018-12-20)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux buster/sid

Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3

locale:
 ?[1] LC_CTYPE=en_US.UTF-8?????? LC_NUMERIC=C
 ?[3] LC_TIME=en_US.UTF-8??????? LC_COLLATE=en_US.UTF-8
 ?[5] LC_MONETARY=en_US.UTF-8??? LC_MESSAGES=en_US.UTF-8
 ?[7] LC_PAPER=en_US.UTF-8?????? LC_NAME=C
 ?[9] LC_ADDRESS=C?????????????? LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats???? graphics? grDevices utils???? datasets methods?? base

loaded via a namespace (and not attached):
[1] compiler_3.5.2


Many thanks,

Micha