REGEN Gridded Dataset Documentation Paper

Introduction

Why do we need the dataset?

Data & Methods

The CCRC dataset was created by acquiring daily station precipitation data from various sources. This raw data was then subjected to an automated quality control process. Subsequently the station data from various sources was merged based on an automated station matching algorithm. Finally the final combined station list was interpolated using oridnary block Kriging. In this section the various raw data sources, automated quality control, automated station matching algorithm and the interpolation method are described.

Raw Gauge data

The raw station data for CCRC dataset came from 4 sources:

  1. the Global Historical Climatological Network - Daily stations hosted by National Oceanic and Atmospheric Administration (NOAA) in USA (Menne et al. 2012).
  2. the raw station data hosted my Global Precipitation Climatology Centre (GPCC), part of the Deutsche Wetterdienst (DWD)
  3. raw station data for Argentina (ask Lisa for source) and
  4. raw station data for Russia (ask Lisa for source)
rawdata <- read.table(file = "/srv/ccrc/data11/z3289452/CCRCGlobalGriddedDailyPrecip/Metadata/VERSION1.1_data/station-yearly-summary.txt")

totRecs <- array(dim = c(10,7,68))
totFlgs <- array(dim = c(10,7,68))
for (Reg in 1:10) {
  for (Src in 1:7) {
    totRecs[Reg,Src,] <- as.numeric(rawdata[(10*2)*(Src - 1) + Reg*2 - 1,])
    totFlgs[Reg,Src,] <- as.numeric(rawdata[(10*2)*(Src - 1) + Reg*2,])
  }
}

rm(rawdata)

# Now the plotting
years = seq(1950,2017)
par(mar = c(1,1,0,0) + 0.1, oma = c(1,1,1,1)+0.1, bg = NA)
layout(mat = matrix(1:2, ncol = 2, byrow = 2), widths = c(1, 0.5))
cols <- c('black', '#e41a1c',"",'#377eb8','#4daf4a','#984ea3','#ff7f00','#e6ab02')
  #c("black","purple", "hotpink", "blue", "darkgreen", "gold","red", "#a6761d")
plot(x = years, y =  totRecs[1,1,], type = "n", xlab = "Year", ylab = "Number of stations used in interpolation", cex.axis = 1.5,
     xlim = c(1950,2018),ylim = c(0,max(totRecs[1,1,])))
for (Reg in c(1:2,4:8)) {
  lines(totRecs[Reg,1,] ~ years,col = cols[Reg], lwd = 4)
}
par(mar = c(0,0,0,0)+0.1)
plot.new()
legend("topleft", legend = c("All", "Africa", "Asia", "Australia", "Europe", "North America", "South America"), 
       lwd = 4,col = cols[c(1,2,4:8)], cex = 2, bty = "n")
Figure 1: Number of stations used for interpolation by CCRC by regions

Figure 1: Number of stations used for interpolation by CCRC by regions

The total number of stations each day, interpolated by CCRC are high with the minimum number of stations per year being above 35,000 and average number of stations per year above 50,000 (see Figure 1). Regionally, the number of stations per year doubles in North America after 2000 and decreases substantially in South America from the late 1990s. There are no Chinese stations in 1950 and there is a large drop in stations in India in 1970. Stations in africa are still sparse throughout, however there are still more stations compared to other daily rainfall products.

rawdata <- read.table(file = "/srv/ccrc/data11/z3289452/CCRCGlobalGriddedDailyPrecip/Metadata/VERSION1.1_data/station-yearly-summary.txt")

totRecs <- array(dim = c(10,7,68))
totFlgs <- array(dim = c(10,7,68))
for (Reg in 1:10) {
  for (Src in 1:7) {
    totRecs[Reg,Src,] <- as.numeric(rawdata[(10*2)*(Src - 1) + Reg*2 - 1,])
    totFlgs[Reg,Src,] <- as.numeric(rawdata[(10*2)*(Src - 1) + Reg*2,])
  }
}

rm(rawdata)

# Plotting
# by sources
years = seq(1950,2017)
totRecsGPCC <- totRecs[1,2,] + totRecs[1,3,] + totRecs[1,4,]
totRecsOther <- totRecs[1,5,] + totRecs[1,6,]
totRecsGHCN <- totRecs[1,7,]

par(mar=  c(5,5,0.5,1) + 0.1, bg = NA)
plot(totRecsGPCC ~ years, type = "n", xlim = c(1950,2018), ylim = c(0,max(totRecsGPCC)),
     xlab = "Year", ylab = "Number of stations used in interpolation", cex.lab = 2, cex.axis = 1.5)
lines(totRecsGPCC ~ years, lwd = 4, col = "black")
lines(totRecsGHCN ~ years, lwd = 4, col = "blue")
lines(totRecsOther ~ years, lwd = 4, col = "darkgreen")
legend("topleft",legend = c("GPCC", "GHCN", "Other"), lwd = 4, col = c("black", "blue", "darkgreen"),
       bty = "n", cex = 2)

Figure 2: Number of stations used for interpolation by source The number of stations per year from the different sources are seen in Figure 2. The majority of station data for CCRC is from the stations hosted by GPCC. Due to the large overlap between the archives, the number of stations from GHCN is higher when fewer stations from GPCC are available. There is a steep decline in stations from GPCC around 2010 and a gradual increase until 1990. The raw data hosted by GPCC could be traced to either the SQL server (“Mirakel”) or to a GPCC file server. The data on the SQL server was of higher quality due to the manual quality control, however, due the time consuming nature of such checks there was little data dating earlier than 1988. The majority of pre-1988 data was available in ASCII format sorted by the country of origin.

# Make a map of in situ stations colour coded by source including a "merged" source label
library(dplyr)
library(maps)

# First import the station list
# create a vector of widths for read.fwf based on the ending position of each column 
end_positions <- c(13, 22, 32, 39, 42, 47, 52, 59, 66, 73, 80, 87, 94, 101, 108, 115, 122, 129, 136, 143, 150)#48, 53, 58, 65, 72, 79, 86, 93, 
                   #100, 107, 114, 121, 128, 135, 142, 149, 156)
widths <- c(end_positions[1], sapply(2:length(end_positions), FUN = function(n) {
  return(end_positions[n] - end_positions[n - 1])
}))

colNames <- c("ID", "Latitude", "Longitude", "Elevation", "CTRY_code",
              "Start_year", "End_year", "numRecs", "numFlgs",
              paste0(rep(c("numRecs_", "numFlgs_"), 6), rep(c("c", "d", "e", "g", "i", "G"), each = 2)))

stationList <- read.fwf("/srv/ccrc/data11/z3289452/CCRCGlobalGriddedDailyPrecip/Metadata/VERSION1.1_data/REGEN-stations_v1-1_extra_short_nonameandwmoid.txt",
                        widths = widths,
                        header = F,
                        col.names = colNames,
                        na.strings = c("-9999", "-999.9"),
                        fileEncoding = "utf8")

# if data from sources c,d or e then from GPCC -> 0
# if data from sources g or i then from other -> 1
# if data from source G then GHCN -> 2
# if a combination of the above then merged -> 3

stationList <- stationList %>% 
  mutate(netGPCCRecs = (numRecs_c - numFlgs_c) + (numRecs_d - numFlgs_d) + (numRecs_e - numFlgs_e),
         netGHCNRecs = (numRecs_G - numFlgs_G),
         netOtherRecs = (numRecs_g - numFlgs_g) + (numRecs_i - numFlgs_i))

stationList <- stationList %>%
  mutate(Source = case_when(
    netGPCCRecs > 0 & netGHCNRecs == 0 & netOtherRecs == 0 ~ "GPCC",
    netGPCCRecs == 0 & netGHCNRecs > 0 & netOtherRecs == 0 ~ "GHCN",
    netGPCCRecs == 0 & netGPCCRecs == 0 & netOtherRecs > 0 ~ "Other",
    netGPCCRecs > 0 & (netGHCNRecs > 0 | netOtherRecs > 0) ~ "mergedGPCC",
    netGHCNRecs > 0 & (netGPCCRecs > 0 | netOtherRecs > 0) ~ "mergedGHCN",
    TRUE ~ "noRecs"
  ))

GPCClatlon <- stationList %>% filter(Source == "GPCC") %>% select(Longitude, Latitude)
GHCNlatlon <- stationList %>% filter(Source == "GHCN") %>% select(Longitude, Latitude)
Otherlatlon <- stationList %>% filter(Source == "Other") %>% select(Longitude, Latitude)
Mergedlatlon <- stationList %>% filter(Source == "mergedGPCC" | Source == "mergedGHCN") %>% 
  select(Longitude, Latitude)

cols <- c('#1b9e77','#d95f02','#7570b3', '#e7298a')
bgcols <- paste0(cols, "bb")

par(mar = c(3,0.5,0.5,0.5) + 0.1, bg = NA)
plot(0, type = "n", xaxt = "n", yaxt = "n", ylim = c(-60, 90), xlim = c(-180, 180), ylab = "", xlab = "")
map("world", interior =  F, ylim = c(-60,90), add = T)
points(GPCClatlon, pch = 19, cex = 0.2, col = bgcols[1]) #, bg = bgcols[1]
points(GHCNlatlon, pch = 19, cex = 0.2, col = bgcols[2])
points(Mergedlatlon, pch = 19, cex = 0.2, col = bgcols[3])
points(Otherlatlon, pch = 19, cex = 0.2, col = bgcols[4])
legend("bottom", horiz = T, bty = "n", col = bgcols, pch = 19, cex = 1.8,
       legend = c("GPCC", "GHCN", "Merged", "Other"), 
       xpd = T, inset = c(0, -0.1))

Quality control

The quality control scripts used for the CCRC dataset was adopted from National Centre for Environmental Information, part of NOAA in USA. The quality control is done by two scripts and two other auxilliary scripts are used to generate the climatologies needed to do the quality control and handle the various data formats from different sources. At the end of the quality control process all data is written in a common format identical to the GHCN-Daily format (see README file, ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/daily/readme.txt).

