

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


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://rsiggeo.2731867.n2.nabble.com/questionsonRasterStackBricktp5553580p5553580.html> Sent from the Rsiggeo mailing list archive at Nabble.com.
>
> _______________________________________________
> RsigGeo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/rsiggeo [[alternative HTML version deleted]]
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo


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 RForge 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
_______________________________________________
RsigGeo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/rsiggeo


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


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 coworked 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


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 (zdimension), 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


> 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


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

