Steve,

The way to start is to write a function that does what you do on a vector
(cell) or matrix (rows are cells). It should then be possible to call that
that function with the 'calc' function.

For example, to get the cumsum of sorted values, you can do

fun <- function(x) cumsum(sort(x)
ss <- calc(s, fun)

The above works because calc finds out that it needs to use 'fun' within a
call to 'apply'.
You can also take the values from the RasterStack (if it is not too large
for memory); or to test with a subset

v <- getValues(s)

do the computations

assign back
# single layer result
r <- raster(s)
values(r) <- v

# multiple layers result
b <- brick(s, values=F)
values(b) <- v

Sometimes this takes some strange tweaks to get it right, but it would be
much easier for me to advice based on an R function with example data than
first having to translate from Excel.

Robert



On Thu, Nov 24, 2011 at 11:10 AM, VanWilgenburg,Steve [Sas] <
Steve.VanWilgenburg@ec.gc.ca> wrote:

> **
>
> Hi Robert
>
> Thanks for the attempt, however your suggestions do not appear to be doing
> what I want.
>
>
> In essence, what I am trying to accomplish is to select a 'Confidence
> region" within the map, for each sample in x, by first assessing the
> likelihood given the data, that the value x[i] could have been generated by
> a process at given location (cell) in the map. I have that part figured out
> (assessing likelihoods). But in order to select a confidence region, I
> essentially need to order the likelihoods from smallest to largest, then
> calculate the cumulative sum of those values. Then using the spline
> function, find the likelihoods at which a threshold cumulative probability
> would be reached to find the Upper portion (most likely) of the
> distribution.
>
>
> Basically I want to select a region of the map (groups of cells (not
> necessarily contiguous) ) that if used in an overlay, would define a region
> such that the selected cells in the likelihood surface (r[i]) would sum to
> the specified amount (e.g. 67%) of HIGHEST probabilities
>
>
> For each r[i], after the following step...
> for (i in 1:length(x)){
> r[i] <- r[[i]] / cellStats(r[[i]],'sum')
> }
>
> I want to assess Cumulative probabilities ( i.e. run a function like
> 'cumsum' after ordering the values of r[i] from smallest to largest). I
> have attached an Excel spreadsheet the illustrates what I want to do, with
> rows in the sheets representing a cell in the raster. After defining a
> confidence region for each sample, I then want to assess the cumulative
> number of samples that were consistent with the data being potentially
> generated at a location in the map, hence the following step..
>
> s<-stack(r)
> i<-maxValue
> rc <- overlay(s, fun=function(i){return(sum(i))} )
>
> For example, in the attached spreadsheet (see sheet r1), the cells with
> the Normalized (to sum) likelihoods >= ~0.0566 would yield the cells
> associated with the upper 67% of probability distribution. Thu the spline
> function would be similar to how I used the polynomial equations in the
> attached to solve for the predicted value of the normalized likelihood when
> (in this case) the cumulative probability = 0.33. The column "likely versus
> unlikely" then reclassifies the "raster cells" based on that predicted
> value.
>
> <<Example of selecting upper Cumulative Probailities.xls>>
> Thus, using the examples from the spreadsheet, the reclassifications would
> be something approximately as follows:
> r[[1]] <-reclass(r[[1]], c(-Inf, 0.0566, 0, 0.0566, Inf, 1))
> r[[2]] <-reclass(r[[2]], c(-Inf, 0.091, 0, 0.091, Inf, 1))
> r[[3]] <-reclass(r[[3]], c(-Inf, 0.065, 0, 0.065, Inf, 1))
> r[[4]] <-reclass(r[[4]], c(-Inf, 0.0674, 0, 0.0674, Inf, 1))
> r[[5]] <-reclass(r[[5]], c(-Inf, 0.15, 0, 0.15, Inf, 1))
> Etc…..
>
> I  suspect it would be difficult (impossible?) to sort the rasters r[i] ?
> Hopefully I am mistaken. If not, would there be way use a loop to extract
> each r[i] to a series of data frames held in memory which could then be
> sorted and cumsum computed and spline function fit to those values? I tried
> creating a loop…
>
> cells<-ncell(r[[1]])
>
> exp<-list()
> for (i in 1:length(x)){
> exp[i] <- cellValues(r[[i]],1:cells)
> }
>
> However, rather then exp[[1]] having 400 values as I expected, it only
> gives one value. Alternatively, perhaps exporting values from the map
> directly to a dataframe then using the loops to calculate a new column(s)
> for each x[i]?
>
> I also tried using extract, but it does not appear to work either.
>
> I hope this clarifies what I am attempting to do.
>
>
> -Steve
>
>
>
> ________________________________
>
> From: Robert J. Hijmans [*mailto:r.hijmans@gmail.com*<r.hijmans@gmail.com>]
>
> Sent: Wednesday, November 23, 2011 6:58 PM
> To: VanWilgenburg,Steve [Sas]
> Cc: r-sig-geo@r-project.org
> Subject: Re: [R-sig-Geo] Finding a way to select upper portion of
> cumulative probabilities in a map of probability estimates
>
> Steve, some suggestions below, Robert
>
> On Wed, Nov 23, 2011 at 12:28 PM, VanWilgenburg,Steve [Sas] <
> Steve.VanWilgenburg@ec.gc.ca> wrote:
>
>         Hi
>
>         I have a problem that I am hoping someone can help me with.
>
>
>         I have been using the raster package and some loops to use a
> normal
>         probability density function to assess likelihoods that a sample
> in a
>         data set could have come from a location on a map (example script
>         below).
>
>         library(raster)
>         ##### example map
>         map <- raster(nrows=20, ncols=20)
>         map[] <-rnorm(ncell(map), -89, 15)
>         xmin(map)<- -5
>         xmax(map)<- 45
>         ymin(map)<- 25
>         ymax(map)<- 75
>
>         plot(map)
>
>         ###### values (data) to be assess likelihoods for
>         x <- rnorm(10,mean=-74, sd= 15)
>
>         ################ Define the variation
>         SD <-11.67114
>
>         ######establish a set of lists to hold the results of subsequent
> loops
>         r <- list()
>         ############ Apply Normal Probability Density Function to Assess
>         Likelihoods ####################
>         for (i in 1:length(x)){
>            r[i] = (1/(sqrt(2*pi)*SD)*exp(-(1/(2*(SD)^2))*(x[i]-map)^2))
>         }
>
>         ############ Normalize by Sum of Values....yielding map of
> probabilities
>         that sums to 1 ####################
>         for (i in 1:length(x)){
>         r[i] <- r[[i]] / cellStats(r[[i]],'sum')
>         }
>
>         ##example plot
>         plot(r[[1]])
>
>
>         This method is working beautifully to assess the likelihoods, but
> I now
>         want to come up with a way of accurately selecting and
> reclassifying the
>         resulting probability maps cells into binary likely (1) versus
> unlikely
>         (0) maps based on selecting all cells in the map above a threshold
>         probability value, such that the sum of the probabilities in the
>         "likely" category is equal to a specified cumulative probability
> (e.g.
>         67%, 95%). I have previously been using some course
> approximations, but
>         now want a highly precise., accurate and efficient method to do
> this,
>         but appear to be stuck.
>
>         If this were a data frame, I could just sort based on the
> estimated
>         probabilities and calculate a new cumulative probability
> variable(s). I
>         could probably then use a spline function to estimate the minimum
>         probability at which the selected cumulative probability of
> interest was
>         reached (similar to the following post
>         (*http://finzi.psych.upenn.edu/R/Rhelp01/archive/13524.html))*<http://finzi.psych.upenn.edu/R/Rhelp01/archive/13524.html))>.
> However, I
>         am not certain how I can do this for raster objects.
>
>         Ideally I would like to apply a spline function to find the
> desired
>         threshold value and directly incorporate the result into a loop
> that
>         would reclass each r[i] into a binary map,
>         e.g. r[[1]] <-reclass(r[[1]], c(-Inf, 0.000694, 0, 0.000694, Inf,
> 1))
>         where 0.000694 was pasted from the lookup spline function
>
>
> It is not entirely clear to me what you are attempting to do, but I am
> under the impression that you are looking for something like this:
>
> s <- stack(r)
> nl <- nlayers(s)
>
> sp <- calc(s, function(x)splinefun(1:nl, x)(1:nl))
>
> or
>
> sp <- calc(s, function(x)splinefun(1:nl, x)(1:nl) > x)
> plot(sp)
>
>
>
>         Then I want to sum the "likely" areas across all samples in the
> data set
>
>
>         E.g.
>         s<-stack(r)
>         i<-maxValue
>         rc <- overlay(s, fun=function(i){return(sum(i))} )
>
>
>
>
>
>
> sp <- calc(s, function(x) sum( splinefun(1:nl, x)(1:nl) > x) )
>
>  ?
>
>
>
>         I would really appreciate it if anyone can provide advice on
>         A) how to estimate the threshold value for each sample
>         B) how I could add a loop and paste the result from A into a list
> of
>         reclass definitions (a separate one for each sample)
>
>         I am hoping to find the most efficient method possible to
> accomplish
>         this, especially since x for some of my applications can be very
> long
>         (e.g. 1000 samples).
>
>         -Steve
>
>
>
>
>
>
>                [[alternative HTML version deleted]]
>
>         _______________________________________________
>         R-sig-Geo mailing list
>         R-sig-Geo@r-project.org
>         *https://stat.ethz.ch/mailman/listinfo/r-sig-geo*<https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>
>
>

	[[alternative HTML version deleted]]