The first quality control scripts does basic integrity checks such as checks for erroneous zeros, multi-day accumulations or accumulations greater than monthly rainfall, duplication or repitition of data and world record exceedances. In addition the first script also checks for outliers by checking for gaps in tails of distributions and checks for climatological outliers. The script also performs some temporal consistency checks by comparing values with consecutive days to look for unrealistic spikes in precipitation. Finally the first script does some megaconsitency checks to check for consistency of cold extremes. The second quality control script does spatial corroboration checks which determines if the value at each station is consistent with the values at neighbouring stations. For further information and detail on the quality control algorithms, refer to Durre et al. 2010.

Merge GHCND and GPCC and other smaller datasets

The merging script was used on multiple ocassions to merge quality controlled data from various sources. First it was used to merge data from various countries on the GPCC file server into a single archive which was then merged with the GPCC “Mirakel” data to create a combined archive of quality controlled GPCC stations. This GPCC archive was then merged with GHCN-Daily archive and subsequently the Argentinian and Russian data respectively.

Two stations are the same if:

  1. The latitude and longitudes match to three decimal places and their altitudes (to the nearest integer) and WMO IDs either match or are missing. Alternatively the stations are also matched if the WMO IDs are non-missing and match and the coordinates match to one decimal place.
  2. If the coordinates are within 1 degree of each other and WMO IDs either match or are missing and the correlation between the time series that overlap is greater than 0.99 and the overlapping time series themselves have atleast 365 elements with a minimum of 10 days with precipitation greater than 1mm.

If stations are not matched but coordinates between two stations are within a degree of each other and correlation between time series is still greater than 0.9, then the stations are saved in a separate “suspicious” list and await further manual review.

The matched station data from various sources is merged according to the following priorities in descending order.

  1. GPCC “Mirakel” data
  2. GPCC ASCII data
  3. Argentina/Russia data
  4. GHCN-Daily data

Overlapping data from a lower priority source is overwritted with data from higher priority sources.

Interpolation Method

The CCRC dataset uses ordinary block Kriging which is identical interpolation method as GPCC-FDD.

Evaluation

Time series of monthly/annual means for CPC, GPCC and CCRC datasets

  • Global
  • Regional
# Comparision with monthly datasets

# We reacreate figures 2.28 and 2.29 and table 2.10 from IPPC AR5 Chapter 2 Observations

library(RNetCDF)
library(plot3D)
library(maps)

# read and format data
# CCRC
# Annual totals
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_1950-2013_1deg_yearsum.nc")
CCRC_annualtot <- var.get.nc(nc, "p")
lon <- var.get.nc(nc, "lon")
lat <- var.get.nc(nc, "lat")
close.nc(nc)
# Climatology
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_1951-2010_1deg_yearsum_timavg.nc")
CCRC_clim <- var.get.nc(nc, "p")
close.nc(nc)
# Anomaly
CCRC_anom <- CCRC_annualtot - array(dim = dim(CCRC_annualtot), data = CCRC_clim)
# Now the averages
# Global
CCRC_anom_globalavg <- apply(CCRC_anom, 3, mean, na.rm = T)
# 60N - 90N
CCRC_anom_60N90Navg <- apply(CCRC_anom[,which(lat > 60 & lat <= 90),], 3, mean, na.rm = T)
# 30N - 60N
CCRC_anom_30N60Navg <- apply(CCRC_anom[,which(lat > 30 & lat <= 60),], 3, mean, na.rm = T)
# 30S - 30N
CCRC_anom_30S30Navg <- apply(CCRC_anom[,which(lat > -30 & lat <= 30),], 3, mean, na.rm = T)
# 60S - 30S
CCRC_anom_60S30Savg <- apply(CCRC_anom[,which(lat > -60 & lat <= -30),], 3, mean, na.rm = T)

# CCRC 40 Yr Long Term
# Annual totals
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/FORTY/REGEN_V1.1_LongTermStns_1950-2013_1deg_yearsum.nc")
CCRC_40yr_annualtot <- var.get.nc(nc, "p")
lon <- var.get.nc(nc, "lon")
lat <- var.get.nc(nc, "lat")
close.nc(nc)
# Climatology
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/FORTY/REGEN_V1.1_LongTermStns_1951-2010_1deg_yearsum_timavg.nc")
CCRC_40yr_clim <- var.get.nc(nc, "p")
close.nc(nc)
# Anomaly
CCRC_40yr_anom <- CCRC_40yr_annualtot - array(dim = dim(CCRC_40yr_annualtot), data = CCRC_40yr_clim)
# Now the averages
# Global
CCRC_40yr_anom_globalavg <- apply(CCRC_40yr_anom, 3, mean, na.rm = T)
# 60N - 90N
CCRC_40yr_anom_60N90Navg <- apply(CCRC_40yr_anom[,which(lat > 60 & lat <= 90),], 3, mean, na.rm = T)
# 30N - 60N
CCRC_40yr_anom_30N60Navg <- apply(CCRC_40yr_anom[,which(lat > 30 & lat <= 60),], 3, mean, na.rm = T)
# 30S - 30N
CCRC_40yr_anom_30S30Navg <- apply(CCRC_40yr_anom[,which(lat > -30 & lat <= 30),], 3, mean, na.rm = T)
# 60S - 30S
CCRC_40yr_anom_60S30Savg <- apply(CCRC_40yr_anom[,which(lat > -60 & lat <= -30),], 3, mean, na.rm = T)

# CRU_TS4.01
# Annual totals
nc <- open.nc("/srv/ccrc/data45/z3289452/MonthlyPRCP/CRU_TS4.01/cru_ts4.01.1950.2013.pre.yearsum.dat.nc")
CRU_annualtot <- var.get.nc(nc, "pre")
lon <- var.get.nc(nc, "lon")
lat <- var.get.nc(nc, "lat")
close.nc(nc)
# Climatology
nc <- open.nc("/srv/ccrc/data45/z3289452/MonthlyPRCP/CRU_TS4.01/cru_ts4.01.1951.2010.climatology.pre.dat.nc")
CRU_clim <- var.get.nc(nc, "pre")
close.nc(nc)
# Anomaly
CRU_anom <- CRU_annualtot - array(dim = dim(CRU_annualtot), data = CRU_clim)
# Now the averages
# Global
CRU_anom_globalavg <- apply(CRU_anom, 3, mean, na.rm = T)
# 60N - 90N
CRU_anom_60N90Navg <- apply(CRU_anom[,which(lat > 60 & lat <= 90),], 3, mean, na.rm = T)
# 30N - 60N
CRU_anom_30N60Navg <- apply(CRU_anom[,which(lat > 30 & lat <= 60),], 3, mean, na.rm = T)
# 30S - 30N
CRU_anom_30S30Navg <- apply(CRU_anom[,which(lat > -30 & lat <= 30),], 3, mean, na.rm = T)
# 60S - 30S
CRU_anom_60S30Savg <- apply(CRU_anom[,which(lat > -60 & lat <= -30),], 3, mean, na.rm = T)

# GPCC Full Monthly V7
# Annual totals
nc <- open.nc("/srv/ccrc/data45/z3289452/MonthlyPRCP/GPCC_Full_Monthly_V7/precip.mon.total.1x1.v7_1950-2013.yearsum.nc")
GPCC_annualtot <- var.get.nc(nc, "precip")
lon <- var.get.nc(nc, "lon")
lat <- var.get.nc(nc, "lat")
close.nc(nc)
# Climatology
nc <- open.nc("/srv/ccrc/data45/z3289452/MonthlyPRCP/GPCC_Full_Monthly_V7/precip.mon.total.1x1.v7_1951-2010.climatology.nc")
GPCC_clim <- var.get.nc(nc, "precip")
close.nc(nc)
# Anomaly
GPCC_anom <- GPCC_annualtot - array(dim = dim(GPCC_annualtot), data = GPCC_clim)
# Now the averages
# Global
GPCC_anom_globalavg <- apply(GPCC_anom, 3, mean, na.rm = T)
# 60N - 90N
GPCC_anom_60N90Navg <- apply(GPCC_anom[,which(lat > 60 & lat <= 90),], 3, mean, na.rm = T)
# 30N - 60N
GPCC_anom_30N60Navg <- apply(GPCC_anom[,which(lat > 30 & lat <= 60),], 3, mean, na.rm = T)
# 30S - 30N
GPCC_anom_30S30Navg <- apply(GPCC_anom[,which(lat > -30 & lat <= 30),], 3, mean, na.rm = T)
# 60S - 30S
GPCC_anom_60S30Savg <- apply(GPCC_anom[,which(lat > -60 & lat <= -30),], 3, mean, na.rm = T)

# GHCN V2
# Ollie's data
# Annual totals
nc <- open.nc("/srv/ccrc/data45/z3289452/MonthlyPRCP/GHCN_V2/precip.mon.total.1950-2013.yearsum.nc")
GHCN_annualtot <- var.get.nc(nc, "precip")
lon <- var.get.nc(nc, "lon")
lat <- var.get.nc(nc, "lat")
close.nc(nc)
# Climatology
nc <- open.nc("/srv/ccrc/data45/z3289452/MonthlyPRCP/GHCN_V2/precip.mon.total.climatology.nc")
GHCN_clim <- var.get.nc(nc, "precip")
close.nc(nc)
# Anomaly
GHCN_anom <- GHCN_annualtot - array(dim = dim(GHCN_annualtot), data = GHCN_clim)
# Now the averages
# Global
GHCN_anom_globalavg <- apply(GHCN_anom, 3, mean, na.rm = T)
# 60N - 90N
GHCN_anom_60N90Navg <- apply(GHCN_anom[,which(lat > 60 & lat <= 90),], 3, mean, na.rm = T)
# 30N - 60N
GHCN_anom_30N60Navg <- apply(GHCN_anom[,which(lat > 30 & lat <= 60),], 3, mean, na.rm = T)
# 30S - 30N
GHCN_anom_30S30Navg <- apply(GHCN_anom[,which(lat > -30 & lat <= 30),], 3, mean, na.rm = T)
# 60S - 30S
GHCN_anom_60S30Savg <- apply(GHCN_anom[,which(lat > -60 & lat <= -30),], 3, mean, na.rm = T)

# Now the timeseries plot
par(mar = c(0,0,0,0), oma = c(4,6,2,2) + 0.1)
layout(mat = matrix(1:5, ncol = 1, byrow = T))
years <- 1950:2013
cols <- c('#e41a1cbb','#377eb8bb','#4daf4abb','#984ea3bb','#ff7f00bb')
        #c('#e41a1cbb','#377eb8bb','#4daf4abb','#984ea3bb')
