I am performing conditional Gaussian simulations following the tutorial. Actually, I am re-running the code I wrote last year. My objective is to generate 6000 simulations (could be 5000 or 4000 as well) at 2050 locations. Last year I did consecutively sets of 1000 simulations (and one run of 3000 simulations) because my computer was not very powerful, but now that I have a new computer I was thinking of running 6000 simulations at once. I get the following message:
The function is interrupted
Traceback: simtub(data.db, grid.db, data.model, neigh, seed = seed, nbsimu = nbsimu, nbtuba = nbtuba)
Error: Error in R_simtub
I try the code with just 100 simulations and it works. Is it possible that the message is due to memory problems?
I include the code I am using bellow.
Thanks in advance for your help,
Mercedes
- Code: Select all
# ### Loop for 6 depths -----------------------------------------------
### Different depths
depths <- c("0_5", "5_15", "15_30", "30_60", "60_100", "100_200")
my.seeds <- c(1984, 1991, 1946, 1948, 2812, 2109)
### Create list to store the outputs
simulations.list <- list()
vgm.models.list <- list()
for(i in 1:length(depths)){
print(depths[[i]])
### Load the residuals for this GSM depth
texture <- read.csv(paste0("table.granulo.data_",depths[[i]],".csv"), header=TRUE, sep=";", dec=",",as.is = T)
### Eliminate the first column
texture <- texture[,-1]
mydata <- texture[, c("x", "y", "id_profil",
paste0("resid.cub.alr.limon_",depths[[i]]),
paste0("resid.cub.alr.argile_",depths[[i]]))]
### Create the dataframe in RGeostat format
data.db <- db.create(mydata, flag.grid=FALSE, ndim=2, autoname=F)
### Perform the Gaussian Anamorphis function
data.anam.limon <- anam.fit(data.db,
paste0("resid.cub.alr.limon_",depths[[i]]),
title="Anamorphosis Function")
print(data.anam.limon)
data.anam.argile <- anam.fit(data.db,
paste0("resid.cub.alr.argile_",depths[[i]]),
title="Anamorphosis Function")
print(data.anam.argile)
### Transform raw data into gaussian
# The next operation converts the raw data into their gaussian transform
# Gaussian.limon, using the anamorphosis function data.anam, by typing:
data.db <- anam.z2y(data.db,paste0("resid.cub.alr.limon_",depths[[i]]), anam=data.anam.limon)
data.db <- anam.z2y(data.db,paste0("resid.cub.alr.argile_",depths[[i]]),anam=data.anam.argile)
### Variables of interes
data.db = db.locate(data.db,paste0("Gaussian.resid.cub.alr.limon_",depths[[i]]),"z",1)
data.db = db.locate(data.db,paste0("Gaussian.resid.cub.alr.argile_",depths[[i]]),"z",2)
### Using the default values
data.vario <- vario.calc(data.db) ### Lag is ~ 67260 m
plot(data.vario, npairdw=TRUE, npairpt=1,title="Experimental Default Variogram")
### Fitting the models
data.model <- model.auto(data.vario,
struct=melem.name(c(1,3)),
verbose=1)
vgm.models.list[i] <- data.model
### Load RMQS sites
load("RMQS_sites.RData")
### Select only the locations
RMQSXY <- RMQS_sites[,c('x_reel','y_reel')]
# Set number of dimensions
ndim = 2
# Create a prediction grid from the RMQS coordinates
grid.db = db.create(RMQSXY[,c('x_reel','y_reel')], flag.grid=F, ndim=2, autoname=F)
# Define the neighbourhood conditions
neigh = neigh.create(ndim=ndim,
type=2,
nmini = 10, # Minimum number of observations
nmaxi = 50, # max number of observations
flag.sector=TRUE, # the Moving Neighborhood search is performed by angular sectors
nsect=8, # Search is done in 8 octants
nsmax = 10,
radius = 100000)
### Now we can perform the conditional Gaussian simulations.
# For this we need the original data, the neighbourhood, the grid, and the variogram model
# Set the seed
seed = my.seeds[[i]]
# A set of 3000 simulations is performed using the Turning Bands algorithm (with 1000 bands)
nbtuba = 1000
nbsimu = 3000
grid.db = simtub(data.db, # the data
grid.db, # the grid where we want to predict
data.model, # the variogram model
neigh, # the moving neighbourhood
seed = seed, # the seed number
nbsimu = nbsimu, # number of simulations
nbtuba = nbtuba # number of bands
)
grid.db <- anam.y2z(grid.db,names=paste0("Simu.Gaussian.resid.cub.alr.limon_",depths[[i]],"*"),anam=data.anam.limon)
grid.db <- anam.y2z(grid.db,names=paste0("Simu.Gaussian.resid.cub.alr.argile_",depths[[i]],"*"),anam=data.anam.argile)
simulations.list[i] <- grid.db
}