[R-sig-Geo] Problem with Sequential Gaussian Simulation

[Micha Silver]

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