# Global
quantile(c(GHCN_anom_globalavg, GPCC_anom_globalavg, CRU_anom_globalavg, CCRC_40yr_anom_globalavg), 
         probs = seq(0, 1, by = 0.05), na.rm = T)
yrange <- c(-75, 50)
plot(years, CCRC_anom_globalavg, ylim = yrange, type = "n",
     xlab = "", ylab = "",
     #xlab = "Years", ylab = "Global Precipitation anomalies (mm/Y)", 
     cex.axis = 2, cex.lab = 2, las = 1,
     xaxt = "n")
lines(years, CCRC_anom_globalavg, lwd = 4, col = cols[1])
lines(years, CCRC_40yr_anom_globalavg, lwd = 4, col = cols[2])
lines(years, GPCC_anom_globalavg, lwd = 4, col = cols[3])
lines(years, CRU_anom_globalavg, lwd = 4, col = cols[4])
lines(years, GHCN_anom_globalavg, lwd = 4, col = cols[5])
text(2010, yrange[2] - yrange[2]/5, "Global", cex = 2)
abline(h = 0, lty = "dashed", lwd = 3)
legend("bottomleft", legend = c("REGEN", "REGEN 40YR", "GPCC", "CRU", "GHCN"), col = cols, lwd = 4, cex = 1.8, bty = "n")
# 60N - 90N
quantile(c(GHCN_anom_60N90Navg, GPCC_anom_60N90Navg, CRU_anom_60N90Navg, CCRC_anom_60N90Navg,
           CCRC_40yr_anom_60N90Navg), probs = seq(0, 1, by = 0.05), na.rm = T)
yrange <- c(-60, 85)
plot(years, CCRC_anom_60N90Navg, ylim = yrange, type = "n",
     xlab = "", ylab = "",
     #xlab = "Years", ylab = "60N - 90N Precipitation anomalies (mm/Y)", 
     cex.axis = 2, cex.lab = 2, las = 1,
     xaxt = "n")
lines(years, CCRC_anom_60N90Navg, lwd = 4, col = cols[1])
lines(years, CCRC_40yr_anom_60N90Navg, lwd = 4, col = cols[2])
lines(years, GPCC_anom_60N90Navg, lwd = 4, col = cols[3])
lines(years, CRU_anom_60N90Navg, lwd = 4, col = cols[4])
lines(years, GHCN_anom_60N90Navg, lwd = 4, col = cols[5])
abline(h = 0, lty = "dashed", lwd = 3)
text(2010, yrange[2] - yrange[2]/5, "60N-90N", cex = 2)
# 30N - 60N
quantile(c(GHCN_anom_30N60Navg, GPCC_anom_30N60Navg, CRU_anom_30N60Navg, CCRC_anom_30N60Navg, 
           CCRC_40yr_anom_30N60Navg), probs = seq(0, 1, by = 0.05), na.rm = T)
yrange <- c(-30, 50)
plot(years, CCRC_anom_30N60Navg, ylim = yrange, type = "n",
     xlab = "", ylab = "", 
     #xlab = "Years", ylab = "30N - 60N Precipitation anomalies (mm/Y)", 
     cex.axis = 2, cex.lab = 2, las = 1,
     xaxt = "n")
lines(years, CCRC_anom_30N60Navg, lwd = 4, col = cols[1])
lines(years, CCRC_40yr_anom_30N60Navg, lwd = 4, col = cols[2])
lines(years, GPCC_anom_30N60Navg, lwd = 4, col = cols[3])
lines(years, CRU_anom_30N60Navg, lwd = 4, col = cols[4])
lines(years, GHCN_anom_30N60Navg, lwd = 4, col = cols[5])
abline(h = 0, lty = "dashed", lwd = 3)
text(2010, yrange[2] - yrange[2]/5, "30N-60N", cex = 2)
# 30S - 30N
quantile(c(GHCN_anom_30S30Navg, GPCC_anom_30S30Navg, CRU_anom_30S30Navg, CCRC_anom_30S30Navg, 
           CCRC_40yr_anom_30S30Navg), probs = seq(0, 1, by = 0.05), na.rm = T)
yrange <- c(-130, 130)
plot(years, CCRC_anom_30S30Navg, ylim = yrange, type = "n",
     xlab = "", ylab = "",
     #xlab = "Years", ylab = "30S - 30N Precipitation anomalies (mm/Y)", 
     cex.axis = 2, cex.lab = 2, las = 1,
     xaxt = "n")
lines(years, CCRC_anom_30S30Navg, lwd = 4, col = cols[1])
lines(years, CCRC_40yr_anom_30S30Navg, lwd = 4, col = cols[2])
lines(years, GPCC_anom_30S30Navg, lwd = 4, col = cols[3])
lines(years, CRU_anom_30S30Navg, lwd = 4, col = cols[4])
lines(years, GHCN_anom_30S30Navg, lwd = 4, col = cols[5])
abline(h = 0, lty = "dashed", lwd = 3)
text(2010, yrange[2] - yrange[2]/5, "30S-30N", cex = 2)
# 60S - 30S
quantile(c(GHCN_anom_60S30Savg, GPCC_anom_60S30Savg, CRU_anom_60S30Savg, CCRC_anom_60S30Savg, 
           CCRC_40yr_anom_60S30Savg), probs = seq(0, 1, by = 0.05), na.rm = T)
yrange <- c(-95, 100)
plot(years, CCRC_anom_60S30Savg, ylim = yrange, type = "n",
     xlab = "", ylab = "",
     #xlab = "Years", ylab = "60S - 30S Precipitation anomalies (mm/Y)", 
     cex.axis = 2, cex.lab = 2, las = 1)
lines(years, CCRC_anom_60S30Savg, lwd = 4, col = cols[1])
lines(years, CCRC_40yr_anom_60S30Savg, lwd = 4, col = cols[2])
lines(years, GPCC_anom_60S30Savg, lwd = 4, col = cols[3])
lines(years, CRU_anom_60S30Savg, lwd = 4, col = cols[4])
lines(years, GHCN_anom_60S30Savg, lwd = 4, col = cols[5])
abline(h = 0, lty = "dashed", lwd = 3)
text(2010, yrange[2] - yrange[2]/5, "60S-30S", cex = 2)
# labels
mtext("Years", side = 1, line = 3, cex = 2)
mtext("Precipitation anomaly (mm/Y)", side = 2, line = 4, cex = 2, adj = -0.5)

add labels

Regional Time series
will probably be removed. No need to show time series from various regions

Spatial patterns between regional datasets and CPC, GPCC and CCRC for an example day

  • USA
  • EOBS
  • Aphrodite
  • AWAP

Total Annual Precipitation

library(RNetCDF)
library(lubridate)
library(plot3D) # for colkey()
library(fields) # for map()

# First the annual total average
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_1950-2013_1deg_yearsum_timavg.nc")
yearsum.avg <- var.get.nc(nc, "p")
lon <- var.get.nc(nc, "lon")
lat <- var.get.nc(nc, "lat")
close.nc(nc)

quantile(c(yearsum.avg), probs = seq(0, 1, by = 0.05), na.rm = T)
brks <- seq(0, 2000, length.out = 6)
zbrks <- brks
zbrks[1] <- min(brks, yearsum.avg, na.rm = T)
zbrks[length(brks)] <- max(brks, yearsum.avg, na.rm = T)
lbrks <- ifelse(brks == zbrks, brks, "")
cols <- c('#ffffcc','#a1dab4','#41b6c4','#2c7fb8','#253494')

# CCRC plot
par(mar = c(0,0,0,0) + 0.1, bg = NA)
layout(mat = matrix(c(1,2), nrow = 2, byrow = T), heights = c(1,0.2))
image(lon, lat, yearsum.avg, xlab = "", ylab = "", breaks = zbrks, 
      col = cols, xaxt = "n", yaxt = "n", cex.lab = 1.5, cex.axis = 1.5, ylim = c(-60,90))
map("world", interior = F, add = T)
par(mar = c(0,8,0,8) + 0.1)
colkey(col = cols, side = 1, clim = range(brks), width = 5, at = brks, 
       labels = lbrks, cex.axis = 1.5)
title("REGEN V1.1 Average total annual precipitation 1950-2013 (mm/year)", line = -1.5, cex.main = 1.5)
# Now the same as above but for long term 40yr data
library(RNetCDF)
library(lubridate)
library(plot3D) # for colkey()
library(fields) # for map()

# First the annual total average
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/FORTY/REGEN_V1.1_LongTermStns_1950-2013_1deg_yearsum_timavg.nc")
yearsum.avg <- var.get.nc(nc, "p")
lon <- var.get.nc(nc, "lon")
lat <- var.get.nc(nc, "lat")
close.nc(nc)

quantile(c(yearsum.avg), probs = seq(0, 1, by = 0.05), na.rm = T)
brks <- seq(0, 2000, length.out = 6)
zbrks <- brks
zbrks[1] <- min(brks, yearsum.avg, na.rm = T)
zbrks[length(brks)] <- max(brks, yearsum.avg, na.rm = T)
lbrks <- ifelse(brks == zbrks, brks, "")
cols <- c('#ffffcc','#a1dab4','#41b6c4','#2c7fb8','#253494')

# CCRC plot
par(mar = c(0,0,0,0) + 0.1, bg = NA)
layout(mat = matrix(c(1,2), nrow = 2, byrow = T), heights = c(1,0.2))
image(lon, lat, yearsum.avg, xlab = "", ylab = "", breaks = zbrks, 
      col = cols, xaxt = "n", yaxt = "n", cex.lab = 1.5, cex.axis = 1.5, ylim = c(-60,90))
map("world", interior = F, add = T)
par(mar = c(0,8,0,8) + 0.1)
colkey(col = cols, side = 1, clim = range(brks), width = 5, at = brks, 
       labels = lbrks, cex.axis = 1.5)
title("REGEN V1.1 Long Term 40yr Average total annual precipitation 1950-2013 (mm/year)", line = -1.5, cex.main = 1.5)

Also show total precipitation from GPCC and CPC and difference plots

library(RNetCDF)
library(lubridate)
library(plot3D) # for colkey()
library(fields) # for map()

