Frequency of occurrence of dry and wet months rasterbrick R

Question

I have a rasterbrick, a made of monthly time series from 1950-2014 (see below). The task I am trying to accomplish on a is to calculate at each grid point:

negative/positive number of dry/wet months per total months of dry/wet events*100,

where

dry has values in the range −1.49 to −1.00 and wet values are in the range 1.49 to 1.00.

The resulting output is a single rasterLayer with dry=negative percentages and wet= positive percentages such that I can map a raster with blue for wet and red for dry.

Sample data can be found HERE

dd=spei03_df
dd[1:2]<-dd[2:1]#swap lat and lon
a=rasterFromXYZ(dd)
dates=seq(as.Date("1950-01-01"), as.Date("2014-12-31"), by="month")
a=setZ(a,dates)

Answer 1

Here's a solution that should work for you. It does not generate positive and negative percentages, but rather the percentage of one of your values. The way percentages work, you'll still be able to map out what you want.

Note: rather than using your dataset, I created a sample data with the same value range. The values will be more or less random sampled (I included a bias for warmer* to get a tendency for one color), so don't expect a pretty pattern.

library(raster)

## create sample data

#values
dry <- seq(-1,-1.49,by=-0.01)
wet <- dry * -1

#raster brick

r <- lapply(seq_len(length(1950:2014)*12),function(x){

  r <- raster()
  r[] <- sample(c(wet,dry),ncell(r),T,c(rep(0.3,50),rep(0.7,50)))
  r
})

r <- do.call(brick,r)

names(r) <- seq(as.Date("1950-01-01"), as.Date("2014-12-31"), by="month")

#output

> r
class       : RasterBrick 
dimensions  : 180, 360, 64800, 780  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 
data source : /tmp/Rtmp3BTUH1/raster/r_tmp_2017-08-09_120141_2859_73743.grd 
names       : X1950.01.01, X1950.02.01, X1950.03.01, X1950.04.01, X1950.05.01, X1950.06.01, X1950.07.01, X1950.08.01, X1950.09.01, X1950.10.01, X1950.11.01, X1950.12.01, X1951.01.01, X1951.02.01, X1951.03.01, ... 
min values  :       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49,       -1.49, ... 
max values  :        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49,        1.49, ... 
  -1.49,   -1.49,    -1.49,    -1.49,    -1.49,    -1.49,    -1.49,    -1.49, ... 
max values  :    1.49,    1.49,    1.49,    1.49,    1.49,    1.49,    1.49,    1.49,    1.49,     1.49,     1.49,     1.49,     1.49,     1.49,     1.49, ... 

And now to calculate the frequencies:

rfreq <- calc(r,fun = function(x){(sum(x > 0) / nlayers(r))*100})

Like this, every pixel will have a value between 0 and 100, indicating a range from all months dry to all months wet, with the equilibrium at 50. So you can easily map this with red and blue. If you still want the percentages, you can just add a legend with the respective labels:

library(rasterVis)

my.at <- c(0,30,50,70,100)

myColorkey <- list(at=my.at,
                   space = 'bottom',
                   labels=list(
                     labels=c('extreme dry','dry','normal','wet','extreme wet'), ## labels
                     at=c(15,40,50,60,85) ## where to print labels
                   ))

levelplot(rfreq, at=my.at,
          colorkey=myColorkey,scales=list(draw=FALSE),margin=FALSE,
          par.settings=RdBuTheme(),
          main = 'Precipitation 1950 - 2014',
          xlab=NULL, ylab=NULL)

HTH

*

Climate change is real