questions on RasterStack/Brick

classic Classic list List threaded Threaded
8 messages Options
Reply | Threaded
Open this post in threaded view
|

questions on RasterStack/Brick

Martin
Dear all,

I have some questions concerning stacks and bricks from the raster package.

1. is it possible to flip the whole stack vertically? I get an error when I use the 'flip' function..

2. is it possible to use the calc function (or something else) on a stack to get a raster, which shows the regression slope of the pixels through the stack (treated like a kind of time series)? And is it possible to compute a correlation between two stacks of the same dimensions?

greetings,
Martin

Martin Brandt

Department of Geography and
Regional Research
UZA II, Althanstr. 14
1090 Wien, AUSTRIA
http://geooekologie.univie.ac.at
Reply | Threaded
Open this post in threaded view
|

Re: questions on RasterStack/Brick

Nikhil Kaza-2
I am getting familiar with raster package. so this may not be the most  
elegant way

flipBrick <- function(x, direction="x"){
  br <- brick(raster(extent(x), nrow=nrow(x), ncol=ncol(x),  
crs=projection(x)))
        for(i in 1:nlayers(x)){
                k <- flip(raster(x,i),direction)
                br <- addLayer(br,k)
                }
return(br)
        }



Nikhil Kaza
Asst. Professor,
City and Regional Planning
University of North Carolina

[hidden email]

On Sep 21, 2010, at 3:49 AM, Martin wrote:

>
> Dear all,
>
> I have some questions concerning stacks and bricks from the raster  
> package.
>
> 1. is it possible to flip the whole stack vertically? I get an error  
> when I
> use the 'flip' function..
>
> 2. is it possible to use the calc function (or something else) on a  
> stack to
> get a raster, which shows the regression slope of the pixels through  
> the
> stack (treated like a kind of time series)? And is it possible to  
> compute a
> correlation between two stacks of the same dimensions?
>
> greetings,
> Martin
>
>
> --
> View this message in context: http://r-sig-geo.2731867.n2.nabble.com/questions-on-RasterStack-Brick-tp5553580p5553580.html
> Sent from the R-sig-geo mailing list archive at Nabble.com.
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo


        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: questions on RasterStack/Brick

Robert Hijmans
In reply to this post by Martin
On Tue, Sep 21, 2010 at 12:49 AM, Martin <[hidden email]> wrote:
>
> Dear all,
>
> I have some questions concerning stacks and bricks from the raster package.
>
> 1. is it possible to flip the whole stack vertically? I get an error when I
> use the 'flip' function..
>

You can use what Nikhil send you or perhaps (for objects that are not
too big), with RasterStack s

x <- as(s, 'SpatialGridDataFrame')
x <- flipVertical(x)
s2 <- stack(x)  # or   b <- brick(x)

or get raster 1.5.9 from R-Forge in a couple of hours and then do

x <- flip(s, direction='y')


> 2. is it possible to use the calc function (or something else) on a stack to
> get a raster, which shows the regression slope of the pixels through the
> stack (treated like a kind of time series)?

You can do things like this:

myfun = function(v, ...) { d <- data.frame(x=1:length(v), y=v);
lm(y~x, data=d)$coefficients[2] }

a <- calc(stack, fun=myfun)


> And is it possible to compute a
> correlation between two stacks of the same dimensions?

Not directly, I think, but if you have two stacks, s1, and s2, and
some patience, you can do:

 r <- raster(s1)
for (i in 1:ncell(r)) {
    r[i] <- cor(as.vector(cellValues(s1, i)), as.vector(cellValues(s2, i)))
}


Robert


> greetings,
> Martin
>
>
> --
> View this message in context: http://r-sig-geo.2731867.n2.nabble.com/questions-on-RasterStack-Brick-tp5553580p5553580.html
> Sent from the R-sig-geo mailing list archive at Nabble.com.
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: questions on RasterStack/Brick

Martin
Dear Nikhil and dear Robert,

thanks a lot for your help.
I first tried Nikhils function, it works well and helps me to learn how these functions work.
Good to see that flipping bricks will be implemented in the new raster package. The raster package is a great piece of work. It makes all those expensive GIS systems almost useless for me, it's exactly what i was looking for for so many years..

I also tried the regression and correlation code and it works well. a huge thanks!