# First the annual total average
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_1988-2013_1deg_yearsum_timavg.nc")
yearsum.avg <- var.get.nc(nc, "p")
lon <- var.get.nc(nc, "lon")
lat <- var.get.nc(nc, "lat")
close.nc(nc)

# Now CPC_global 0.5 deg (regridded to 1 degree) 1979-2013
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/CPC_0.5X0.5_Daily_Unified_Global_Precipitation/CPC_global_1deg_1988-2013_yearsum_timavg.nc")
CPC_yearsum.avg <- var.get.nc(nc, "precip")
CPC.lon <- var.get.nc(nc, "lon")
CPC.lat <- var.get.nc(nc, "lat")
close.nc(nc)

# Now GPCC_FDD 1 deg 1988-2013
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/GPCC_FDD/full_data_daily_1988-2013_yearsum_timavg.nc")
GPCC_yearsum.avg <- var.get.nc(nc, "p")
GPCC.lon <- var.get.nc(nc, "lon")
GPCC.lat <- var.get.nc(nc, "lat")
close.nc(nc)

# Now REGEN (CCRC) long term 40yr data
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/FORTY/REGEN_V1.1_LongTermStns_1988-2013_1deg_yearsum_timavg.nc")
REGEN40yr_yearsum.avg <- var.get.nc(nc, "p")
close.nc(nc)

# Difference fields
# CCRC - GPCC
CCRC_GPCC.annualavg <- (yearsum.avg - GPCC_yearsum.avg)/yearsum.avg*100
# CCRC - CPC
CCRC_CPC.annualavg <- (yearsum.avg - CPC_yearsum.avg)/yearsum.avg*100
# REGEN - REGEN long term
CCRC_REGEN40yr.annualavg <- (yearsum.avg - REGEN40yr_yearsum.avg)/yearsum.avg*100
quantile(c(yearsum.avg), probs = seq(0, 1, by = 0.05), na.rm = T)
brks <- seq(0, 2000, length.out = 6)
zbrks <- brks
zbrks[1] <- min(brks, yearsum.avg, GPCC_yearsum.avg, CPC_yearsum.avg, CCRC_REGEN40yr.annualavg, na.rm = T)
zbrks[6] <- max(brks, yearsum.avg, GPCC_yearsum.avg, CPC_yearsum.avg, CCRC_REGEN40yr.annualavg, na.rm = T)
lbrks <- ifelse(brks == zbrks, brks, "")
cols <- c('#ffffcc','#a1dab4','#41b6c4','#2c7fb8','#253494') 

# CCRC plot
par(mar = c(0,0,0,0) + 0.1, bg= NA)
layout(mat = matrix(c(1,2), nrow = 2, byrow = T), heights = c(1,0.2))
image(lon, lat, yearsum.avg, xlab = "", ylab = "", breaks = zbrks, 
      col = cols, xaxt = "n", yaxt = "n", cex.lab = 1.5, cex.axis = 1.5, ylim = c(-60, 90))
map("world", interior = F, add = T)
# title(xlab = "CCRC Averaged total annual precipitation 1988-2013 (mm)", line = 1, cex.lab = 1.5)
par(mar = c(0,8,0,8) + 0.1)
colkey(col = cols, side = 1, clim = range(brks), width = 5, at = brks, 
       labels = lbrks, cex.axis = 1.5)
title("REGEN V1.1 Averaged total annual precipitation 1988-2013 (mm/year)", line = -1.5, cex.main = 1.5)
# Now the difference plots
quantile(c(CCRC_GPCC.annualavg, CCRC_CPC.annualavg), probs = seq(0, 1, by = 0.05), na.rm = T)
brks <- seq(-30, 30, length.out = 7)

# CCRC - GPCC
zbrks <- brks
zbrks[1] <- min(brks, CCRC_GPCC.annualavg, na.rm = T)
zbrks[length(brks)] <- max(brks, CCRC_GPCC.annualavg, na.rm = T)
lbrks <- ifelse(brks == zbrks, paste0(brks, "%"), "")
cols <- c('#b2182b','#ef8a62','#fddbc7','#d1e5f0','#67a9cf','#2166ac')
par(mar = c(0,0,0,0) + 0.1, bg = NA)
layout(mat = matrix(c(1,2), nrow = 2, byrow = T), heights = c(1,0.2))
image(lon, lat, CCRC_GPCC.annualavg, xlab = "", ylab = "", breaks = zbrks, 
      col = cols, xaxt = "n", yaxt = "n", cex.lab = 1.5, cex.axis = 1.5, ylim = c(-60, 90))
map("world", interior = F, add = T)
par(mar = c(0,8,0,8) + 0.1)
colkey(col = cols, side = 1, clim = range(brks), width = 5, at = brks, 
       labels = lbrks, cex.axis = 1.5)
title("Percentage diff between REGEN and GPCC-FDD Avg total prcp 1988-2013", line = -1.5, cex.main = 1.5)

# CCRC - CPC
zbrks <- brks
zbrks[1] <- min(brks, CCRC_CPC.annualavg, na.rm = T)
zbrks[length(brks)] <- max(brks, CCRC_CPC.annualavg, na.rm = T)
lbrks <- ifelse(brks == zbrks, paste0(brks, "%"), "")
cols <- c('#b2182b','#ef8a62','#fddbc7','#d1e5f0','#67a9cf','#2166ac')
par(mar = c(0,0,0,0) + 0.1, bg = NA)
layout(mat = matrix(c(1,2), nrow = 2, byrow = T), heights = c(1,0.2))
image(lon, lat, CCRC_CPC.annualavg, xlab = "", ylab = "", breaks = zbrks, 
      col = cols, xaxt = "n", yaxt = "n", cex.lab = 1.5, cex.axis = 1.5, ylim = c(-60, 90))
map("world", interior = F, add = T)
par(mar = c(0,8,0,8) + 0.1)
colkey(col = cols, side = 1, clim = range(brks), width = 5, at = brks, 
       labels = lbrks, cex.axis = 1.5)
title("Percentage diff between REGEN and CPC_Glb Avg total prcp 1988-2013", line = -1.5, cex.main = 1.5)

# REGEN - REGEN40yr
zbrks <- brks
zbrks[1] <- min(brks, CCRC_REGEN40yr.annualavg, na.rm = T)
zbrks[length(brks)] <- max(brks, CCRC_REGEN40yr.annualavg, na.rm = T)
lbrks <- ifelse(brks == zbrks, paste0(brks, "%"), "")
cols <- c('#b2182b','#ef8a62','#fddbc7','#d1e5f0','#67a9cf','#2166ac')
par(mar = c(0,0,0,0) + 0.1, bg = NA)
layout(mat = matrix(c(1,2), nrow = 2, byrow = T), heights = c(1,0.2))
image(lon, lat, CCRC_REGEN40yr.annualavg, xlab = "", ylab = "", breaks = zbrks, 
      col = cols, xaxt = "n", yaxt = "n", cex.lab = 1.5, cex.axis = 1.5, ylim = c(-60, 90))
map("world", interior = F, add = T)
par(mar = c(0,8,0,8) + 0.1)
colkey(col = cols, side = 1, clim = range(brks), width = 5, at = brks, 
       labels = lbrks, cex.axis = 1.5)
title("Percentage diff between REGEN and REGEN Long Term 40yr Avg total prcp 1988-2013", line = -01.5, cex.main = 1.5)

Maximum Annual Precipitation

library(RNetCDF)
library(lubridate)
library(plot3D) # for colkey()
library(fields) # for map()

# First the annual total average
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_1950-2013_1deg_p_yearmax_timavg.nc")
yearmax.avg <- var.get.nc(nc, "p")
lon <- var.get.nc(nc, "lon")
lat <- var.get.nc(nc, "lat")
close.nc(nc)

quantile(c(yearmax.avg), probs = seq(0, 1, by = 0.05), na.rm = T)
brks <- seq(0, 75, length.out = 6)
zbrks <- brks
zbrks[1] <- min(brks, yearmax.avg, na.rm = T)
zbrks[length(brks)] <- max(brks, yearmax.avg, na.rm = T)
lbrks <- ifelse(brks == zbrks, brks, "")
cols <- c('#feebe2','#fbb4b9','#f768a1','#c51b8a','#7a0177')

# CCRC plot
par(mar = c(0,0,0,0) + 0.1, bg = NA)
layout(mat = matrix(c(1,2), nrow = 2, byrow = T), heights = c(1,0.2))
image(lon, lat, yearmax.avg, xlab = "", ylab = "", breaks = zbrks, 
      col = cols, xaxt = "n", yaxt = "n", cex.lab = 1.5, cex.axis = 1.5, ylim = c(-60,90))
map("world", interior = F, add = T)
par(mar = c(0,8,0,8) + 0.1)
colkey(col = cols, side = 1, clim = range(brks), width = 5, at = brks, 
       labels = lbrks, cex.axis = 1.5)
title("REGEN V1.1 Average annual maximum precipitation 1950-2013 (mm/year)", line = -1.5, cex.main = 1.5)
# Now do the same for Long term 40yr data
library(RNetCDF)
library(lubridate)
library(plot3D) # for colkey()
library(fields) # for map()

# First the annual total average
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/FORTY/REGEN_V1.1_LongTermStns_1950-2013_1deg_p_yearmax_timavg.nc")
yearmax.avg <- var.get.nc(nc, "p")
lon <- var.get.nc(nc, "lon")
lat <- var.get.nc(nc, "lat")
close.nc(nc)

quantile(c(yearmax.avg), probs = seq(0, 1, by = 0.05), na.rm = T)
brks <- seq(0, 75, length.out = 6)
zbrks <- brks
zbrks[1] <- min(brks, yearmax.avg, na.rm = T)
zbrks[length(brks)] <- max(brks, yearmax.avg, na.rm = T)
lbrks <- ifelse(brks == zbrks, brks, "")
cols <- c('#feebe2','#fbb4b9','#f768a1','#c51b8a','#7a0177')

# CCRC plot
par(mar = c(0,0,0,0) + 0.1, bg = NA)
layout(mat = matrix(c(1,2), nrow = 2, byrow = T), heights = c(1,0.2))
image(lon, lat, yearmax.avg, xlab = "", ylab = "", breaks = zbrks, 
      col = cols, xaxt = "n", yaxt = "n", cex.lab = 1.5, cex.axis = 1.5, ylim = c(-60,90))
