Cleans the Berkeley Earth BEST temperature anomaly data and generates a raster of dated, temperature estimates

clean_BEST_data.r

returnBESTtemp <- function(tempFile) {
  
  # tempFile: complete path to a netcdf BEST temperature file of daily records [need not be geographic subset]
  
  ave <- brick(tempFile, var = "climatology")
  anomalies <- brick(tempFile, var = "temperature")
  
  # -- deal with absence of leap year in anomaly series
  ave.dates <- as_date(as.numeric(unlist(lapply(strsplit(names(ave), "X"), function(x) x[[2]]))), origin = "1949-12-31")
  which(ave.dates == "1950-02-28")   # 59th doy
  
  leap.ave <- stack(ave[[1:59]], ave[[59]], ave[[60:365]])   # leap year version duplicates climatology from feb 28 for leap day
  leap.ave.bd <- format(as_date(1:366, origin = "1951-12-31"), format = "%b-%d")
  
  clima <- stack(ave, ave, leap.ave, ave)   # since starting with 1950, non-leap, non-leap, leap, non-leap, and then recycle. 
  clima.fullsize <- stack(replicate(17, clima))  # 17 iterations of this 4 year sequence, now comparable in dimensions to anomalies object  
  
  dates <- as_date(as.numeric(unlist(lapply(strsplit(names(anomalies), "\\."), function(x) x[[2]]))), origin = "1949-12-31")
  years <- unique(year(dates)) 
  month.day <- format(dates, format="%b-%d")
  
  # -- drop 2018, which is problematic because it's only a partial year
  anomalies <- subset(anomalies, which(!year(dates) %in% 2018), drop = T) 
  dates <- dates[!year(dates) %in% 2018]
  
  print(length(dates) == dim(anomalies)[3])
  temp <- anomalies + clima  # recycles to length of anomalies 
  
  # -- perform spot checks to ensure values are properly summed 
  test.layers <- floor(runif(20) * dim(anomalies)[3]) %>% unique()
  anomalies.test <- anomalies[[test.layers]]
  monthday.test <- month.day[test.layers]
  indices <- unlist(lapply(monthday.test, function(x) which(leap.ave.bd == x)))   # julian day values 
  
  # -- these all check out -- satisfied it's taking the correct sums -- # 
  sum.random <- anomalies[[test.layers]] + leap.ave[[indices]]  # manual combination of anomalies layer and climatology
  true.values <- temp[[test.layers]]   # our combined temperature values layer 
  min(true.values == sum.random)

  
  temp <- setZ(temp, dates)
  names(temp) <- as.character(dates)
  
  temp
}