Martin
Martin Brandt

Department of Geography and
Regional Research
UZA II, Althanstr. 14
1090 Wien, AUSTRIA
http://geooekologie.univie.ac.at
Reply | Threaded
Open this post in threaded view
|

Re: questions on RasterStack/Brick

Logan
In reply to this post by Robert Hijmans
Dear all,

Robert suggested that, with a little patience, the following code could be used to correlate two raster stacks:

 r <- raster(s1)
for (i in 1:ncell(r)) {
    r[i] <- cor(as.vector(cellValues(s1, i)), as.vector(cellValues(s2, i)))
}

Not having that much patience, my co-worked Pieter Beck and I wrote a function to correlate two raster stacks of same x, y, and z extents. The function is considerably faster than the code above:

stack.correlation <- function(stack1, stack2, cor.method, na.val, na.rm = T){
  # output template
  cor.map <- raster(stack1)
  # combine stacks
  T12 <- stack(stack1,stack2)
  NAvalue(T12) = na.val
 
  # the function takes a vector, partitions it in half, then correlates
  # the two sections, returning the correlation coefficient.
  stack.sequence.cor <- function(myvec,na.rm=T){
    myvecT1<-myvec[1:(length(myvec)/2)]
    myvecT2<-myvec[(length(myvec)/2+1):length(myvec)]
    return(cor(myvecT1,myvecT2, method =  cor.method))
    }

  # apply the function above to each cell and write the correlation
  # coefficient to the output template.
  cor.map <- stackApply(T12, indices = rep(1, length(jul_mckenney)),
    fun = stack.sequence.cor, na.rm = TRUE)
 
  return(cor.map)
}

example call:
my.cor.map <- stack.correlation(tmp.jul, tmp.aug, cor.method = "spearman", na.val = -9999)

Hope this helps,

Logan Berner
lberner@whrc.org
The Woods Hole Research Center
Logan Berner
The Woods Hole Research Center
www.whrc.org
Reply | Threaded
Open this post in threaded view
|

Re: questions on RasterStack/Brick

Logan
Apologies-- there was an error in the stack correlation function of my last post. Below is an updated version of the function.  

stack.correlation <- function(stack1, stack2, cor.method, na.val, na.rm = T){
  # output template
  cor.map <- raster(stack1)
  # combine stacks
  combined.stack <- stack(stack1,stack2)
  NAvalue(combined.stack) = na.val
 
  # The function takes a vector of cell values (z-dimension), partitions it in half, then correlates
  # the two sections, returning the correlation coefficient.
  stack.sequence.cor <- function(full.vec,na.rm=T){
    vec1<-full.vec[1:(length(full.vec)/2)]
    vec2<-full.vec[(length(full.vec)/2+1):length(full.vec)]
    return(cor(vec1,vec2, method=cor.method))
    }

  # Apply the function above to each xy cell and write the correlation
  # coefficient to the output template.
  cor.map <- stackApply(combined.stack, indices = rep(1, length(combined.stack)),
    fun = stack.sequence.cor, na.rm = TRUE)
  return(cor.map)
}

example call:
my.cor.map <- stack.correlation(august_temp, july_temp, "spearman", -9999)
#plot(my.cor.map)


Logan
Logan Berner
The Woods Hole Research Center
www.whrc.org
Reply | Threaded
Open this post in threaded view
|

Re: questions on RasterStack/Brick

Robert Hijmans

> Apologies-- there was an error in the stack correlation function of my last post. Below is an updated
> version of the function.  

Thanks Logan,
That is very useful, but I would not use stackApply, but the simpler calc. Below a reworked version of your function:


stackcor <- function(s1, s2, method='spearman') {
        mycor <- function(v) {
                x <- v[1:split]
                y <- v[(split+1):(2*split)]
                cor(x, y, method=method)
        }
        s <- stack(s1, s2)
        split <- nlayers(s)/2
        calc(s, fun=mycor )
}

Best, Robert
Reply | Threaded
Open this post in threaded view
|

Re: questions on RasterStack/Brick

ppanday
Hi Robert

To create a map of significance of the correlation, I tried substituing the command cor with

cor.test(x,y, method=method)$p.value

However, I get an error asying "Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
  cannot use this function"

Is it possible to use this function to do that so I get a significance of each pixel correlation?

Prajjwal