map("world", interior = F, add = T)
par(mar = c(0,8,0,8) + 0.1)
colkey(col = cols, side = 1, clim = range(brks), width = 5, at = brks, 
       labels = lbrks, cex.axis = 1.5)
title("REGEN V1.1 Long Term 40yr Average annual maximum precipitation 1950-2013 (mm/year)", line = -1.5, cex.main = 1.5)

Also show RX1DAY from GPCC and CPC and difference plots

library(RNetCDF)
library(lubridate)
library(plot3D) # for colkey()
library(fields) # for map()

# First the annual total average
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_1988-2013_1deg_p_yearmax_timavg.nc")
yearmax.avg <- var.get.nc(nc, "p")
lon <- var.get.nc(nc, "lon")
lat <- var.get.nc(nc, "lat")
close.nc(nc)

# Now CPC_global 0.5 deg (regridded to 1 degree) 1979-2013
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/CPC_0.5X0.5_Daily_Unified_Global_Precipitation/CPC_global_1deg_1988-2013_yearmax_timavg.nc")
CPC_yearmax.avg <- var.get.nc(nc, "precip")
CPC.lon <- var.get.nc(nc, "lon")
CPC.lat <- var.get.nc(nc, "lat")
close.nc(nc)

# Now GPCC_FDD 1 deg 1988-2013
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/GPCC_FDD/full_data_daily_1988-2013_yearmax_timavg.nc")
GPCC_yearmax.avg <- var.get.nc(nc, "p")
GPCC.lon <- var.get.nc(nc, "lon")
GPCC.lat <- var.get.nc(nc, "lat")
close.nc(nc)

# REGEN40yr
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/FORTY/REGEN_V1.1_LongTermStns_1988-2013_1deg_p_yearmax_timavg.nc")
REGEN40yr_yearmax.avg <- var.get.nc(nc, "p")
close.nc(nc)

# Difference fields
# CCRC - GPCC
CCRC_GPCC.annualmax <- (yearmax.avg - GPCC_yearmax.avg)/yearmax.avg*100
# CCRC - CPC
CCRC_CPC.annualmax <- (yearmax.avg - CPC_yearmax.avg)/yearmax.avg*100
# REGEN - REGEN40yr
CCRC_REGEN40yr.annualmax <- (yearmax.avg - REGEN40yr_yearmax.avg)/yearmax.avg*100
quantile(c(yearmax.avg), probs = seq(0, 1, by = 0.05), na.rm = T)
brks <- seq(0, 75, length.out = 6)
zbrks <- brks
zbrks[1] <- min(brks, yearmax.avg, GPCC_yearmax.avg, CPC_yearmax.avg, REGEN40yr_yearmax.avg, na.rm = T)
zbrks[6] <- max(brks, yearmax.avg, GPCC_yearmax.avg, CPC_yearmax.avg, REGEN40yr_yearmax.avg, na.rm = T)
lbrks <- ifelse(brks == zbrks, brks, "")
cols <- c('#feebe2','#fbb4b9','#f768a1','#c51b8a','#7a0177')

# CCRC plot
par(mar = c(0,0,0,0) + 0.1, bg= NA)
layout(mat = matrix(c(1,2), nrow = 2, byrow = T), heights = c(1,0.2))
image(lon, lat, yearmax.avg, xlab = "", ylab = "", breaks = zbrks, 
      col = cols, xaxt = "n", yaxt = "n", cex.lab = 1.5, cex.axis = 1.5, ylim = c(-60, 90))
map("world", interior = F, add = T)
# title(xlab = "CCRC Averaged total annual precipitation 1988-2013 (mm)", line = 1, cex.lab = 1.5)
par(mar = c(0,8,0,8) + 0.1)
colkey(col = cols, side = 1, clim = range(brks), width = 5, at = brks, 
       labels = lbrks, cex.axis = 1.5)
title("REGEN V1.1 Averaged annual maximum precipitation 1988-2013 (mm/year)", line = -1.5, cex.main = 1.5)
# Now the difference plots
quantile(c(CCRC_GPCC.annualmax, CCRC_CPC.annualmax, CCRC_REGEN40yr.annualmax), probs = seq(0, 1, by = 0.05), na.rm = T)
brks <- seq(-45, 45, length.out = 7)

# CCRC - GPCC
zbrks <- brks
zbrks[1] <- min(brks, CCRC_GPCC.annualmax, na.rm = T)
zbrks[length(brks)] <- max(brks, CCRC_GPCC.annualmax, na.rm = T)
lbrks <- ifelse(brks == zbrks, paste0(brks, "%"), "")
cols <- c('#b2182b','#ef8a62','#fddbc7','#d1e5f0','#67a9cf','#2166ac')
par(mar = c(0,0,0,0) + 0.1, bg = NA)
layout(mat = matrix(c(1,2), nrow = 2, byrow = T), heights = c(1,0.2))
image(lon, lat, CCRC_GPCC.annualmax, xlab = "", ylab = "", breaks = zbrks, 
      col = cols, xaxt = "n", yaxt = "n", cex.lab = 1.5, cex.axis = 1.5, ylim = c(-60, 90))
map("world", interior = F, add = T)
par(mar = c(0,8,0,8) + 0.1)
colkey(col = cols, side = 1, clim = range(brks), width = 5, at = brks, 
       labels = lbrks, cex.axis = 1.5)
title("Percentage diff between REGEN and GPCC-FDD Avg max prcp 1988-2013", line = -1.5, cex.main = 1.5)

# CCRC - CPC
zbrks <- brks
zbrks[1] <- min(brks, CCRC_CPC.annualmax, na.rm = T)
zbrks[length(brks)] <- max(brks, CCRC_CPC.annualmax, na.rm = T)
lbrks <- ifelse(brks == zbrks, paste0(brks, "%"), "")
cols <- c('#b2182b','#ef8a62','#fddbc7','#d1e5f0','#67a9cf','#2166ac')
par(mar = c(0,0,0,0) + 0.1, bg = NA)
layout(mat = matrix(c(1,2), nrow = 2, byrow = T), heights = c(1,0.2))
image(lon, lat, CCRC_CPC.annualmax, xlab = "", ylab = "", breaks = zbrks, 
      col = cols, xaxt = "n", yaxt = "n", cex.lab = 1.5, cex.axis = 1.5, ylim = c(-60, 90))
map("world", interior = F, add = T)
par(mar = c(0,8,0,8) + 0.1)
colkey(col = cols, side = 1, clim = range(brks), width = 5, at = brks, 
       labels = lbrks, cex.axis = 1.5)
title("Percentage diff between REGEN and CPC_Glb Avg max prcp 1988-2013", line = -01.5, cex.main = 1.5)

# REGEN - REGEN40yr
zbrks <- brks
zbrks[1] <- min(brks, CCRC_REGEN40yr.annualmax, na.rm = T)
zbrks[length(brks)] <- max(brks, CCRC_REGEN40yr.annualmax, na.rm = T)
lbrks <- ifelse(brks == zbrks, paste0(brks, "%"), "")
cols <- c('#b2182b','#ef8a62','#fddbc7','#d1e5f0','#67a9cf','#2166ac')
par(mar = c(0,0,0,0) + 0.1, bg = NA)
layout(mat = matrix(c(1,2), nrow = 2, byrow = T), heights = c(1,0.2))
image(lon, lat, CCRC_REGEN40yr.annualmax, xlab = "", ylab = "", breaks = zbrks, 
      col = cols, xaxt = "n", yaxt = "n", cex.lab = 1.5, cex.axis = 1.5, ylim = c(-60, 90))
map("world", interior = F, add = T)
par(mar = c(0,8,0,8) + 0.1)
colkey(col = cols, side = 1, clim = range(brks), width = 5, at = brks, 
       labels = lbrks, cex.axis = 1.5)
title("Percentage diff between REGEN and REGEN Long Term 40yr Avg max prcp 1988-2013", line = -1.5, cex.main = 1.5)

The following will be replaced by timeseries averaged over grid boxes of specific events.
Europe

# Time series of significant rainfall events in regions with regional datasets 
# averaged over an area over a few days

library(RNetCDF)
library(lubridate)
library(fields)
library(plot3D)

# Europe
# REGEN
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_19680908-19680922_1deg.nc")
# There was heavy rainfall on 15th Sept 1968 that caused floods in the Great Floods of 1968 in UK and France
dates <- seq(ymd("19680908"), ymd("19680922"), by = 1)
p <- var.get.nc(nc, "p", start = c(173, 138, NA), count = c(15, 8, NA))
s <- var.get.nc(nc, "s", start = c(173, 138, NA), count = c(15, 8, NA))
lon <- var.get.nc(nc, "lon")[173:187] # -7.5:6.5
lat <- var.get.nc(nc, "lat")[138:145] # 47.5:54.5
close.nc(nc)
# E-OBS
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/E-OBS_1950-2017_v16.0/rr_0.25deg_reg_1965-1979_v16.0.nc")
Eobs.lon <- var.get.nc(nc, "longitude")[132:188] # -7.625:6.375
Eobs.lat <- var.get.nc(nc, "latitude")[89:117] # 47.375:54.375
date <- ymd("19680915")
n <- which(date == seq(ymd("19650101"), ymd("19791231"), by = 1))
Eobs.p <- var.get.nc(nc, "rr", start = c(132, 89, n-7), count = c(57, 29, 15))
Eobs.p <- Eobs.p/10
close.nc(nc)
# E-OBS 1 deg
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/E-OBS_1950-2017_v16.0/rr_1deg_reg_19689808-1980922_7.5W_6.5E_47.5N_54.5N_v16.0.nc")
Eobs1deg.p <- var.get.nc(nc, "rr")/10
close.nc(nc)
# REGEN Long Term 40yr
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/FORTY/REGEN_V1.1_LongTermStns_19680908-19680922_1deg.nc")
# There was heavy rainfall on 15th Sept 1968 that caused floods in the Great Floods of 1968 in UK and France
REGEN40yr.p <- var.get.nc(nc, "p", start = c(173, 138, NA), count = c(15, 8, NA))
close.nc(nc)


