CODEWITHSUNDEEP
rmapstime-seriesraster
code123
code123

Reputation: 2146

Convert natural year to hydrological year time series Rasterbrick R

Given a rasterbrick with natural years (Jan-Dec), how can one convert it to have hydrological years (starting on the 1st October and ending on the following 30th September)? Something similar (but not exact) to the question here.

Sample data:

     nl <- 768
     s <- brick(nrows = 510,  ncols = 1068,
                   xmn = -180, xmx = 180, ymn = -90, ymx = 90,
                   crs = "+proj=longlat +datum=WGS84",
                   nl = nl,values=TRUE)
        dates <- seq(as.Date("1950-01-01"), as.Date("2013-12-31"), by = "month")
        s <- setZ(s, dates)
vals <- 1:ncell(s)
s[]=vals

Upvotes: 1

Views: 263

Answers (1)

mikeck
mikeck

Reputation: 3788

I wrote a function for this a while ago.

to_wateryear <- function(thedate){
  thedate = as.Date(thedate)
  year = as.integer(substr(thedate, 1, 4))
  cutoff =  cutoff = as.Date(paste(year, '-10-01', sep=''))
  return(ifelse(thedate < cutoff, 
           year*1000 + thedate - as.Date(paste(year - 1, '-09-30', sep='')),
           (year + 1)*1000 + thedate - as.Date(paste(year, '-09-30', sep=''))))
}

Examples:

to_wateryear(as.Date(Sys.time()))
to_wateryear('2013-10-01')
# leap year behavior
to_wateryear('2012-09-30')
to_wateryear('2013-09-30')
to_wateryear('2012-02-29')
to_wateryear('2013-03-31')

EDIT From your comments it sounds like what you actually want is to split your rasterbrick by water year, rather than representing dates in terms of hydrologic year.

library(raster)
nl <- 768
s <- brick(nrows = 510,  ncols = 1068, 
  xmn = -180, xmx = 180, ymn = -90, ymx = 90,
  crs = "+proj=longlat +datum=WGS84",
  nl = nl,values=TRUE)

dates <- seq(as.Date("1950-01-01"), as.Date("2013-12-31"), by = "month")
# use water year
s <- setZ(s, to_wateryear(dates))
vals <- 1:ncell(s)
s[]=vals

# split by water year
hyears = unique(getZ(s) %/% 1000)
res = vector("list", length(hyears))
for(i in 1:length(hyears))
  res[[i]] = s[[which(getZ(s) %/% 1000 == hyears[i])]]
names(res) = paste0("HY", hyears)

Upvotes: 1

Related Questions