Diann,

A good way to start would have been to provide a function that works for
one cell, I had to do some guessing now, but I think the below should point
into the direction of a solution.

library(raster)
library(rgdal)
library(mc2d) #load package

# you can directly read the ArcGRID files using the 'raster' function (and
the rgdal package)
# no need to export them to ascii; but here I use your matrices
Pmn <- raster(matrix(c(3.5, 2.1, 2.9, 3.2, 3.5, 2.1, 2.9, 3.2), nrow=2))
Psd <- raster(matrix(c(0.108, 0.117, 0.092, 0.125, 0.108, 0.117, 0.092,
0.125), nrow=2))
Wab <- raster(matrix(c(47, 33, 56, 44, 47, 33, 56, 44), nrow=2))

# combine the RasterLayer objects into a RasterStack
s <- stack(Pmn, Psd, Wab)

# Function that works for a vector of three numbers
ndunc(1001)
fun <- function(x) {
dPmn <- mcdata(x[1], type="0")
dPsd <- mcdata(x[2], type="0")
dWab <- mcdata(x[3], type="0")
stPte <- mcstoc(rnorm, type="U", dPmn, dPsd, rtrunc=TRUE, linf=0)
stWab <- mcstoc(rtriang, type="U", min=0, mode=dWab, max=75)
risk <- stPte * stWab
# now what output do you want?
# there are many options, e.g.:
# quantile(risk[])
        # or simply:
mean(risk)
}

x <- calc(s, fun)
plot(x)
x <- writeRaster(x, 'mc.tif')  # for ArcGIS
as.matrix(x)


# this is equivalent to
d <- getValues(s)
z <- apply(d, 1, fun)
matrix(z, nrow=2, byrow=TRUE)

Robert

On Mon, Mar 12, 2012 at 3:57 PM, Diann Prosser <dprosser@usgs.gov> wrote:

> Hi all,
>
> I hope you forgive me for this cross-post. I requested help a few days ago
> on the general R-help board but haven't received a response. I think that a
> solution can most likely be answered by someone familiar with working with
> large spatial datasets.
>
> I have raster datasets exported from ArcGRID format into ASCII text files,
> read successfully into R. They are large - 4848 x 4053 cells. I am running
> monte-carlo simulations using the mc2d package (mcstoc function) and am
> having a difficult time applying (or looping) the simulations for each cell
> in my matrix.
>
> Below is the code I've been using to test the process on reduced (2x2)
> matrices, followed by error messages.
> Below that I post the link to (and complete text) from my initial inquiry
> on
> Saturday.
>
> Thank you kindly for any help you might be able to provide.
>
> ________________________________________________________________________
> #Code - using 2x2 example matrices for each raster:
> #terms in bold are the rasters from which I want to pull data for the
> mcstoc
> function
>
> #Construct reduced (2x2) matrices for testing code
> Pmn<- matrix(c(3.5, 2.1, 2.9, 3.2), nrow=2)
> Psd<- matrix(c(0.108, 0.117, 0.092, 0.125), nrow=2)
> Wab<- matrix(c(47, 33, 56, 44), nrow=2)
> m<- matrix(nrow=2, ncol=2)
>
> library(mc2d) #load package
>
> #create mcnodes for use in mcstoc (since array arguments not permitted as
> determined by Error 1)
> dPmn<- mcdata(Pmn, type="0")
> dPsd<- mcdata(Psd, type="0")
> dWab<- mcdata(Wab, type="0")
>
> ndunc(1001)
> for (i in m) {
> stPte <- mcstoc(rnorm, type="U", *dPmn*, *dPsd*, rtrunc=TRUE, linf=0)
> stWab <- mcstoc(rtriang, type="U", min=0, mode=*dWab*, max=75)
> risk <- stPte * stWab
> }
>
> ____________________________________________________________________________
>
> Error 1:
> Error in function (argsd, typethismc)  :
>  Array prohibited in mcstoc as parameter. Use an mcnode instead (in
> response to this, I have created the  code above, but now get Error 2
> below)
>
> Error 2:
> Error in mcdata(Pmn, type = "0") :
>  The array size is not compatible with the node dimension. Should be of
> dim: 1 1 1
>
>
> ____________________________________________________________________________
> Initial post can be found at:
>
> http://r.789695.n4.nabble.com/Draw-values-from-multiple-data-sets-as-inputs-to-a-Monte-Carlo-function-then-apply-across-entire-matx-td4462874.html
>
> 3/10/2012
> Hi all,
>
> I am trying to implement a Monte-Carlo simulation for each cell in a
> spatial matrix (using mcd2 package) .
> I have figured out how to conduct the simulation using data from a single
> location (where I manually input distribution parameters into the R code),
> but am having trouble (a) adjusting the code to pull input variables from
> my various data sets and then (b) applying the entire process across each
> of the cells of the matrices.
>
> I have been doing a lot of reading about loops (a big no-no?), apply, and
> ddply, but can not quite figure it out.
>
> Here is the situation:
>
> Data:
> I have (for simplicity) 3 spatial raster data sets (each 4848 x 4053
> cells) as ASCII files:
> -Poultry density (mean value in each cell)
> -Poultry density (standard deviation in each cell)
> -Wild bird density (single estimate in each cell)
> I read them into R using read.table. The data look correct:
> Pmn <- read.table("D:/Data/PoultryMeans.txt")
> Psd <- read.table("D:/Data/PoultryStDev.txt")
> Wde <- read.table("D:/Data/WildBirdDensity.txt")
>
> The Model:
> In the Monte-Carlo simulation, Poultry and Wild birds have different
> distributions (normal and triangle, respectively).
> Below are the 2 lines of code that use the mcstoc function to draw the
> samples.
> The values in bold are ones that I would like to draw from the data tables
> I read in above.
> For example, 3.5 would be cell (i,j) in the Poultry MEAN density table;
> 0.108 would be cell (i,j) in the Poultry STDEV table; and 47 the single
> estimate for cell (i,j) of the Wild bird density table.
>
> Poultry <- mcstoc(rnorm, type="U", 3.5, 0.108, rtrunc=TRUE, linf=0)
> Wild <- mcstoc(rtriang, type="U", min=0, mode=47, max=75)
> Risk <- Poultry * Wild        #this is the risk function the MC is applied
> to
>
> Questions:
> 1) How can I edit the Poultry and Wild variables above to read the data
> values directly from the 3 input tables (i.e., replacing 3.5, 0.108, and
> 47 with some variable name for the data table and using a loop?? Or
> somehow use apply or ddply?)
> 2) Have the entire process be run for every cell in the 4848 x 4053
> matrix?
>
> Thank you for any help you can provide to get me moving forward!
> Diann
>
> 3/12/2012
> Dear all,
>
> I have developed a loop but receive an error (error 1 below) because the
> mcstoc function prohibits arrays. I tried creating an mcnode to represent
> the data frame, but received Error 2 below (array size not compatible with
> node dimension).
>
> I'm still at a loss on how to apply the Monte-Carlo simulation for each
> cell
> in a matrix. There is an mcapply function, but I did not have success
> applying it. Any insight would be much appreciated.
>
> Error 1:
> Error in function (argsd, typethismc)  :
>  Array prohibited in mcstoc as parameter. Use an mcnode instead
>
> Error 2:
> Error in mcdata(Pmn, type = "0") :
>  The array size is not compatible with the node dimension. Should be of
> dim: 1 1 1
>
> ________________________________________________________________________
> #Code - using 2x2 example matrices for each raster:
>
> #Construct reduced (2x2) matrices for testing code
> Pmn<- matrix(c(3.5, 2.1, 2.9, 3.2), nrow=2)
> Psd<- matrix(c(0.108, 0.117, 0.092, 0.125), nrow=2)
> Wab<- matrix(c(47, 33, 56, 44), nrow=2)
> m<- matrix(nrow=2, ncol=2)
>
> library(mc2d) #load package
>
> #create mcnodes for use in mcstoc (since array arguments not permitted)
> dPmn<- mcdata(Pmn, type="0")
> dPsd<- mcdata(Psd, type="0")
> dWab<- mcdata(Wab, type="0")
>
> ndunc(1001)
> for (i in m) {
> stPte <- mcstoc(rnorm, type="U", dPmn, dPsd, rtrunc=TRUE, linf=0)
> stWab <- mcstoc(rtriang, type="U", min=0, mode=dWab, max=75)
> risk <- stPte * stWab
> }
>
> ____________________________________________________________________________
>
>
>
> --
> View this message in context:
> http://r-sig-geo.2731867.n2.nabble.com/apply-monte-carlo-simulation-for-each-cell-in-a-matrix-originally-raster-tp7367094p7367094.html
> Sent from the R-sig-geo mailing list archive at Nabble.com.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]