# Spatial average to get the time series
p.ts <- apply(p, 3, mean, na.rm = T)
Eobs.ts <- apply(Eobs.p, 3, mean, na.rm = T)
Eobs1deg.ts <- apply(Eobs1deg.p, 3, mean, na.rm = T)
REGEN40yr.ts <- apply(REGEN40yr.p, 3, mean, na.rm = T)
yrange <- range(p.ts, Eobs.ts, Eobs1deg.ts)
n <- which(date == dates)
dates.ts <- dates[(n - 7):(n + 7)]
cols <- c('#e41a1cdd','#377eb8dd','#4daf4add','#984ea3dd','#ff7f00dd')
par(mar = c(6,5,0,0) + 0.1, bg = NA)
plot(dates.ts, p.ts, type = "l", xlab = "", ylab = "Average precipitation (mm/day)", 
     cex.lab = 2, cex.axis = 2, lwd = 6, las = 1, xaxt = "n", ylim = yrange, col = cols[1])
lines(dates.ts, REGEN40yr.ts, col = cols[2], lwd = 6)
lines(dates.ts, Eobs.ts, col = cols[3], lwd = 6)
lines(dates.ts, Eobs1deg.ts, col = cols[4], lwd = 6)
axis(side = 1, at = dates.ts, labels = F)
text(dates.ts, labels = paste(month(dates.ts, label = T), day(dates.ts)), 
     par("usr")[3], srt = 45, cex = 1.5, xpd = T, adj = c(1,1.2))
title(main = "", xlab = "Date", line = 5, cex.lab = 2)
legend("topright", legend = c("REGEN V1.1", "REGEN Long term 40yr", "E-OBS V16 0.25 deg", "E-OBS V16 1 deg"), 
       col = cols[1:4], lwd = 6, bty = "n", cex = 2)
# Try doing the same for the entire month of Dec 2002
date <- ymd("20021201")
dates <- seq(ymd("20021201"), ymd("20021231"), 1)

# REGEN
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/CCRC_V1.0_20021201-20021231_1deg_p.nc")
lon <- var.get.nc(nc, "lon")[169:210] 
lat <- var.get.nc(nc, "lat")[125:150]
p <- var.get.nc(nc, "p", start = c(169, 125, NA), count = c(211-169, 151-125, NA))
close.nc(nc)

# REGEN Long Term 40yr
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.0/FORTY/CCRC_V1.0_20021201-20021231_1deg_longterm40yr_p.nc")
lon <- var.get.nc(nc, "lon")[169:210] 
lat <- var.get.nc(nc, "lat")[125:150]
REGEN40yr.p <- var.get.nc(nc, "p", start = c(169, 125, NA), count = c(211-169, 151-125, NA))
close.nc(nc)

# EOBS
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/E-OBS_1950-2017_v16.0/rr_0.25deg_reg_Dec2002_v16.0.nc")
Eobs.lon <- var.get.nc(nc, "longitude")[116:280] # -11.5:29.5
Eobs.lat <- var.get.nc(nc, "latitude")[37:137] # 34.5:59.5
n <- which(date == seq(ymd("19950101"), ymd("20101231"), by = 1))
Eobs.p <- var.get.nc(nc, "rr", start = c(116, 37, NA), count = c(281-116, 138-37, NA))/10
close.nc(nc)
# GPCC-FDD
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/GPCC_FDD/full_data_daily_Dec2002.nc")
Gpcc.lon <- var.get.nc(nc, "lon")[169:210] 
Gpcc.lat <- var.get.nc(nc, "lat")[125:150]
Gpcc.p <- var.get.nc(nc, "p", start = c(169, 125, NA), count = c(211-169, 151-125, NA))
close.nc(nc)
# E-Obs 1 deg
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/E-OBS_1950-2017_v16.0/rr_1deg_reg_Dec2002_v16.0.nc")
Eobs1deg.lon <- var.get.nc(nc, "longitude")[30:71]
Eobs1deg.lat <- var.get.nc(nc, "latitude")[10:35]
Eobs1deg.p <- var.get.nc(nc, "rr", start = c(30, 10, NA), count = c(72-30, 36-10, NA))/10
close.nc(nc)
#E-Obs 0.5 deg
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/E-OBS_1950-2017_v16.0/rr_0.5deg_20021201-20021231_11.5W_29.5E_34.5N_59.5N.nc")
Eobs0.5deg.p <- var.get.nc(nc, "rr")/10
Eobs0.5deg.lon <- var.get.nc(nc, "longitude")
Eobs0.5deg.lat <- var.get.nc(nc, "latitude")
close.nc(nc)

# Spatial average to get the time series
p.ts <- apply(p, 3, mean, na.rm = T)
Eobs.ts <- apply(Eobs.p, 3, mean, na.rm = T)
Gpcc.ts <- apply(Gpcc.p, 3, mean, na.rm = T)
REGEN40yr.ts <- apply(REGEN40yr.p, 3, mean, na.rm = T)
Eobs1deg.ts <- apply(Eobs1deg.p, 3, mean, na.rm = T)
Eobs0.5deg.ts <- apply(Eobs0.5deg.p, 3, mean, na.rm = T)
yrange <- range(p.ts, Eobs.ts, Gpcc.ts, Eobs1deg.ts, Eobs0.5deg.ts)
dates.ts <- dates
cols <- c('#e41a1cdd','#377eb8dd','#4daf4add','#984ea3dd','#ff7f00dd') 
        #c('#1b9e77dd','#d95f02dd','#7570b3dd', '#e7298add', '#66a61edd')
par(mar = c(6,5,0,0) + 0.1, bg = NA)
plot(dates.ts, p.ts, type = "l", xlab = "", ylab = "Average precipitation (mm/day)", 
     cex.lab = 2, cex.axis = 2, lwd = 6, las = 1, xaxt = "n", ylim = yrange, col = cols[1])
lines(dates.ts, REGEN40yr.ts, col = cols[2], lwd = 6)
lines(dates.ts, Gpcc.ts, col = cols[3], lwd = 6)
lines(dates.ts, Eobs.ts, col = cols[4], lwd = 6)
#lines(dates.ts, Eobs1deg.ts, col = cols[5], lwd = 6)
#lines(dates.ts, Eobs0.5deg.ts, col = cols[5], lwd = 6)
axis(side = 1, at = dates.ts, labels = F)
text(dates.ts, labels = paste(month(dates.ts, label = T), day(dates.ts)), 
     par("usr")[3], srt = 45, cex = 1.5, xpd = T, adj = c(1,1.2))
title(main = "", xlab = "Date", line = 5, cex.lab = 2)
legend("topright", legend = c("REGEN V1.0", "REGEN Long term 40yr", "GPCC-FDD", "E-OBS V16 0.25 deg", "E-OBS V16 1 deg", "E-OBS V16 0.5 deg")[1:4], 
       col = cols[1:4], lwd = 6, bty = "n", cex = 2)
par(new = T, bg = NA)
map("world", xlim = c(-11.5, 29.5), ylim = c(34.5, 59.5))
# AUSTRALIA Instead
# Cyclone YASI 2011
library(RNetCDF)
library(lubridate)

nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_20110127-20110210_1deg.nc")
lon <- var.get.nc(nc, "lon")[292:335]
lat <- var.get.nc(nc, "lat")[46:79]
dates.ts <- seq(ymd("20110127"), ymd("20110210"), by = 1)
p <- var.get.nc(nc, "p", start = c(292, 46, NA), count = c(336-292, 80-46, NA))
close.nc(nc)

# Long term 40yr date
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/FORTY/REGEN_V1.1_LongTermStns_20110127-20110210_1deg.nc")
REGEN40yr.p <- var.get.nc(nc, "p", start = c(292, 46, NA), count = c(336-292, 80-46, NA))
close.nc(nc)

# Shift AWAP for same day comparison
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/AWAP/AWAP_PRCP_20110128-20110211_Yasi_0.5deg_remapcon2_land_noneg.nc")
awap.lon <- var.get.nc(nc, "lon")
awap.lat <- rev(var.get.nc(nc, "lat"))
awap.p <- var.get.nc(nc, "pre")[,length(awap.lat):1,]
close.nc(nc)

# Spatial average to get the time series
p.ts <- apply(p, 3, mean, na.rm = T)
awap.ts <- apply(awap.p, 3, mean, na.rm = T)
REGEN40yr.ts <- apply(REGEN40yr.p, 3, mean, na.rm = T)
yrange <- range(p.ts, awap.ts)
cols <- c('#1b9e77dd','#d95f02dd','#7570b3dd', '#e7298add', '#66a61edd')
par(mar = c(6,5,0,0) + 0.1, bg = NA)
plot(dates.ts, p.ts, type = "l", xlab = "", ylab = "Average precipitation (mm/day)", 
     cex.lab = 2, cex.axis = 2, lwd = 6, las = 1, xaxt = "n", ylim = yrange, col = cols[1])
lines(dates.ts, REGEN40yr.ts, col = cols[2], lwd = 6)
lines(dates.ts, awap.ts, col = cols[3], lwd = 6)
axis(side = 1, at = dates.ts, labels = F)
text(dates.ts, labels = paste(month(dates.ts, label = T), day(dates.ts)), 
     par("usr")[3], srt = 45, cex = 1.5, xpd = T, adj = c(1,1.2))
title(main = "", xlab = "Date", line = 5, cex.lab = 2)
legend("topleft", legend = c("REGEN V1.1", "REGEN Long term 40yr", "AWAP 0.5deg"), 
       col = cols[1:3], lwd = 6, bty = "n", cex = 2)
# APHRODITE

library(RNetCDF)
library(lubridate)

dates <- seq(ymd("19911025"), ymd("19911115"), by = 1)
which(dates == ymd("19911105"))

nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_19911025-19911115_1deg.nc")
lon <- var.get.nc(nc, "lon")[300:307] # 119.5 126.5
lat <- var.get.nc(nc, "lat")[96:110] # 5.5 19.5
p <- var.get.nc(nc, "p", start = c(300, 96, NA), count = c(308-300, 111-96, NA))
s <- var.get.nc(nc, "s", start = c(300, 96, NA), count = c(308-300, 111-96, NA))
close.nc(nc)

# Long Term 40yr data
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/FORTY/REGEN_V1.1_LongTermStns_19911025-19911115_1deg.nc")
REGEN40yr.p <- var.get.nc(nc, "p", start = c(300, 96, NA), count = c(308-300, 111-96, NA))
#REGEN40yr.s <- var.get.nc(nc, "s", start = c(300, 96, NA), count = c(308-300, 111-96, NA))
close.nc(nc)

nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/APHRODITE/APHRO_MA_025deg_V1101.19911025-19911115.nc")
Aphro.lon <- var.get.nc(nc, "longitude")[240:268] # 119.875 126.875
Aphro.lat <- var.get.nc(nc, "latitude")[81:140] # 5.125 19.875
Aphro.p <- var.get.nc(nc, "precip", start = c(240, 81, NA), count = c(269-240, 141-81, NA))
Aphro.s <- var.get.nc(nc, "rstn", start = c(240, 81, NA), count = c(269-240, 141-81, NA))
close.nc(nc)

find.Xtent.index <- function(Xtent, lon, lat) {
# Returns the x and y indices of vectors of longitude and latitude given a 
# raster::extent object created for example by raster:drawExtent()
  return(list(x = c(which.min(abs(Xtent[1] - lon)), 
                    which.min(abs(Xtent[2] - lon))),
              y = c(which.min(abs(Xtent[3] - lat)), 
                    which.min(abs(Xtent[4] - lat)))))
}

# Spatial average to get the time series
p.ts <- apply(p, 3, mean, na.rm = T)
aphro.ts <- apply(Aphro.p, 3, mean, na.rm = T)
REGEN40yr.ts <- apply(REGEN40yr.p, 3, mean, na.rm = T)
yrange <- range(p.ts, aphro.ts)
cols <- c('#1b9e77dd','#d95f02dd','#7570b3dd', '#e7298add', '#66a61edd')
par(mar = c(6,5,0,0) + 0.1, bg = NA)
plot(dates, p.ts, type = "l", xlab = "", ylab = "Average precipitation (mm/day)", 
     cex.lab = 2, cex.axis = 2, lwd = 6, las = 1, xaxt = "n", ylim = yrange, col = cols[1])
lines(dates, REGEN40yr.ts, col = cols[2], lwd = 6)
lines(dates, aphro.ts, col = cols[3], lwd = 6)
axis(side = 1, at = dates, labels = F)
text(dates, labels = paste(month(dates, label = T), day(dates)), 
     par("usr")[3], srt = 45, cex = 1.5, xpd = T, adj = c(1,1.2))
title(main = "", xlab = "Date", line = 5, cex.lab = 2)
legend("topright", legend = c("REGEN V1.1", "REGEN Long term 40yr","APHRODITE 0.25deg"), 
       col = cols, lwd = 6, bty = "n", cex = 2)
# Now do the same for US
library(RNetCDF)
library(lubridate)
library(raster)
library(maps)

dates <- seq(ymd("19780723"), ymd("19780806"), 1)

# Import REGEN data
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_19780723-19780806_1deg.nc")
lon <- var.get.nc(nc, "lon")[54:115]
lat <- var.get.nc(nc, "lat")[116:134]
p <- var.get.nc(nc, "p", start = c(54, 116, NA), count = c(116-54, 135-116, NA))
close.nc(nc)

# Long term 40yr data
# Import REGEN data
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/FORTY/REGEN_V1.1_LongTermStns_19780723-19780806_1deg.nc")
lon <- var.get.nc(nc, "lon")[54:115]
lat <- var.get.nc(nc, "lat")[116:134]
REGEN40yr.p <- var.get.nc(nc, "p", start = c(54, 116, NA), count = c(116-54, 135-116, NA))
close.nc(nc)

# Import CPC data
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/CPC_.25x.25_Daily_US_Unified_Precipitation/precip.V1.0.192780723-19780806.classic.nc")
# CPC longitudes go from 0 to 360, we want -180 to 180
cpc.lon <- var.get.nc(nc, "lon") - 360
cpc.lat <- var.get.nc(nc, "lat")
cpc.p <- var.get.nc(nc, "precip")
close.nc(nc)

