On Thu, Nov 24, 2011 at 2:31 PM, VanWilgenburg,Steve [Sas] <
Steve.VanWilgenburg@ec.gc.ca> wrote:

> Hi Robert
>
> I only intended that spreadsheet as another way illustrate what it was I
> was trying to accomplish. I have attached two ascii rasters and data and
> functions in another file if that helps.
>
>
What would help me most is a reproducible example in the message (if you
can not do it with Raster objects, use a data.frame or matrix) that shows
(in code, or pseudo-code) what you want to do. As you did in part in your
first message, and see below.



> I am still trying to figure this out, sorry if I am rather slow on the
> uptake...
>
>
> If x contains the data values I am assessing and r is the likelihood
> function applied to the map, should this
> fun <- function(x) cumsum(sort(x)
> ss <- calc(s, fun)
>
> Instead be...
> fun <- function(r) cumsum(sort(r)
> ss <- calc(s, fun)
>
>
no, 'x' is a local variable, that refers to the layers in s

since I want the cells in s[i] to be a cumulative sum of the cells in
> r[i]?
>
>
yes, but s is the stack of layers in r

s <- stack(r)



> When I tried this on some simulated data, it does not appear to work, as
>
> > ss <- calc(r[[1]], fun)
> > ss
> class       : RasterLayer
> dimensions  : 20, 20, 400  (nrow, ncol, ncell)
> resolution  : 2.5, 2.5  (x, y)
> extent      : -5, 45, 25, 75  (xmin, xmax, ymin, ymax)
> projection  : +proj=longlat +datum=WGS84
> values      : in memory
> min value   : 2.425674e-14
> max value   : 0.009195606
>
> Has a Maximum value that does not equal 1 (as it should), and in fact is
>
exactly equal to the likelihood surface it was supposed to have
> calculated the cumulative sum on i.e.
>


It should not be 1 because it is a cumsum for a singe layer.

cumsum(0.5) == 0.5


>
> > r[[1]]
> class       : RasterLayer
> dimensions  : 20, 20, 400  (nrow, ncol, ncell)
> resolution  : 2.5, 2.5  (x, y)
> extent      : -5, 45, 25, 75  (xmin, xmax, ymin, ymax)
> projection  : +proj=longlat +datum=WGS84
> values      : in memory
> min value   : 2.425674e-14
> max value   : 0.009195606
>
>
> As opposed to what happens with a vector of simulated data
> > p<-dnorm(x,mean=-74,sd=15)/sum(p)
> > p
>  [1] 0.04198945 0.04163704 0.07339848 0.05509149 0.18052523 0.17121966
> 0.11465071 0.17160485 0.06263885 0.08724424
> > cumx<-cumsum(p)
> > cump<-cumsum(sort(p))
> > cump
>  [1] 0.04163704 0.08362649 0.13871798 0.20135684 0.27475531 0.36199956
> 0.47665026 0.64786992 0.81947477 1.00000000
>
>
Yes, but here you use a vector (=multiple layers, and you also divide by
sum(p)
Building on your example:


library(raster)
map <- raster(nrows=20, ncols=20, xmn=-5, xmx=45, ymn=25, ymx=75)
map[] <-rnorm(ncell(map), -89, 15)
x <- rnorm(10,mean=-74, sd= 15)
SD <-11.67114
r <- list()
for (i in 1:length(x)){
r[i] = (1/(sqrt(2*pi)*SD)*exp(-(1/(2*(SD)^2))*(x[i]-map)^2))
    r[i] <- r[[i]] / cellStats(r[[i]],'sum')
}
s <- stack(r)
fun <- function(x) cumsum(sort(x)/sum(x))
ss <- calc(s, fun)

cell1 <- s[1]
fun(cell1)

this is the same as

ss[1]


is this what you are after:

fun2 <- function(x) {
a <- cumsum(sort(x)/sum(x))
a <- splinefun(1:nl, a)(1:(10*nl)/10)
which.max(a > threshold) / (10 * nl)
}

nl <- nlayers(s)
threshold = 0.67
ss2 <- calc(s, fun2)




In addition, ideally if we can get a cumsum function to work on the
> raster objects, then I would like to apply it within a loop...
>
> > z<-list()
> > for (i in 1:length(x)){
> + z[i] <- calc(r[[i]], fun)
> + }
>
>
I doubt that. Rather create a stack and compute on that.

you said:

> 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

So if the above is not correct, what would you do with data.frame d,
created below?

library(raster)
map <- raster(nrows=20, ncols=20, xmn=-5, xmx=45, ymn=25, ymx=75)
map[] <-rnorm(ncell(map), -89, 15)
plot(map)
x <- rnorm(10,mean=-74, sd= 15)
SD <-11.67114
r <- list()
for (i in 1:length(x)){
r[i] = (1/(sqrt(2*pi)*SD)*exp(-(1/(2*(SD)^2))*(x[i]-map)^2))
    r[i] <- r[[i]] / cellStats(r[[i]],'sum')
}
s <- stack(r)
d <- data.frame(values(s))

If you can get this to work:

fun <- function(x) {
... ?
... ?
}

res <- fun(d)


Then the next step is to do (perhaps after a few changes)

x <- calc(s, fun)


Hope this gets you closer,
Robert



>
> -Steve
>
> ________________________________
>
> From: Robert J. Hijmans [mailto:r.hijmans@gmail.com]
> Sent: Thursday, November 24, 2011 1:38 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,
>
> 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
>
>
>
>

	[[alternative HTML version deleted]]

