When dealing with raster images, I found the raster and rts packages in R to be very flexible for analysing my raster data unlike ArcGIS. In my learning efforts, I have had to do some stacking/bricking of individual images to form a continuous time series.
Now, I would like to extract some values of interest from the stacks just as one will do with a time series object in R.
My raster data (A
) is rainfall every 30mins. So layers in my rasterstack are separated by
30minstime step. For comparison with another raster data (
B) set which is at 6-hourly time step, I summed
A` to 6-hourly as follows:
dates1=seq(as.POSIXct("2015-04-01 00:00:00"), as.POSIXct("2015-11-30 23:59:59"), by="6 hours",tz="GMT")
index=rep(seq(1,976,by=1),each=6) #11712/12=976 nlayers/# of intervals for accumulation
ras <- setZ(A,dates1)
sixhrly <- stackApply(ras,indices=index, fun=sum)
Output from a pixel of sixhrly looks like:
`time` `value`
2015-11-28 00:00:00 5.452091e-02
2015-11-28 06:00:00 4.030154e-01
2015-11-28 12:00:00 3.107421e-01
2015-11-28 18:00:00 8.851760e-02
2015-11-29 00:00:00 9.171816e-02
2015-11-29 06:00:00 8.803389e-02
2015-11-29 12:00:00 3.286928e-02
2015-11-29 18:00:00 0.000000e+00
2015-11-30 00:00:00 5.635915e-02
2015-11-30 06:00:00 1.299293e-03
2015-11-30 12:00:00 3.354174e-02
2015-11-30 18:00:00 6.680167e-03
1) using stackApply
(which I find to have the most speed on my system), how can I find the maximum value per pixel month? The output is a rasterstack
, with layer1=April
, layer2=May
etc...
2) How can I obtain the timestamps for each max value? This is important if I want to know what time the max rainfall occurred and if A
and B
.
A more generic way for doing this is most welcome as well since data can be on hourly and minutes time step over many years and one has to extract some values of interest.