# Convert cpc.p into a raster brick
cpc.p.raster <- flip(t(brick(cpc.p, xmn = min(cpc.lat), xmx = max(cpc.lat), ymn = min(cpc.lon), ymx = max(cpc.lon),
                             crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")), direction = "y")

# Below shapefile downloaded from https://www.census.gov/geo/maps-data/data/cbf/cbf_nation.html
# but the above wrld_simpl polygon will also work
uspolygon <- crop(shapefile("/home/z3289452/Downloads/cb_2016_us_nation_20m/cb_2016_us_nation_20m.cpg"),
                  extent(-125,-66,25,50))
# Transform the CRS of the uspolygon shapefile
uspolygon <- spTransform(uspolygon, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))

# Convert REGEN prcp into raster brick
p.raster <- flip(t(brick(p, xmn = min(lat), xmx = max(lat), ymn = min(lon), ymx = max(lon),
                  crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")), direction = "y")
# mask the REGEN data by US spatial coastline polygon
p.mask <- mask(p.raster, uspolygon)
# Crop the CPC data by the extend of REGEN data
cpc.p.crop <- crop(cpc.p.raster, extent(p.mask))

# Convert REGEN40yr prcp into raster brick
REGEN40yr.p.raster <- flip(t(brick(REGEN40yr.p, xmn = min(lat), xmx = max(lat), ymn = min(lon), ymx = max(lon),
                  crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")), direction = "y")
# mask the REGEN data by US spatial coastline polygon
REGEN40yr.p.mask <- mask(REGEN40yr.p.raster, uspolygon)


# Calculate the spatial means of each raster layer
cpc.ts <- cellStats(cpc.p.crop, 'mean')
# checking if the values look reasonable
# apply(cpc.p, 3, mean, na.rm = T)
p.ts <- cellStats(p.mask, 'mean')
REGEN40yr.p.ts <- cellStats(REGEN40yr.p.mask, 'mean')
yrange <- range(p.ts, cpc.ts)
n <- which(date == dates)
dates.ts <- dates
cols <- c('#1b9e77dd','#d95f02dd','#7570b3dd', '#e7298add', '#66a61edd')
par(mar = c(6,5.5,0,0) + 0.1, bg = NA)
plot(dates.ts, p.ts, type = "l", xlab = "", ylab = "", 
     cex.lab = 2, cex.axis = 2, lwd = 6, las = 1, xaxt = "n", ylim = yrange, col = cols[1])
lines(dates.ts, REGEN40yr.p.ts, col = cols[2], lwd = 6)
lines(dates.ts, cpc.ts, col = cols[3], lwd = 6)
axis(side = 1, at = dates.ts, labels = F)
text(dates.ts, labels = paste(month(dates.ts, label = T), day(dates.ts)), 
     par("usr")[3], srt = 45, cex = 1.5, xpd = T, adj = c(1,1.2))
title(main = "", xlab = "Date", line = 5, cex.lab = 2)
title(ylab = "Average precipitation (mm/day)", line = 4, cex.lab = 2)
legend("topleft", legend = c("REGEN V1.1", "REGEN Long term 40yr","CPC CONUS 0.25deg"), 
       col = cols[1:3], lwd = 6, bty = "n", cex = 2)

Data sparse regions

  • Siberia
  • South America
  • Africa
  • India

There will be a TS averaged over central africa where stations are variable

# Timeseries averaged over central Africa

library(RNetCDF)
library(plot3D)
library(maps)
library(lubridate)

# Seasonal TS
# done with cdo
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.0/ALL/CCRC_V1.0_1950-2013_1deg_AllStations_CentralAfrica_fldmean_seasmean.nc")
CentralAfrica.Seas.TS <- var.get.nc(nc, "p")
CentralAfrica.DJF.TS <- CentralAfrica.Seas.TS[seq(1, length(CentralAfrica.Seas.TS), 4)]
CentralAfrica.MAM.TS <- CentralAfrica.Seas.TS[seq(2, length(CentralAfrica.Seas.TS), 4)]
CentralAfrica.JJA.TS <- CentralAfrica.Seas.TS[seq(3, length(CentralAfrica.Seas.TS), 4)]
CentralAfrica.SON.TS <- CentralAfrica.Seas.TS[seq(4, length(CentralAfrica.Seas.TS), 4)]

par(mar = c(4,4.5,2,2) + 0.1, bg = NA)
cols <- c('#e41a1cbb','#377eb8bb','#4daf4abb','#984ea3bb')
years <- 1950:2013
quantile(c(CentralAfrica.DJF.TS[1:(length(CentralAfrica.DJF.TS) - 1)], CentralAfrica.MAM.TS, CentralAfrica.JJA.TS, CentralAfrica.SON.TS),
         probs = seq(0, 1, by = 0.05))
yrange <- c(1.15, 4.85)
plot(years, CentralAfrica.DJF.TS[1:(length(CentralAfrica.DJF.TS) - 1)], ylim = yrange, type = "n", 
     las = 1, cex = 2, cex.axis = 2, cex.lab = 2,
     xlab = "Year", ylab = "Average Seasonal Precipitation (mm/day)")
lines(years, CentralAfrica.DJF.TS[1:(length(CentralAfrica.DJF.TS) - 1)], lwd = 4, col = cols[1])
lines(years, CentralAfrica.MAM.TS, lwd = 4, col = cols[2])
lines(years, CentralAfrica.JJA.TS, lwd = 4, col = cols[3])
lines(years, CentralAfrica.SON.TS, lwd = 4, col = cols[4])
legend("topright", c("DJF", "MAM", "JJA", "SON"), lwd = 4, cex = 2, bty = "n", col = cols)
text(1960, yrange[2] - yrange[2]/50, "Central Africa", cex = 2)

There will also be some fields of a significant event in tropical central africa
2008 Benin floods

library(RNetCDF)
library(lubridate)
library(raster)
library(maps)
library(maptools)

# import Africa quality control mask and map it to see where to average
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_1950-2013_1deg_QualityMask.nc")
lon <- var.get.nc(nc, "lon")
lat <- var.get.nc(nc, "lat")
REGEN_africa_mask <- var.get.nc(nc, "p")
close.nc(nc)

# Create a mask over Benin
Benin_mask <- REGEN_africa_mask
Benin_mask.raster <- flip(t(raster(Benin_mask, xmn = min(lat), xmx = max(lat), ymn = min(lon), ymx = max(lon),
                                   crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")), direction = "y")
Benin <- map("world", "Benin", fill = T, plot = F)
Benin.sp <- map2SpatialPolygons(Benin, IDs = Benin$names, 
                                proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
Benin_mask.raster <- mask(Benin_mask.raster, Benin.sp)

# write out to netcdf
writeRaster(Benin_mask.raster,
            filename = "/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/BeningMask.nc",
            format = "CDF",
            varname = "p", varunit = "Boolean",
            longname = "A mask for the country Benin in West Africa",
            xname = "lon", yname = "lat", overwrite = T)

# using cdo invert the latitude of the raster and then setgrid to the same as central Africa REGEN and 
# then mask central africa REGEN with the mask
# average the masked Benin data spatially

Run the steps above to produce the file loaded below for the plots

# import the spatially averaged timeseries
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_1950-2013_1deg_Benin_fldmean.nc")
Benin.ts <- var.get.nc(nc, "p")
dates <- seq(ymd("19500101"), ymd("20131231"), 1)
close.nc(nc)

# let's make some plots
# first the entire time series
par(mar = c(2.5,4.5,0.2,0.2) + 0.1, bg = NA)
col <- "#386cb0"
fadedCol <- "#386cb088"
col2 <- "#d95f02"
fadedCol2 <- "#d95f0288"
plot(dates, Benin.ts, col = fadedCol, lwd = 2, xlab = "", ylab = "Daily precipitation rate (mm/day)", type = "l",
     cex.axis = 1.5, cex.lab = 2)
subdates <- function(year) {which(year(dates) == year)}
lines(dates[subdates(1957)], Benin.ts[subdates(1957)], col = col, lwd = 2)
lines(dates[subdates(1962)], Benin.ts[subdates(1962)], col = col, lwd = 2)
lines(dates[subdates(2008)], Benin.ts[subdates(2008)], col = col, lwd = 2)

# Now make the timeseries plot for just 1957 and 2008 for comparison
plot(dates[subdates(1957)], Benin.ts[subdates(1957)], col = col, lwd = 2, xlab = "", ylab = "Daily precipitation rate (mm/day)", type = "l",
     cex.axis = 1.5, cex.lab = 2)
xx <- c(dates[subdates(1957)], rev(dates[subdates(1957)]))
yy1957 <- c(Benin.ts[subdates(1957)], rep(0, length(dates[subdates(1957)])))
yy2008 <- c(Benin.ts[subdates(2008)][2:366], rep(0, length(dates[subdates(1957)])))
polygon(xx, yy1957, col = fadedCol, border = NA)
lines(dates[subdates(1957)], Benin.ts[subdates(2008)][2:366], col = col2, lwd = 2)
polygon(xx, yy2008, col = fadedCol2, border = NA)
legend("topright", legend = c("1957", "2008"), border = c(col, col2), fill = c(fadedCol, fadedCol2), cex = 2,
       bty = "n")

Temporal correlation

  • matrix of spatial correlations
  • taylor diagrams

CCRC vs GPCC correlation
Field Correlation

# Field correlations
library(RNetCDF)
library(lubridate)

nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_CPC_1979-2013_1deg_fldcor.nc")
CCRC_CPC.fldcor <- var.get.nc(nc, "p")
close.nc(nc)
CCRC_CPC.dates <- seq(ymd("19790101"), ymd("20131231"), 1)

nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_GPCC_1988-2013_1deg_fldcor.nc")
CCRC_GPCC.fldcor <- var.get.nc(nc, "p")
close.nc(nc)
CCRC_GPCC.dates <- seq(ymd("19880101"), ymd("20131231"), 1)

nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_LongTermStns_1950-2013_1deg_fldcor.nc")
REGEN_REGEN40yr.fldcor <- var.get.nc(nc, "p")
close.nc(nc)
REGEN_REGEN40yr.dates <- seq(ymd("19500101"), ymd("20131231"), 1)


par(mar = c(5,5,4,2)+0.1, bg= NA)
plot(CCRC_CPC.dates, CCRC_CPC.fldcor, type = 'n', xlab = "", ylab = "REGEN vs CPC Correlation",
     ylim = c(0.15,1), cex.axis = 2, cex.lab = 2)
lines(CCRC_CPC.dates, CCRC_CPC.fldcor, lwd = 2, col = "#1b9e77")

par(mar = c(5,5,4,2)+0.1, bg = NA)
plot(CCRC_GPCC.dates, CCRC_GPCC.fldcor, type = 'n', xlab = "", ylab = "REGEN vs GPCC Correlation",
     ylim = c(0.15,1), cex.axis = 2, cex.lab = 2)
lines(CCRC_GPCC.dates, CCRC_GPCC.fldcor, lwd = 2, col = '#d95f02')

par(mar = c(5,5,4,2)+0.1, bg = NA)
plot(REGEN_REGEN40yr.dates, REGEN_REGEN40yr.fldcor, type = 'n', xlab = "", ylab = "REGEN vs REGEN Long term 40yr  Correlation",
     ylim = c(0.15,1), cex.axis = 2, cex.lab = 1.5)
lines(REGEN_REGEN40yr.dates, REGEN_REGEN40yr.fldcor, lwd = 2, col = '#7570b3')

Temporal Correlation

library(plot3D)
# Temporal correlation 
nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_GPCC_1988-2013_1deg_timcor.nc")
CCRC_GPCC.timcor <- var.get.nc(nc, "p")
lon <- var.get.nc(nc, "lon")
lat <- var.get.nc(nc, "lat")
close.nc(nc)

quantile(CCRC_GPCC.timcor, probs = seq(0, 1, 0.05), na.rm = T)
brks <- seq(0, 1, length.out = 6)
zbrks <- brks
zbrks[1] <- min(brks, CCRC_GPCC.timcor, na.rm = T)
zbrks[length(brks)] <- max(brks, CCRC_GPCC.timcor, na.rm = T)
lbrks <- ifelse(brks == zbrks, brks, "")
cols <- c('#ffffcc','#a1dab4','#41b6c4','#2c7fb8','#253494')

par(mar = c(0,0,0,0) + 0.1, bg = NA)
layout(mat = matrix(c(1,2), nrow = 2, byrow = T), heights = c(1,0.15))
image(lon, lat, CCRC_GPCC.timcor, xlab = "", ylab = "", breaks = zbrks, 
      col = cols, xaxt = "n", yaxt = "n", cex.lab = 1.5, cex.axis = 1.5, ylim = c(-60, 90))
map("world", interior = F, add = T)
par(mar = c(0,8,0,8) + 0.1)
colkey(col = cols, side = 1, clim = range(brks), width = 5, at = brks, 
       labels = paste(lbrks), cex.axis = 1.5)
title("REGEN V1.1 vs GPCC-FDD temporal correlation 1988-2013", line = -1.5, cex.main = 1.5)

Temporal Correlation

# Temporal correlation 

nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_CPC_1979-2013_1deg_timcor.nc")
CCRC_CPC.timcor <- var.get.nc(nc, "p")
lon <- var.get.nc(nc, "lon")
lat <- var.get.nc(nc, "lat")
close.nc(nc)

quantile(CCRC_CPC.timcor, probs = seq(0, 1, 0.05), na.rm = T)
brks <- c(-1, seq(0, 1, length.out = 6))
zbrks <- brks
zbrks[1] <- min(brks, CCRC_CPC.timcor, na.rm = T)
zbrks[length(brks)] <- max(brks, CCRC_CPC.timcor, na.rm = T)
lbrks <- ifelse(brks == zbrks, brks, "")
cols <- c('#fa9fb5','#ffffcc','#a1dab4','#41b6c4','#2c7fb8','#253494')

par(mar = c(0,0,0,0) + 0.1, bg = NA)
layout(mat = matrix(c(1,2), nrow = 2, byrow = T), heights = c(1,0.15))
image(lon, lat, CCRC_CPC.timcor, xlab = "", ylab = "", breaks = zbrks, 
      col = cols, xaxt = "n", yaxt = "n", cex.lab = 1.5, cex.axis = 1.5, ylim = c(-60, 90))
map("world", interior = F, add = T)
par(mar = c(0,8,0,8) + 0.1)
colkey(col = cols, side = 1, clim = c(1,7), width = 5, at = 1:7, 
       labels = paste(lbrks), cex.axis = 1.5)
title("REGEN V1.1 vs CPC_Glb temporal correlation 1988-2013", line = -1.5, cex.main = 1.5)

Field correlation based on only those grids where ALL and Long term quality masks overlap: 

# Temporal correlation 
library(maps)
library(plot3D)

nc <- open.nc("/srv/ccrc/data11/z3289452/DailyPRCP_netCDFs/CCRCGlobalGriddedDailyPrecip/VERSION1.1/ALL/REGEN_V1.1_ALLStns_LongTermStns_1950-2013_1deg_timcor.nc")
REGEN_REGEN40yr.timcor <- var.get.nc(nc, "p")
lon <- var.get.nc(nc, "lon")
lat <- var.get.nc(nc, "lat")
close.nc(nc)

quantile(REGEN_REGEN40yr.timcor, probs = seq(0, 1, 0.05), na.rm = T)
brks <- c(-1, seq(0, 1, length.out = 6))
zbrks <- brks
zbrks[1] <- min(brks, REGEN_REGEN40yr.timcor, na.rm = T)
zbrks[length(brks)] <- max(brks, REGEN_REGEN40yr.timcor, na.rm = T)
lbrks <- ifelse(brks == zbrks, brks, "")
cols <- c('#fa9fb5','#ffffcc','#a1dab4','#41b6c4','#2c7fb8','#253494')

par(mar = c(0,0,0,0) + 0.1, bg = NA)
layout(mat = matrix(c(1,2), nrow = 2, byrow = T), heights = c(1,0.15))
image(lon, lat, REGEN_REGEN40yr.timcor, xlab = "", ylab = "", breaks = zbrks, 
      col = cols, xaxt = "n", yaxt = "n", cex.lab = 1.5, cex.axis = 1.5, ylim = c(-60, 90))
map("world", interior = F, add = T)
par(mar = c(0,8,0,8) + 0.1)
colkey(col = cols, side = 1, clim = c(1,7), width = 5, at = 1:7, 
       labels = paste(lbrks), cex.axis = 1.5)
title("REGEN vs REGEN Long Term 40yr temporal correlation 1988-2013", line = -1.5, cex.main = 1.5)

Temporal correlation based on all global land grids:
Temporal correlation based on only those grids where ALL and Long term quality masks overlap: No need to correlations with monthly datasets