Intensification of daily light, moderate and extreme rainfall across Australia - Plots and code

This document consists of all the code to create the plots in the manuscript. In most places the .RData data files have been provided. Descriptions of how the .Rdata files were generated have been provided as well. In some instances the code to generate the .Rdata files have also been provided.

Frequency changes in wet days

First we set up our working directory and initial variables and libs.

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

#  setwd 
dir = "/srv/ccrc/data45/z3289452/PhD_repository/"
setwd("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/")

# Data import
# Before importing data, both AWAP and GPCC_FDD were regridded onto the same 1X1 degree grid using cdo remapcon2
dates <- seq(as.Date("1958-01-01"), as.Date("2013-12-31"),1)
lon <- seq(112, by=0.5, length.out = 88)
lat <- rev(seq(-10, by=-0.5, length.out = 68))

** Initial data import from .nc files** Make sure the libs and initial variables are set up. **AWAP**

# Import AWAP data
nc.dir <- "/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/AWAP/" 
nc.name <- "AWAP_PRCP_1900-2015_0.5deg_remapcon2_land_noneg.nc"
nc <- open.nc(paste0(nc.dir, nc.name))
PRCP_AWAP <- var.get.nc(nc, "pre", start=c(1,1,21185), count=c(88,68,length(dates)))
close.nc(nc); rm(nc); gc()
PRCP_AWAP <- PRCP_AWAP[,dim(PRCP_AWAP)[2]:1,]
#Limit analysis to 2013 because 2014-15 data for Natgrid and Barnes not available
PRCP_AWAP_1 <- PRCP_AWAP[,,1:10227] # 10227->1985-12-31 10592->1986-12-31
PRCP_AWAP_2 <- PRCP_AWAP[,,10228:dim(PRCP_AWAP)[3]] # 20454->2013-12-31

AWAP wet days data This is not needed for wet days plots.

# remove precip under 1mm
WET_AWAP_1 <- PRCP_AWAP_1
WET_AWAP_1[which(WET_AWAP_1 < 1.0)] <- NA
WET_AWAP_2 <- PRCP_AWAP_2
WET_AWAP_2[which(WET_AWAP_2 < 1.0)] <- NA

BOA (BRNS)

# Import Brns data
nc.dir <- "/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/NatgridBarnes/"
nc.name <- "Prcp_daily_barnesgrid_0.5x0.5deg_1958-2013_Req1_land.nc"
nc <- open.nc(paste0(nc.dir,nc.name))
PRCP_BRNS <- var.get.nc(nc, "Prcp")
close.nc(nc)
# Natgrid and Barnes wet prcp for two periods
PRCP_BRNS_1 <- PRCP_BRNS[,,1:10227] # 10227->1985-12-31 10592->1986-12-31
PRCP_BRNS_2 <- PRCP_BRNS[,,10228:dim(PRCP_BRNS)[3]] # 20454->2013-12-31

** BOA wet days data** This is not needed for wet days plots.

# remove precip under 1mm
WET_BRNS_1 <- PRCP_BRNS_1
WET_BRNS_1[which(WET_BRNS_1 < 1.0)] <- NA
WET_BRNS_2 <- PRCP_BRNS_2
WET_BRNS_2[which(WET_BRNS_2 < 1.0)] <- NA

** NNI (NATN)**

# Import Natn data
nc.dir <- "/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/NatgridBarnes/"
nc.name <- "Prcp_daily_natgrid_0.5x0.5deg_1958-2013_land.nc"
nc <- open.nc(paste0(nc.dir,nc.name))
PRCP_NATN <- var.get.nc(nc, "Prcp")
close.nc(nc)
PRCP_NATN_1 <- PRCP_NATN[,,1:10227] # 10227->1985-12-31 10592->1986-12-31
PRCP_NATN_2 <- PRCP_NATN[,,10228:dim(PRCP_NATN)[3]] # 20454->2013-12-31

** NNI wet days data** This is not needed for wet days plots.

# remove precip under 1mm
WET_NATN_1 <- PRCP_NATN_1
WET_NATN_1[which(WET_NATN_1 < 1.0)] <- NA
WET_NATN_2 <- PRCP_NATN_2
WET_NATN_2[which(WET_NATN_2 < 1.0)] <- NA

Load masks for AWAP, BOA and NNI

# Load AWAP BRNS and NATN masks
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/AWAP_BRNS_NATN_mask.RData")
# Get rid of Giles
AWAP_mask[31:36,36:40] <- 0 
BRNS_mask[31:36,36:40] <- 0 

# Import and land sea mask from a netcdf file
nc.dir <- "/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/"
LSM.nc <- open.nc(paste0(nc.dir, "AWAP_Land-Sea-Mask_0.5deg_latflipped.nc"))
LSM <- var.get.nc(LSM.nc, "LSM")
rm(nc.dir, LSM.nc); gc()

#Create missing masks for plotting
#AWAP
AWAP.missing <- AWAP_mask
AWAP.missing[is.na(AWAP.missing)] <- 0
AWAP.missing[!LSM] <- NA
AWAP.missing[which(AWAP.missing == 1)] <- NA
AWAP.missing <- !AWAP.missing
#BRNS
BRNS.missing <- BRNS_mask
BRNS.missing[is.na(BRNS.missing)] <- 0
BRNS.missing[!LSM] <- NA
BRNS.missing[which(BRNS.missing == 1)] <- NA
BRNS.missing <- !BRNS.missing

Function to select colours given breaks.

create.col.list <- function(brks, null.point) {
  require(graphics)
  RdOl <- colorRampPalette(c("darkred","orchid","yellow"))
  OlBl <- colorRampPalette(c("lightgreen", "slateblue", "blue4"))
  num.red <- which.min(abs(brks - null.point))-1
  num.blue <- (length(brks)-1)-num.red
  if(num.red > num.blue){
    return(c(RdOl(num.red),OlBl(num.red)[1:num.blue]))
  } else {
    return(c(RdOl(num.blue)[(num.blue-num.red+1):num.blue],OlBl(num.blue)))
  }
}

AWAP, BOA and NNI annual change in frequency of wet days data

# Either run the calc.wetratiodiff function or load prerun and saved data below
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/AWAP_FDD_quantile_analysis_wetratiodiff_updated.RData")

The above data was created using the following function.

# Plot difference in proportions of wet days between the two periods
calc.wetratio <- function (PRCP_1, PRCP_2, dates1, dates2, mask) {
  WetRatio_1 <- array(dim=dim(PRCP_1)[1:2])
  WetRatio_2 <- array(dim=dim(PRCP_1)[1:2])
  for(i in 1:length(lon)){
    for(j in 1:length(lat)){
      WetRatio_1[i,j] <- length(which(PRCP_1[i,j,] >= 1.0))/length(dates1)
      WetRatio_2[i,j] <- length(which(PRCP_2[i,j,] >= 1.0))/length(dates2)
      if(WetRatio_1[i,j] == 0 & WetRatio_2[i,j] == 0) {
        WetRatio_1[i,j] <- NA
        WetRatio_2[i,j] <- NA
      }
    }
  }
  WetRatio_2 <- WetRatio_2*100
  WetRatio_1 <- WetRatio_1*100
  WetRatio_diff <- WetRatio_2 - WetRatio_1
  WetRatio_diff[which(mask==0)]<-NA
  invisible(WetRatio_diff)
}

# get the dates for the two periods
index.split <- which(dates == "1985-12-31")
dates.1 <- dates[1:index.split]
dates.2 <- dates[(index.split + 1): length(dates)]

Evaluate the function calc.wetratio for AWAP, BOA and NNI.

AWAP.wetratiodiff <- calc.wetratio(PRCP_AWAP_1, PRCP_AWAP_2, dates.1, dates.2, AWAP_mask)
# Memory management
rm(PRCP_AWAP, PRCP_AWAP_1, PRCP_AWAP_2); gc()
BRNS.wetratiodiff <- calc.wetratio(PRCP_BRNS_1, PRCP_BRNS_2, dates.1, dates.2, BRNS_mask)
# Memory management
rm(PRCP_BRNS, PRCP_BRNS_1, PRCP_BRNS_2); gc()
NATN.wetratiodiff <- calc.wetratio(PRCP_NATN_1, PRCP_NATN_2, dates.1, dates.2, NATN_mask)
# Memory management
rm(PRCP_NATN, PRCP_NATN_1, PRCP_NATN_2); gc()

Here PRCP_AWAP_1/2 are raw data arrays with dimensions c(num(lon), num(lat), num(dates.1/2)) (see initial data import code). Save .Rdata file with AWAP/BRNS/NATN.calc.wetratio

save(AWAP.wetratiodiff, BRNS.wetratiodiff, NATN.wetratiodiff, file = "RData/AWAP_FDD_quantile_analysis_wetratiodiff_updated.RData")

AWAP, BOA and NNI seasonal change in frequency of wet days AWAP Seasonal data import.

nc.dir <- "/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/AWAP/" 
nc.name.pre <- "AWAP_PRCP_1958-2013_0.5deg_remapcon2_land_noneg_"
season.names <- c("DJF","MAM","JJA","SON")
nc.name <- sapply(season.names, function (x) paste0(nc.name.pre,x,".nc"))
PRCP_AWAP_Season <- list()
nc <- open.nc(paste0(nc.dir,"AWAP_PRCP_1958-2014_0.5deg_remapcon2_land_noneg_DJF.nc"))
PRCP_AWAP_Season[[1]] <- var.get.nc(nc,"pre")
close.nc(nc);rm(nc);gc()
PRCP_AWAP_Season[[1]] <- PRCP_AWAP_Season[[1]][,dim(PRCP_AWAP_Season[[1]])[2]:1,]
# remove the first two months (1958 Jan/Feb) and the last month (2014 Dec)
PRCP_AWAP_Season[[1]] <- PRCP_AWAP_Season[[1]][,,60:5113]
for (i in 2:length(season.names)) {
  nc <- open.nc(paste0(nc.dir, nc.name[i]))
  PRCP_AWAP_Season[[i]] <- var.get.nc(nc, "pre")
  close.nc(nc); rm(nc); gc()
  PRCP_AWAP_Season[[i]] <- PRCP_AWAP_Season[[i]][,dim(PRCP_AWAP_Season[[i]])[2]:1,]
}
rm(i,nc.name.pre,nc.name,nc.dir)
names(PRCP_AWAP_Season) <- season.names

dates.season.AWAP <- list()
dates.djf <- seq(from=as.Date("1958-03-01"), to=as.Date("2014-02-28"),1)
dates.season.AWAP[[1]] <- dates.djf[which(month(dates.djf) %in% c(12,1,2))]; rm(dates.djf)
dates.season.AWAP[[2]] <- dates[which(month(dates) %in% c(3,4,5))]
dates.season.AWAP[[3]] <- dates[which(month(dates) %in% c(6,7,8))]
dates.season.AWAP[[4]] <- dates[which(month(dates) %in% c(9,10,11))]
names(dates.season.AWAP) <- season.names
dates.season.AWAP.1 <- list()
dates.season.AWAP.2 <- list()
index.split <- which(dates.season.AWAP[[1]] == "1986-02-28")
dates.season.AWAP.1[[1]] <- dates.season.AWAP[[1]][1:index.split]
dates.season.AWAP.2[[1]] <- dates.season.AWAP[[1]][(index.split + 1):length(dates.season.AWAP[[1]])]
for (s in 2:4) {
  index.split <- max(which(year(dates.season.AWAP[[s]]) == 1985))
  dates.season.AWAP.1[[s]] <- dates.season.AWAP[[s]][1:index.split]
  dates.season.AWAP.2[[s]] <- dates.season.AWAP[[s]][(index.split + 1):length(dates.season.AWAP[[s]])]
}

WET_AWAP_Season_1 <- list()
WET_AWAP_Season_2 <- list()
# do DJF separately because it is complicated
index.split <- which(dates.season.AWAP[[1]] == "1986-02-28")
WET_AWAP_Season_1[[1]] <- PRCP_AWAP_Season[[1]][,,1:index.split]
WET_AWAP_Season_2[[1]] <- PRCP_AWAP_Season[[1]][,,(index.split + 1):dim(PRCP_AWAP_Season[[1]])[3]]
WET_AWAP_Season_1[[1]][which(WET_AWAP_Season_1[[1]] < 1)] <- NA
WET_AWAP_Season_2[[1]][which(WET_AWAP_Season_2[[1]] < 1)] <- NA
for (i in 2:length(season.names)) {
  index.split <- max(which(year(dates.season.AWAP[[i]])==1985))
  WET_AWAP_Season_1[[i]] <- PRCP_AWAP_Season[[i]][,,1:index.split]
  WET_AWAP_Season_2[[i]] <- PRCP_AWAP_Season[[i]][,,(index.split+1):dim(PRCP_AWAP_Season[[i]])[3]]
  WET_AWAP_Season_1[[i]][which(WET_AWAP_Season_1[[i]] < 1)] <- NA
  WET_AWAP_Season_2[[i]][which(WET_AWAP_Season_2[[i]] < 1)] <- NA
}
names(WET_AWAP_Season_1) <- season.names
names(WET_AWAP_Season_2) <- season.names

rm(i,index.split,nc.dir,PRCP_AWAP_Season);gc()
# save RData image of seasonal data import
save(WET_AWAP_Season_1, WET_AWAP_Season_2, dates.season.AWAP.1, dates.season.AWAP.2, lon, lat, season.names,
     file="RData/AWAP_FDD_quantile_analysis_AWAPseasonaldataimport.RData")

Load .Rdata file with pre-run seasonal data import ^

load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/AWAP_FDD_quantile_analysis_AWAPseasonaldataimport.RData")

Evaluate calc.wetratio for each season in a loop

AWAP.wetratiodiff.Seas <- list()
for (s in 1:4) {
  writeLines(paste("Season",season.names[[s]]))
  AWAP.wetratiodiff.Seas[[s]] <- calc.wetratio(WET_AWAP_Season_1[[s]], WET_AWAP_Season_2[[s]], 
                                               dates.season.AWAP.1[[s]], dates.season.AWAP.2[[s]], AWAP_mask)
}

BOA (BRNS) Seasonal data import.

nc.dir <- "/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/NatgridBarnes/" 
nc.name.pre <- "Prcp_daily_barnesgrid_0.5x0.5deg_1958-2013_Req1_land_"
season.names <- c("DJF","MAM","JJA","SON")
nc.name <- sapply(season.names, function (x) paste0(nc.name.pre,x,".nc"))
PRCP_BRNS_Season <- list()
nc <- open.nc(paste0(nc.dir,nc.name[1]))
PRCP_BRNS_Season[[1]] <- var.get.nc(nc,"Prcp")
close.nc(nc);rm(nc);gc()
# remove the first two months (1958 Jan/Feb) and the last month (2013 Dec)
# Since we don't have 2014 data, our period 2 will be 1 season shorter than period 1
PRCP_BRNS_Season[[1]] <- PRCP_BRNS_Season[[1]][,,60:5023]
for (i in 2:length(season.names)) {
  nc <- open.nc(paste0(nc.dir, nc.name[i]))
  PRCP_BRNS_Season[[i]] <- var.get.nc(nc, "Prcp")
  close.nc(nc); rm(nc); gc()
}
rm(i,nc.name.pre,nc.name,nc.dir)
names(PRCP_BRNS_Season) <- season.names


dates.season.BRNS <- list()
dates.djf <- seq(from=as.Date("1958-03-01"), to=as.Date("2013-02-28"),1)
dates.season.BRNS[[1]] <- dates.djf[which(month(dates.djf) %in% c(12,1,2))]; rm(dates.djf)
dates.season.BRNS[[2]] <- dates[which(month(dates) %in% c(3,4,5))]
dates.season.BRNS[[3]] <- dates[which(month(dates) %in% c(6,7,8))]
dates.season.BRNS[[4]] <- dates[which(month(dates) %in% c(9,10,11))]
names(dates.season.BRNS) <- season.names
dates.season.BRNS.1 <- list()
dates.season.BRNS.2 <- list()
index.split <- which(dates.season.BRNS[[1]]=="1986-02-28")
dates.season.BRNS.1[[1]] <- dates.season.BRNS$DJF[1:index.split]
dates.season.BRNS.2[[1]] <- dates.season.BRNS$DJF[(index.split+1):length(dates.season.BRNS[[1]])]
for (s in 2:4) {
  index.split <- max(which(year(dates.season.BRNS[[s]]) == 1985))
  dates.season.BRNS.1[[s]] <- dates.season.BRNS[[s]][1:index.split]
  dates.season.BRNS.2[[s]] <- dates.season.BRNS[[s]][(index.split+1):length(dates.season.BRNS[[s]])]
}

WET_BRNS_Season_1 <- list()
WET_BRNS_Season_2 <- list()
# do DJF separately because it is complicated
index.split <- which(dates.season.BRNS[[1]]=="1986-02-28")
WET_BRNS_Season_1[[1]] <- PRCP_BRNS_Season[[1]][,,1:index.split]
WET_BRNS_Season_2[[1]] <- PRCP_BRNS_Season[[1]][,,(index.split+1):dim(PRCP_BRNS_Season[[1]])[3]]
WET_BRNS_Season_1[[1]][which(WET_BRNS_Season_1[[1]] < 1)] <- NA
WET_BRNS_Season_2[[1]][which(WET_BRNS_Season_2[[1]] < 1)] <- NA
for (i in 2:length(season.names)) {
  index.split <- max(which(year(dates.season.BRNS[[i]])==1985))
  WET_BRNS_Season_1[[i]] <- PRCP_BRNS_Season[[i]][,,1:index.split]
  WET_BRNS_Season_2[[i]] <- PRCP_BRNS_Season[[i]][,,(index.split+1):dim(PRCP_BRNS_Season[[i]])[3]]
  WET_BRNS_Season_1[[i]][which(WET_BRNS_Season_1[[i]] < 1)] <- NA
  WET_BRNS_Season_2[[i]][which(WET_BRNS_Season_2[[i]] < 1)] <- NA
}
names(WET_BRNS_Season_1) <- season.names
names(WET_BRNS_Season_2) <- season.names

rm(i,index.split,nc.dir,PRCP_BRNS_Season);gc()
# save RData image of seasonal data import
save(WET_BRNS_Season_1, WET_BRNS_Season_2, dates.season.BRNS.1, dates.season.BRNS.2, lon, lat, season.names,
     file="RData/AWAP_FDD_quantile_analysis_BRNSseasonaldataimport.RData")

Load .Rdata file with pre-run seasonal data import ^

load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/AWAP_FDD_quantile_analysis_BRNSseasonaldataimport.RData")

Evaluate calc.wetratio for each season in a loop

BRNS.wetratiodiff.Seas <- list()
for (s in 1:4) {
  writeLines(paste("Season",season.names[[s]]))
  BRNS.wetratiodiff.Seas[[s]] <- calc.wetratio(WET_BRNS_Season_1[[s]], WET_BRNS_Season_2[[s]], 
                                               dates.season.BRNS.1[[s]], dates.season.BRNS.2[[s]], BRNS_mask)
}

NNI (NATN) Seasonal data import.

nc.dir <- "/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/NatgridBarnes/" 
nc.name.pre <- "Prcp_daily_natgrid_0.5x0.5deg_1958-2013_land_"
season.names <- c("DJF","MAM","JJA","SON")
nc.name <- sapply(season.names, function (x) paste0(nc.name.pre,x,".nc"))
PRCP_NATN_Season <- list()
nc <- open.nc(paste0(nc.dir,nc.name[1]))
PRCP_NATN_Season[[1]] <- var.get.nc(nc,"Prcp")
close.nc(nc);rm(nc);gc()
# remove the first two months (1958 Jan/Feb) and the last month (2013 Dec)
# Since we don't have 2014 data, our period 2 will be 1 season shorter than period 1
PRCP_NATN_Season[[1]] <- PRCP_NATN_Season[[1]][,,60:5023]
for (i in 2:length(season.names)) {
  nc <- open.nc(paste0(nc.dir, nc.name[i]))
  PRCP_NATN_Season[[i]] <- var.get.nc(nc, "Prcp")
  close.nc(nc); rm(nc); gc()
}
rm(i,nc.name.pre,nc.name,nc.dir)
names(PRCP_NATN_Season) <- season.names


dates.season.NATN <- list()
dates.djf <- seq(from=as.Date("1958-03-01"), to=as.Date("2013-02-28"),1)
dates.season.NATN[[1]] <- dates.djf[which(month(dates.djf) %in% c(12,1,2))]; rm(dates.djf)
dates.season.NATN[[2]] <- dates[which(month(dates) %in% c(3,4,5))]
dates.season.NATN[[3]] <- dates[which(month(dates) %in% c(6,7,8))]
dates.season.NATN[[4]] <- dates[which(month(dates) %in% c(9,10,11))]
names(dates.season.NATN) <- season.names
dates.season.NATN.1 <- list()
dates.season.NATN.2 <- list()
index.split <- which(dates.season.NATN[[1]]=="1986-02-28")
dates.season.NATN.1[[1]] <- dates.season.NATN$DJF[1:index.split]
dates.season.NATN.2[[1]] <- dates.season.NATN$DJF[(index.split+1):length(dates.season.NATN[[1]])]
for (s in 2:4) {
  index.split <- max(which(year(dates.season.NATN[[s]]) == 1985))
  dates.season.NATN.1[[s]] <- dates.season.NATN[[s]][1:index.split]
  dates.season.NATN.2[[s]] <- dates.season.NATN[[s]][(index.split+1):length(dates.season.NATN[[s]])]
}

WET_NATN_Season_1 <- list()
WET_NATN_Season_2 <- list()
# do DJF separately because it is complicated
index.split <- which(dates.season.NATN[[1]] == "1986-02-28")
WET_NATN_Season_1[[1]] <- PRCP_NATN_Season[[1]][,,1:index.split]
WET_NATN_Season_2[[1]] <- PRCP_NATN_Season[[1]][,,(index.split+1):dim(PRCP_NATN_Season[[1]])[3]]
WET_NATN_Season_1[[1]][which(WET_NATN_Season_1[[1]] < 1)] <- NA
WET_NATN_Season_2[[1]][which(WET_NATN_Season_2[[1]] < 1)] <- NA
for (i in 2:length(season.names)) {
  index.split <- max(which(year(dates.season.NATN[[i]])==1985))
  WET_NATN_Season_1[[i]] <- PRCP_NATN_Season[[i]][,,1:index.split]
  WET_NATN_Season_2[[i]] <- PRCP_NATN_Season[[i]][,,(index.split+1):dim(PRCP_NATN_Season[[i]])[3]]
  WET_NATN_Season_1[[i]][which(WET_NATN_Season_1[[i]] < 1)] <- NA
  WET_NATN_Season_2[[i]][which(WET_NATN_Season_2[[i]] < 1)] <- NA
}
names(WET_NATN_Season_1) <- season.names
names(WET_NATN_Season_2) <- season.names

#load the mask
load("RData/AWAP_BRNS_NATN_mask.RData")

rm(i,index.split,nc.dir,PRCP_NATN_Season);gc()
# save RData image of seasonal data import
save(WET_NATN_Season_1, WET_NATN_Season_2, dates.season.NATN.1, dates.season.NATN.2, lon, lat, season.names,
     file="RData/AWAP_FDD_quantile_analysis_NATNseasonaldataimport.RData")

Load .Rdata file with pre-run seasonal data import ^

load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/AWAP_FDD_quantile_analysis_NATNseasonaldataimport.RData")

Evaluate calc.wetratio for each season in a loop

NATN.wetratiodiff.Seas <- list()
for (s in 1:4) {
  writeLines(paste("Season",season.names[[s]]))
  NATN.wetratiodiff.Seas[[s]] <- calc.wetratio(WET_NATN_Season_1[[s]], WET_NATN_Season_2[[s]], 
                                               dates.season.NATN.1[[s]], dates.season.NATN.2[[s]], NATN_mask)
}

Create the panel plot with both annual and seasonal wet days data

plot.wetratiodiff function This function creates a map plot of the wet ratio difference data

plot.wetratiodiff <- function (wetratiodiff, brks, col, zrange, mask) {
  zbrks <- brks; zbrks[1] <- zrange[1]; zbrks[length(zbrks)] <- zrange[2]
  image(lon,lat, wetratiodiff, breaks=zbrks, col=col, xaxt='n', yaxt='n', xlab="", ylab="")
  image(lon,lat, mask, col="gray", add = T)
  map("world", "Australia", add=T)
}

Before we go ahead make sure the masks are ready (run masks code chunk).

Plot NNI panel

# first we set up the layout
layout(mat = matrix(c(1,1,2:5,6,6), nrow = 4, ncol = 2, byrow = T), heights = c(0.35, 0.3, 0.3, 0.05), TRUE)

# determine brks and col for plotting
zrange <- range(AWAP.wetratiodiff.Seas, BRNS.wetratiodiff.Seas, NATN.wetratiodiff.Seas, na.rm = T)
brks <- seq(-2, 2, length.out = 9)
        #seq(-2,2,0.2)
zbrks <- brks; zbrks[1] <- zrange[1]; zbrks[length(zbrks)] <- zrange[2]
col <- c('#8c510a','#bf812d','#dfc27d','#f6e8c3','#c7eae5','#80cdc1','#35978f','#01665e')
      #create.col.list(brks, 0)


par(mar=c(1,9,1,9), oma=c(1,5,0.5,1.5))
plot.wetratiodiff(NATN.wetratiodiff, brks, col, zrange, AWAP.missing)
par(mar=c(0,0,0,0))
for (s in 1:4) {
  plot.wetratiodiff(NATN.wetratiodiff.Seas[[s]], brks, col, zrange, AWAP.missing)
  mtext(season.names[[s]], side = 3,line=-2, adj = 0.1)
}
colkey(col=col, side=1, clim=range(brks), width=8, at=brks, labels = paste(round(zbrks,1)))
mtext("(a)", side = 3, line = -3, adj = -0.12, cex = 1.7, outer = T)
mtext("(b)", side = 3, line = -25, adj = -0.12, cex = 1.7, outer = T)

Plot AWAP panel

# first we set up the layout
layout(mat = matrix(c(1,1,2:5,6,6), nrow = 4, ncol = 2, byrow = T), heights = c(0.35, 0.3, 0.3, 0.05), TRUE)

# determine brks and col for plotting
zrange <- range(AWAP.wetratiodiff.Seas, BRNS.wetratiodiff.Seas, NATN.wetratiodiff.Seas, na.rm = T)
brks <- seq(-2, 2, length.out = 9)
        #seq(-2,2,0.2)
zbrks <- brks; zbrks[1] <- zrange[1]; zbrks[length(zbrks)] <- zrange[2]
col <- c('#8c510a','#bf812d','#dfc27d','#f6e8c3','#c7eae5','#80cdc1','#35978f','#01665e')
      #create.col.list(brks, 0)


par(mar=c(1,9,1,9), oma=c(1,5,0.5,1.5))
plot.wetratiodiff(AWAP.wetratiodiff, brks, col, zrange, AWAP.missing)
par(mar=c(0,0,0,0))
for (s in 1:4) {
  plot.wetratiodiff(AWAP.wetratiodiff.Seas[[s]], brks, col, zrange, AWAP.missing)
  mtext(season.names[[s]], side = 3,line=-2, adj = 0.1)
}
colkey(col=col, side=1, clim=range(brks), width=8, at=brks, labels = paste(round(zbrks,1)))
mtext("(a)", side = 3, line = -3, adj = -0.12, cex = 1.7, outer = T)
mtext("(b)", side = 3, line = -25, adj = -0.12, cex = 1.7, outer = T)

Plot BOA panel

# first we set up the layout
layout(mat = matrix(c(1,1,2:5,6,6), nrow = 4, ncol = 2, byrow = T), heights = c(0.35, 0.3, 0.3, 0.05), TRUE)

# determine brks and col for plotting
zrange <- range(AWAP.wetratiodiff.Seas, BRNS.wetratiodiff.Seas, NATN.wetratiodiff.Seas, na.rm = T)
brks <- seq(-2, 2, length.out = 9)
        #seq(-2,2,0.2)
zbrks <- brks; zbrks[1] <- zrange[1]; zbrks[length(zbrks)] <- zrange[2]
col <- c('#8c510a','#bf812d','#dfc27d','#f6e8c3','#c7eae5','#80cdc1','#35978f','#01665e')
      #create.col.list(brks, 0)


par(mar=c(1,9,1,9), oma=c(1,5,0.5,1.5))
plot.wetratiodiff(BRNS.wetratiodiff, brks, col, zrange, BRNS.missing)
par(mar=c(0,0,0,0))
for (s in 1:4) {
  plot.wetratiodiff(BRNS.wetratiodiff.Seas[[s]], brks, col, zrange, BRNS.missing)
  mtext(season.names[[s]], side = 3,line=-2, adj = 0.1)
}
colkey(col=col, side=1, clim=range(brks), width=8, at=brks, labels = paste(round(zbrks,1)))
mtext("(a)", side = 3, line = -3, adj = -0.12, cex = 1.7, outer = T)
mtext("(b)", side = 3, line = -25, adj = -0.12, cex = 1.7, outer = T)

Quantile Ratio and Normalised Difference based scaling calculations

my.qqplot function This function is needed inside the calc.quantsignif function

# My qqplot function
my.qqplot <- function(x, y, minlen, plot.it = TRUE, xlab = deparse(substitute(x)),
                      ylab = deparse(substitute(y)), ...)
{
  sx <- sort(x)
  sy <- sort(y)
  lenx <- length(sx)
  leny <- length(sy)
  sx <- approx(1L:lenx, sx, n = minlen)$y
  sy <- approx(1L:leny, sy, n = minlen)$y
  if (plot.it)
    plot(sx, sy, xlab = xlab, ylab = ylab, ...)
  invisible(list(x = sx, y = sy))
}

calc.quantsignif function This is the main function that: 1. calculates the Normalised difference quantile ratios for each quantile at each grid cell 2. Does the bootstrapping to also give the confidence intervals for each ND statistic Normalised difference is defined by:
$$ ND = \frac{QR_2 - QR_1}{mean(QR_1, QR_2)} $$
The confidence intervals are given by the 5th and 95th percentiles of the bootstrapped distribution.

Bootstrapping is done by resampling in 2 yearly chunks a 1000 times to create 1000 estimates of 301 quantiles.

calc.quantSignif <- function(WET_1, WET_2, numQuantiles, dates) {
  dims <- dim(WET_1)
  if (!identical(dims, dim(WET_2))) warning("WET_1 and WET_2 dimensions different")
  dayYear <- year(dates)
  years <- unique(dayYear)
  yearpairs <- c(1:(length(years)/2))
  quantRatios <- array(dim=c(dims[1],dims[2],numQuantiles))
  confints.QR <- array(dim=c(dims[1],dims[2],numQuantiles,2)) # 2 for 0.05 and 0.95 confints
  #first calculate the quantRatios using the my.qqplot function
  for (i in 1:dims[1]){
    for (j in 1:dims[2]){
      if (j == 1) {writeLines(paste("Starting Longitude column: ", i))}
      if (sum(!is.na(WET_1[i,j,])) >= numQuantiles && sum(!is.na(WET_2[i,j,])) >= numQuantiles) {
        qx <- my.qqplot(WET_1[i,j,], WET_2[i,j,], plot.it = F, minlen = numQuantiles)
        quantRatios[i,j,] <- ifelse(!(qx$y - qx$x), 0, (qx$y - qx$x)/(apply(rbind(qx$x, qx$y), 2, mean, na.rm = T)))
        # Now calculate the confints
        #ptm1 <- proc.time()
        NullDist <- array(dim = c(1000,numQuantiles))
        for (b in 1:1000) {
          enoughSamples <- FALSE
          tries <- 1
          while(!enoughSamples & tries <= 10) {
            yearPairSample <- sample(yearpairs, replace = T)
            yearSample <- years[c(rbind(2*yearPairSample-1,2*yearPairSample))]
            daySample1 <- unlist(sapply(yearSample[1:(length(years)/2)], function(YSVector, dYVector) {return(which(dYVector==YSVector))}, dYVector=dayYear))
            daySample2 <- unlist(sapply(yearSample[((length(years)/2)+1):length(years)], function(YSVector, dYVector) {return(which(dYVector==YSVector))}, dYVector=dayYear))
            valSample1 <- c(WET_1[i,j,],WET_2[i,j,])[daySample1]
            valSample2 <- c(WET_1[i,j,],WET_2[i,j,])[daySample2]
            if (sum(!is.na(valSample1)) >= numQuantiles && sum(!is.na(valSample2)) >= numQuantiles) enoughSamples = TRUE
            tries <- tries + 1
          }
          vxy <- my.qqplot(valSample1, valSample2, plot.it = F, minlen = numQuantiles)
          NullDist[b,] <- ifelse(!(vxy$y - vxy$x), 0, (vxy$y - vxy$x)/(apply(rbind(vxy$x, vxy$y), 2, mean, na.rm = T)))
        }
        #ptm2 <- proc.time() # 11.725 seconds
        confints.QR[i,j,,] <- t(apply(NullDist, 2, function(d) quantile(d, probs=c(0.05,0.95))))
      } else {
        if (sum(!is.na(WET_1[i,j,])) > 0 || sum(!is.na(WET_2[i,j,])) > 0) {
          warning(paste0("Not enough non NA values to calculate ", numQuantiles, " quantiles: skipping [", i,",", j,"]th grid cell"))
        }
      }
    }
  }
  rm(dims,dayYear,years,yearpairs,qx,NullDist,yearPairSample,yearSample,daySample1,daySample2,valSample1,valSample2,vxy); gc()
  invisible(list(quantRatios=quantRatios,confints.QR=confints.QR))
}

Now we evaluate the above function for AWAP, BOA and NNI ** AWAP ** This is an overnight process.

AWAP.quantSignif <- calc.quantSignif(WET_AWAP_1, WET_AWAP_2, numQuantiles = 301, dates)

** BOA **

BRNS.quantSignif <- calc.quantSignif(WET_BRNS_1, WET_BRNS_2, 301, dates)

** NNI **

NATN.quantSignif <- calc.quantSignif(WET_NATN_1, WET_NATN_2, 301, dates)

Load pre-run results

load("RData/AWAP_FDD_quantile_analysis_AWAPquantSignig.RData")
load("RData/AWAP_FDD_quantile_analysis_BRNSquantSignif.RData")
load("RData/AWAP_FDD_quantile_analysis_NATNquantSignif.RData")

Calculate AWAP BOA and NNI gt1lt1ratio AWAP

xvals <- seq(0,1,length.out = 301)
numQuantiles = 301
AWAP.QR.signif <- AWAP.quantSignif$quantRatios
for (i in 1:length(lon)) {
  for (j in 1:length(lat)) {
    for (q in 1:numQuantiles) {
      if (!is.na(AWAP.QR.signif[i,j,q])) {
        if (AWAP.QR.signif[i,j,q] <= AWAP.quantSignif$confints.QR[i,j,q,2] && AWAP.QR.signif[i,j,q] >= AWAP.quantSignif$confints.QR[i,j,q,1]) {
          AWAP.QR.signif[i,j,q] <- NA
        }
      }
    }
  }
}
AWAP.QR.signif.gt1lt1ratio <- apply(AWAP.QR.signif, 3, function(x) length(which(x>1))/length(which(x<1)))

BOA

xvals <- seq(0,1,length.out = 301)
numQuantiles = 301
BRNS.QR.signif <- BRNS.quantSignif$quantRatios
for (i in 1:length(lon)) {
  for (j in 1:length(lat)) {
    for (q in 1:numQuantiles) {
      if (!is.na(BRNS.QR.signif[i,j,q])) {
        if (BRNS.QR.signif[i,j,q] <= BRNS.quantSignif$confints.QR[i,j,q,2] && BRNS.QR.signif[i,j,q] >= BRNS.quantSignif$confints.QR[i,j,q,1]) {
          BRNS.QR.signif[i,j,q] <- NA
        }
      }
    }
  }
}
BRNS.QR.signif.gt1lt1ratio <- apply(BRNS.QR.signif, 3, function(x) length(which(x>1))/length(which(x<1)))

NNI

xvals <- seq(0,1,length.out = 301)
numQuantiles = 301
NATN.QR.signif <- NATN.quantSignif$quantRatios
for (i in 1:length(lon)) {
  for (j in 1:length(lat)) {
    for (q in 1:numQuantiles) {
      if (!is.na(NATN.QR.signif[i,j,q])) {
        if (NATN.QR.signif[i,j,q] <= NATN.quantSignif$confints.QR[i,j,q,2] && NATN.QR.signif[i,j,q] >= NATN.quantSignif$confints.QR[i,j,q,1]) {
          NATN.QR.signif[i,j,q] <- NA
        }
      }
    }
  }
}
NATN.QR.signif.gt1lt1ratio <- apply(NATN.QR.signif, 3, function(x) length(which(x>1))/length(which(x<1)))

Calculate Normalised difference based scaling ND is calculated locally at each grid cell and then divided locally by a mean temperacture change and then averaged spatially.

First calculate the mean temperature change.

# Import temperature time series
path.TS <- "/srv/ccrc/data11/z3289452/Temp_timeseries/"
Temp.timeseries <- read.table(paste0(path.TS, "20170324.txt"), header = F)$V2
years <- seq(1910,2016)

# Split time series into two periods
temp.period1 <- Temp.timeseries[which(years >= 1958 & years <= 1985)] 
years.period1 <- years[which(years >= 1958 & years <= 1985)]

temp.period2 <- Temp.timeseries[which(years >= 1986 & years <= 2013)] 
years.period2 <- years[which(years >= 1986 & years <= 2013)]

# Find mean temperature change between two periods
mean.temp.change <- mean(temp.period2) - mean(temp.period1) # 0.42 K

Now do the calculation for AWAP, BOA and NNI. Run the initial data import code chunks. AWAP

probs = 301
# calculate quantiles locally first
qx <- apply(WET_AWAP_1, c(1,2), quantile, probs = probs, na.rm = T)
qy <- apply(WET_AWAP_2, c(1,2), quantile, probs = probs, na.rm = T)
# 
ND.data <- mapply(function(x,y) {
  return((y - x)/mean(c(x,y)))
}, qx, qy)
WET_AWAP.ND.mapply <- array(data = ND.data, dim = dim(qx))
ND.scaling <- WET_AWAP.ND.mapply/mean.temp.change
WET_AWAP.ND.scaling.spatavg <- apply(ND.scaling, 1, mean, na.rm = T)
save(WET_AWAP.ND.scaling.spatavg, file = "RData/AWAP_FDD_quantile_analysis_WETAWAPNDscalingSpatavg.RData")

BOA

probs = 301
# calculate quantiles locally first
qx <- apply(WET_BRNS_1, c(1,2), quantile, probs = probs, na.rm = T)
qy <- apply(WET_BRNS_2, c(1,2), quantile, probs = probs, na.rm = T)
# 
ND.data <- mapply(function(x,y) {
  return((y - x)/mean(c(x,y)))
}, qx, qy)
WET_BRNS.ND.mapply <- array(data = ND.data, dim = dim(qx))
ND.scaling <- WET_BRNS.ND.mapply/mean.temp.change
WET_BRNS.ND.scaling.spatavg <- apply(ND.scaling, 1, mean, na.rm = T)
save(WET_BRNS.ND.scaling.spatavg, file = "RData/AWAP_FDD_quantile_analysis_WETBRNSNDscalingSpatavg.RData")

NNI

probs = 301
# calculate quantiles locally first
qx <- apply(WET_NATN_1, c(1,2), quantile, probs = probs, na.rm = T)
qy <- apply(WET_NATN_2, c(1,2), quantile, probs = probs, na.rm = T)
# 
ND.data <- mapply(function(x,y) {
  return((y - x)/mean(c(x,y)))
}, qx, qy)
WET_NATN.ND.mapply <- array(data = ND.data, dim = dim(qx))
ND.scaling <- WET_NATN.ND.mapply/mean.temp.change
WET_NATN.ND.scaling.spatavg <- apply(ND.scaling, 1, mean, na.rm = T)
save(WET_NATN.ND.scaling.spatavg, file = "RData/AWAP_FDD_quantile_analysis_WETNATNNDscalingSpatavg.RData")

load pre-run results for scaling

numQuantiles <- 301
quantiles <- seq(0,1,length.out = numQuantiles)
load("RData/AWAP_FDD_quantile_analysis_WETAWAPNDscalingSpatavg.RData")
load("RData/AWAP_FDD_quantile_analysis_WETBRNSNDscalingSpatavg.RData")
load("RData/AWAP_FDD_quantile_analysis_WETNATNNDscalingSpatavg.RData")

Now we plot both the gt1lt1ratio and scaling plots on the same plot Besides the data directly above you also need to load data from section "Calculate AWAP BOA and NNI gt1lt1ratio"

zrange <- range(AWAP.QR.signif.gt1lt1ratio, BRNS.QR.signif.gt1lt1ratio, NATN.QR.signif.gt1lt1ratio)
layout(mat = 1)
par(mar=c(4.1,5.1,1.5,4.1))
plot(xvals,AWAP.QR.signif.gt1lt1ratio, ylim = zrange, pch=19, cex=0.5,xlab="Quantiles", 
     ylab="Dots: Ratio of number of grids 
     with QR>1 to number of grids with QR<1 ")
points(xvals,BRNS.QR.signif.gt1lt1ratio, col="red", pch=19, cex=0.5)
points(xvals,NATN.QR.signif.gt1lt1ratio, col="blue", pch=19, cex=0.5)
abline(h=1)
par(new = T)
plot(quantiles, WET_AWAP.ND.scaling.spatavg, type = "n", axes = F,
     xlab = NA, ylab = NA)
axis(side = 4)
mtext(side = 4, line = 3, "Lines: Mean temperature scaling (1/K)")
lines(quantiles, WET_AWAP.ND.scaling.spatavg, lwd = 2)
lines(quantiles, WET_BRNS.ND.scaling.spatavg, lwd = 2, col = "red")
lines(quantiles, WET_NATN.ND.scaling.spatavg, lwd = 2, col = "blue")
abline(h = 0.07, lty = 3)


legend("topleft",c('AWAP','BRNS','NATN'),pch = 19,col=c("black","red","blue"))

Supplementary Materials

Map of difference in total rainfall between periods

Import or load NNI wet day data

load("RData/AWAP_FDD_quantile_analysis_InitNCImport.RData")
WET_AWAP_2 <- WET_AWAP_2[,,2:dim(WET_AWAP_2)[3]] # correcting for an indexing error on lines 145, 154, 162
WET_BRNS_2 <- WET_BRNS_2[,,2:dim(WET_BRNS_2)[3]]
WET_NATN_2 <- WET_NATN_2[,,2:dim(WET_NATN_2)[3]]

and run the "masks" code chunk.
Calculate the difference in total rainfall between period 1 and period 2 (2 - 1).

WET_NATN_1.total <- apply(WET_NATN_1, c(1,2), FUN = sum, na.rm = T)
WET_NATN_2.total <- apply(WET_NATN_2, c(1,2), FUN = sum, na.rm = T)
WET_NATN.total.diff <- WET_NATN_2.total - WET_NATN_1.total
# Apply mask
WET_NATN.total.diff[!is.na(AWAP.missing)] <- NA

Do the plotting.

# Plot
# determine brks and col for plotting
zrange <- range(WET_NATN.total.diff, na.rm = T)
brks <- seq(-8000,12000,1000)
zbrks <- brks; zbrks[1] <- zrange[1]; zbrks[length(zbrks)] <- zrange[2]
col <- create.col.list(brks, 0)

layout(mat = c(1,2), heights = c(0.85, 0.15))
mar = par()$mar; mar[1] <- 4; mar[3] <- 1.1
par(mar = mar, oma = c(1.1, 0, 0, 0))
image(lon, lat, WET_NATN.total.diff, breaks = zbrks, col = col, xlab = "Longitude", ylab = "Latitude")
image(lon, lat, AWAP.missing, col = "gray", add = T)
map("world", "Australia", add = T)
colkey(col = col, side = 1, clim = range(brks), width = 5, at = brks, labels = paste(round(zbrks,1)))
mtext("Difference in total rainfall (period 2 - period 1) (mm)", side = 1, line = 4)

REVISIONS

The biggest common criticism from both reviewers was the use of WET day quantiles as opposed to all day quantiles. On the recommendation of one of the reviewers (Christoph Schar) we decided to redo the quantile change calculations. We implemented the following changes:

  1. Replaced the quantile ratio statistic by the relative change (or normalised difference) statistic.
  2. Calculate the relative quantile change using
    1. the maximum number of non zero days (> 0 mm) (f) between the two periods and then calculating the change between the top (f) data values.
    2. the minimum number of wet days (>= 1mm) (f) between the two periods and then calculating the change between the top (f) data values.

Below we present the code to redo the quantile change calculation and panel plot, and the scaling calculation and panel plot.

# Let's mask the two arrays
mask <- function(x, mask) {
  stopifnot(identical(dim(x), dim(mask)))
  x[mask] <- NA
  return(x)
}

calc.quantiles <- function(ts1, ts2, numQuantiles = NA) {
  stopifnot(length(ts1) == length(ts2))
  ts1.sorted <- sort(ts1)
  len1 <- length(ts1.sorted)
  ts2.sorted <- sort(ts2)
  len2 <- length(ts2.sorted)
  f1 <- length(which(ts1 > 0))
  f2 <- length(which(ts2 > 0))
  f <- max(f1, f2)
  if (is.na(numQuantiles)) {
    warning("Number of quantiles not supplied. Calculating maximum possible")
    numQuantiles <- f
  }
  ts1.include <- ts1.sorted[(len1 - f + 1):len1]
  ts2.include <- ts2.sorted[(len2 - f + 1):len2]
  x <- approx(1L:f, ts1.include, n = numQuantiles)$y
  y <- approx(1L:f, ts2.include, n = numQuantiles)$y
  return(list(x = x, y = y))
}

calc.wetquantiles <- function(ts1, ts2, numQuantiles = NA) {
  stopifnot(length(ts1) == length(ts2))
  ts1.sorted <- sort(ts1)
  len1 <- length(ts1.sorted)
  ts2.sorted <- sort(ts2)
  len2 <- length(ts2.sorted)
  f1 <- length(which(ts1 >= 1))
  f2 <- length(which(ts2 >= 1))
  f <- min(f1, f2)
  if (is.na(numQuantiles)) {
    warning("Number of quantiles not supplied. Calculating maximum possible")
    numQuantiles <- f
  }
  ts1.include <- ts1.sorted[(len1 - f + 1):len1]
  ts2.include <- ts2.sorted[(len2 - f + 1):len2]
  x <- approx(1L:f, ts1.include, n = numQuantiles)$y
  y <- approx(1L:f, ts2.include, n = numQuantiles)$y
  return(list(x = x, y = y))
}

rel.change <- function(x, y) {
  stopifnot(length(x) == length(y))
  return(mapply(function(a, b) {
    return((b - a)*2/(a + b))
  }, x, y))
}

create.col.list <- function(brks, null.point) {
  require(graphics)
  RdOl <- colorRampPalette(c("darkred","orchid","yellow"))
  OlBl <- colorRampPalette(c("lightgreen", "slateblue", "blue4"))
  num.red <- which.min(abs(brks - null.point))-1
  num.blue <- (length(brks)-1)-num.red
  if(num.red > num.blue){
    return(c(RdOl(num.red),OlBl(num.red)[1:num.blue]))
  } else {
    return(c(RdOl(num.blue)[(num.blue-num.red+1):num.blue],OlBl(num.blue)))
  }
}

calc.quantChange <- function(quantSignif) {
  dim <- dim(quantSignif[[1]])
  stopifnot(length(dim) == 3)
  quantChange.signif <- quantSignif[[1]]
  for (i in 1:dim[1]) {
    for (j in 1:dim[2]) {
      for (q in 1:dim[3]) {
        if (!is.na(quantChange.signif[i,j,q])) {
          if ((quantChange.signif[i,j,q] <= quantSignif[[2]][i,j,q,2] && quantChange.signif[i,j,q] >= quantSignif[[2]][i,j,q,1]) || any(is.na(quantSignif[[2]][i,j,q,]))) {
            quantChange.signif[i,j,q] <- NA
          }
        }
      }
    }
  }
  QCsignif.gt1lt1ratio <- apply(quantChange.signif, 3, function(x) length(which(x > 0))/length(which(x < 0)))
  return(list(quantChange.signif = quantChange.signif, QCsignif.gt1lt1ratio = QCsignif.gt1lt1ratio))
}

plot.quantChangePanel <- function(quantChange, lon, lat, mask) {
  dim = dim(quantChange)
  xvals <- seq(0,1,length.out = dim[3])
  zrange <- range(quantChange, na.rm = T)
  brks <- seq(-0.6,1,0.2)
  zbrks <- brks; zbrks[1] <- min(brks[1],zrange[1]); zbrks[length(brks)] <- max(brks[length(brks)],zrange[2]) # since 0.3 < zrange[1]
  col <- create.col.list(brks,0)
  col[1] <- "indianred1"
  quant.index <- c(seq(31,by = 30,length.out = 9),seq(286,301,by = 3)) # c(seq(1,by = 15,length.out = 20),295,seq(298,301))
  layout(mat = matrix(c(seq(1,15),rep(16,5)), nrow = 4,ncol = 5,byrow = T), heights = c(1,1,1,0.5))
  par(mar = c(0,0,0,0),oma = c(1,1.5,1,1.5))
  for (i in 1:15) {
    image(lon,lat,quantChange[,,quant.index[i]], 
          breaks = zbrks, xaxt = 'n', yaxt = 'n', xlab = "",col = col)
    image(lon, lat, mask, col="gray", add = T)
    map("world", "Australia", add = T)
    if (i < 21) {
      mtext(sprintf("%.2f",xvals[quant.index[i]]),side = 3,line = -1.8,adj = 0.05, cex = 1.3)
    } else {
      mtext(sprintf("%.3f",xvals[quant.index[i]]),side = 3,line = -1.8,adj = 0.05, cex = 1.3)
    }
  }
  colkey(col = col, side = 3, clim = range(brks), width = 5, at = brks, labels = paste(round(zbrks,1)), cex.axis = 1.6)
}

# Now the scaling
# calculate quantiles locally first
calc.scaling <- function(PRCP1, PRCP2, numQuantiles, WET = FALSE, mean.temp.change) {
  if (missing(mean.temp.change)) stop("Must supply the mean temperature change for scaling")
  stopifnot(identical(dim(PRCP1), dim(PRCP2)))
  dim = dim(PRCP1)
  stopifnot(length(dim) == 3)
  quantiles1 <- array(dim = c(dim[1], dim[2], numQuantiles))
  quantiles2 <- array(dim = c(dim[1], dim[2], numQuantiles))
  for (i in 1:dim[1]) {
    writeLines(paste("Processing column", i))  
    for (j in 1:dim[2]) {
      if (!sum(is.na(PRCP1[i,j,])) & !sum(is.na(PRCP2[i,j,]))) {
        if (WET) {
          quantiles <- calc.wetquantiles(PRCP1[i,j,],PRCP2[i,j,], numQuantiles = numQuantiles)
        } else {
          quantiles <- calc.quantiles(PRCP1[i,j,],PRCP2[i,j,], numQuantiles = numQuantiles)
        }
        quantiles1[i,j,] <- quantiles$x
        quantiles2[i,j,] <- quantiles$y
      }
    }
  }
  relChange <- rel.change(quantiles1, quantiles2)
  relChange.mapply <- array(data = relChange, dim = dim(quantiles1))
  relChange.scaling <- relChange.mapply/mean.temp.change
  relChange.scaling.spatavg <- apply(relChange.scaling, 3, mean, na.rm = T)
  return(relChange.scaling.spatavg)
}

# Plot
plot.scaling <- function(scaling) {
  par(mar = c(5,5,4,2)+0.1)
  plot(xvals, scaling, xlab = "Quantiles", ylab = "Mean scaling (1/K)", 
       type = "n", cex.lab = 1.75, cex.axis = 1.75)
  lines(xvals, scaling, lwd = 2)
  abline(h = 0.07, lty = 3)
}

We calculate the relative quantile change and confidence intervals in the following script using bootstrapping.

# R Script to calculate relative change in quantiles
# and confidence intervals for non zero values and 
# for all wet day values

library(lubridate)
library(RNetCDF)

sink("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/calcRelQuantChangewConfints_AWAP.R",
     append = F, split = T)

# Data import
# Before importing data, both AWAP and GPCC_FDD were regridded onto the same 1X1 degree grid using cdo remapcon2
dates <- seq(as.Date("1958-01-01"), as.Date("2013-12-31"),1)
lon <- seq(112, by=0.5, length.out = 88)
lat <- rev(seq(-10, by=-0.5, length.out = 68))

# Import Natn data
nc.dir <- "/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/NatgridBarnes/"
nc.name <- "Prcp_daily_natgrid_0.5x0.5deg_1958-2013_land.nc"
nc <- open.nc(paste0(nc.dir,nc.name))
PRCP_NATN <- var.get.nc(nc, "Prcp")
close.nc(nc)
PRCP_NATN_1 <- PRCP_NATN[,,1:10227] # 10227->1985-12-31 10592->1986-12-31
PRCP_NATN_2 <- PRCP_NATN[,,10228:dim(PRCP_NATN)[3]] # 20454->2013-12-31
rm(PRCP_NATN); gc()

# Load AWAP BRNS and NATN masks
load("RData/AWAP_BRNS_NATN_mask.RData")
# Get rid of Giles
AWAP_mask[31:36,36:40] <- 0 
BRNS_mask[31:36,36:40] <- 0
# First calculate AWAP.missing mask
AWAP.missing <- AWAP_mask
AWAP.missing[is.na(AWAP.missing)] <- 0
AWAP.missing[which(AWAP.missing == 1)] <- NA
AWAP.missing <- !AWAP.missing

ptm1 <- proc.time()
NATN.quantSignif <- calc.quantSignif(PRCP1 = PRCP_NATN_1, PRCP2 = PRCP_NATN_2, 301, dates)
ptm2 <- proc.time()
writeLines(paste("Time taken for NATN.quantsignif calculation:", ptm2 - ptm1))

ptm1 <- proc.time()
NATN.wetquantSignif <- calc.wetquantSignif(PRCP1 = PRCP_NATN_1, PRCP2 = PRCP_NATN_2, 301, dates)
ptm2 <- proc.time()
writeLines(paste("Time taken for NATN.wetquantsignif calculation:", ptm2 - ptm1))

save(NATN.quantSignif, NATN.wetquantSignif, file = "/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/calcRelQuantChangewConfints_AWAP.RData")

Run the following to set things up.

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

#  setwd 
dir = "/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/"
setwd(dir)

# Data import
# Before importing data, both AWAP and GPCC_FDD were regridded onto the same 1X1 degree grid using cdo remapcon2
lon <- seq(112, by=0.5, length.out = 88)
lat <- rev(seq(-10, by=-0.5, length.out = 68))

#Create missing masks for plotting
# Load AWAP BRNS and NATN masks
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/AWAP_BRNS_NATN_mask.RData")
# Get rid of Giles
AWAP_mask[31:36,36:40] <- 0 
BRNS_mask[31:36,36:40] <- 0 

#LSM
# Import and land sea mask from a netcdf file
nc.dir <- "/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/"
LSM.nc <- open.nc(paste0(nc.dir, "AWAP_Land-Sea-Mask_0.5deg_latflipped.nc"))
LSM <- var.get.nc(LSM.nc, "LSM")
rm(nc.dir, LSM.nc); gc()

#AWAP
AWAP.missing <- AWAP_mask
AWAP.missing[is.na(AWAP.missing)] <- 0
AWAP.missing[!LSM] <- NA
AWAP.missing[which(AWAP.missing == 1)] <- NA
AWAP.missing <- !AWAP.missing
#BRNS
BRNS.missing <- BRNS_mask
BRNS.missing[is.na(BRNS.missing)] <- 0
BRNS.missing[!LSM] <- NA
BRNS.missing[which(BRNS.missing == 1)] <- NA
BRNS.missing <- !BRNS.missing

First NNI

#load the data from file that we saved in the previous code chunk
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/calcRelQuantChangewConfints_NATN.RData")

NATN.quantChange <- calc.quantChange(NATN.quantSignif)
NATN.wetquantChange <- calc.quantChange(NATN.wetquantSignif)
plot.quantChangePanel(NATN.quantChange[[1]], lon, lat, AWAP.missing)
plot.quantChangePanel(NATN.wetquantChange[[1]], lon, lat, AWAP.missing)
# Now all grids, not just significant ones
plot.quantChangePanel(NATN.quantSignif$relQuantChange, lon, lat, AWAP.missing)
plot.quantChangePanel(NATN.wetquantSignif$relQuantChange, lon, lat, AWAP.missing)

Now AWAP

#load the data from file that we saved in the previous code chunk
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/calcRelQuantChangewConfints_AWAP.RData")

AWAP.quantChange <- calc.quantChange(AWAP.quantSignif)
AWAP.wetquantChange <- calc.quantChange(AWAP.wetquantSignif)
plot.quantChangePanel(AWAP.quantChange[[1]], lon, lat, AWAP.missing)
plot.quantChangePanel(AWAP.wetquantChange[[1]], lon, lat, AWAP.missing)
# Now all grids, not just significant ones
plot.quantChangePanel(AWAP.quantSignif$relQuantChange, lon, lat, AWAP.missing)
plot.quantChangePanel(AWAP.wetquantSignif$relQuantChange, lon, lat, AWAP.missing)

Now BRNS

#load the data from file that we saved in the previous code chunk
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/calcRelQuantChangewConfints_BRNS.RData")

BRNS.quantChange <- calc.quantChange(BRNS.quantSignif)
BRNS.wetquantChange <- calc.quantChange(BRNS.wetquantSignif)
plot.quantChangePanel(BRNS.quantChange[[1]], lon, lat, BRNS.missing)
plot.quantChangePanel(BRNS.wetquantChange[[1]], lon, lat, BRNS.missing)
# Now all grids, not just significant ones
plot.quantChangePanel(BRNS.quantSignif$relQuantChange, lon, lat, BRNS.missing)
plot.quantChangePanel(BRNS.wetquantSignif$relQuantChange, lon, lat, BRNS.missing)

Now the scaling plots:

First the data import and calculation:

library(lubridate)
library(RNetCDF)

#  setwd 
dir = "/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/"
setwd(dir)

# Data import
# Before importing data, both AWAP and GPCC_FDD were regridded onto the same 1X1 degree grid using cdo remapcon2
lon <- seq(112, by=0.5, length.out = 88)
lat <- rev(seq(-10, by=-0.5, length.out = 68))
dates <- seq(as.Date("1958-01-01"), as.Date("2013-12-31"),1)

# Load AWAP BRNS and NATN masks
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/AWAP_BRNS_NATN_mask.RData")
# Get rid of Giles
AWAP_mask[31:36,36:40] <- 0 
BRNS_mask[31:36,36:40] <- 0

# First calculate AWAP.missing mask
AWAP.missing <- AWAP_mask
AWAP.missing[is.na(AWAP.missing)] <- 0
AWAP.missing[which(AWAP.missing == 1)] <- NA
AWAP.missing <- !AWAP.missing
# First calculate BRNS.missing mask
BRNS.missing <- BRNS_mask
BRNS.missing[is.na(BRNS.missing)] <- 0
BRNS.missing[which(BRNS.missing == 1)] <- NA
BRNS.missing <- !BRNS.missing

# Import AWAP data
nc.dir <- "/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/AWAP/" 
nc.name <- "AWAP_PRCP_1900-2015_0.5deg_remapcon2_land_noneg.nc"
nc <- open.nc(paste0(nc.dir, nc.name))
PRCP_AWAP <- var.get.nc(nc, "pre", start=c(1,1,21185), count=c(88,68,length(dates)))
close.nc(nc); rm(nc); gc()
PRCP_AWAP <- PRCP_AWAP[,dim(PRCP_AWAP)[2]:1,]
#Limit analysis to 2013 because 2014-15 data for Natgrid and Barnes not available
PRCP_AWAP_1 <- PRCP_AWAP[,,1:10227] # 10227->1985-12-31 10592->1986-12-31
PRCP_AWAP_2 <- PRCP_AWAP[,,10228:dim(PRCP_AWAP)[3]] # 20454->2013-12-31
# Mask the data
for (t in 1:10227) {
  PRCP_AWAP_1[,,t] <- mask(PRCP_AWAP_1[,,t], AWAP.missing)
  PRCP_AWAP_2[,,t] <- mask(PRCP_AWAP_2[,,t], AWAP.missing)
}

# Calculate Scaling
# You need mean temperature change for this
# Run code chunk above!
AWAP.scaling <- calc.scaling(PRCP_AWAP_1, PRCP_AWAP_2, numQuantiles = 301, WET = F, mean.temp.change = mean.temp.change)
AWAP.wetscaling <- calc.scaling(PRCP_AWAP_1, PRCP_AWAP_2, numQuantiles = 301, WET = T, mean.temp.change = mean.temp.change)

#Cleanup
rm(PRCP_AWAP, PRCP_AWAP_1, PRCP_AWAP_2); gc()

# Import Natn data
nc.dir <- "/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/NatgridBarnes/"
nc.name <- "Prcp_daily_natgrid_0.5x0.5deg_1958-2013_land.nc"
nc <- open.nc(paste0(nc.dir,nc.name))
PRCP_NATN <- var.get.nc(nc, "Prcp")
close.nc(nc)
PRCP_NATN_1 <- PRCP_NATN[,,1:10227] # 10227->1985-12-31 10592->1986-12-31
PRCP_NATN_2 <- PRCP_NATN[,,10228:dim(PRCP_NATN)[3]] # 20454->2013-12-31
# Mask the data
for (t in 1:10227) {
  PRCP_NATN_1[,,t] <- mask(PRCP_NATN_1[,,t], AWAP.missing)
  PRCP_NATN_2[,,t] <- mask(PRCP_NATN_2[,,t], AWAP.missing)
}

# Calculate Scaling
# You need mean temperature change for this
# Run code chunk above!
NATN.scaling <- calc.scaling(PRCP_NATN_1, PRCP_NATN_2, numQuantiles = 301, WET = F, mean.temp.change = mean.temp.change)
NATN.wetscaling <- calc.scaling(PRCP_NATN_1, PRCP_NATN_2, numQuantiles = 301, WET = T, mean.temp.change = mean.temp.change)

#Cleanup
rm(PRCP_NATN, PRCP_NATN_1, PRCP_NATN_2); gc()

# Import Brns data
nc.dir <- "/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/NatgridBarnes/"
nc.name <- "Prcp_daily_barnesgrid_0.5x0.5deg_1958-2013_Req1_land.nc"
nc <- open.nc(paste0(nc.dir,nc.name))
PRCP_BRNS <- var.get.nc(nc, "Prcp")
close.nc(nc)
# Natgrid and Barnes wet prcp for two periods
PRCP_BRNS_1 <- PRCP_BRNS[,,1:10227] # 10227->1985-12-31 10592->1986-12-31
PRCP_BRNS_2 <- PRCP_BRNS[,,10228:dim(PRCP_BRNS)[3]] # 20454->2013-12-31
# Mask the data
for (t in 1:10227) {
  PRCP_BRNS_1[,,t] <- mask(PRCP_BRNS_1[,,t], AWAP.missing)
  PRCP_BRNS_2[,,t] <- mask(PRCP_BRNS_2[,,t], AWAP.missing)
}

# Calculate Scaling
# You need mean temperature change for this
# Run code chunk above!
BRNS.scaling <- calc.scaling(PRCP_BRNS_1, PRCP_BRNS_2, numQuantiles = 301, WET = F, mean.temp.change = mean.temp.change)
BRNS.wetscaling <- calc.scaling(PRCP_BRNS_1, PRCP_BRNS_2, numQuantiles = 301, WET = T, mean.temp.change = mean.temp.change)

#Cleanup
rm(PRCP_BRNS, PRCP_BRNS_1, PRCP_BRNS_2); gc()

Now we do the plotting:

plot.scaling <- function(AWAP.scaling, BRNS.scaling, NATN.scaling) {
  xvals <- seq(0, 1, length.out = length(AWAP.scaling))
  par(mar = c(5,5,4,2) + 0.1)
  plot(xvals, AWAP.scaling, xlab = "Quantiles", ylab = "Mean scaling (1/K)", 
       type = "n", cex.lab = 1.75, cex.axis = 1.75,
       ylim = range(AWAP.scaling, BRNS.scaling, NATN.scaling, na.rm = T))
  lines(xvals, AWAP.scaling, lwd = 2, col = "black")
  lines(xvals, BRNS.scaling, lwd = 2, col = "red")
  lines(xvals, NATN.scaling, lwd = 2, col = "blue")
  abline(h = 0.07, lty = 3)
}

plot.scaling(AWAP.scaling, BRNS.scaling, NATN.scaling)
plot.scaling(AWAP.wetscaling, BRNS.wetscaling, NATN.wetscaling)

For the lower quantiles, often the other quantile is zero which results in a relative change of either -2 (-x/(x/2)) or +2. Furthermore, NNI exhibits some huge changes of around -25 at the lower end and +82 on the upper end. As such we have a look at the range of values for the significant quantile changes in wet days. For NNI these are [-0.9, 1.2] and for BOA and AWAP these are [-1, 1.3]. So to see the scaling curve without these "outliers", we restrict the quantile changes to [-1, 1.3] as it is inclusive of the NNI range. After doing this we get the following curve:

Meeting with Christoph Schar

  1. Stop mentioning that daily rainfall has intensified. Rephrase it as when it rains, the magnitudes are more intense. Or in other words, intensification of conditional rainfall or something like that.
  2. Figure 2 should also show maps based on all day distributions (based on his method) and text should explain the differences. He argued that it is confusing to a hydrologists who could try and draw return periods based on my wet day distributions.
  3. Scaling stuff should be removed because scaling cannot be calculated using in-consistent all-day quantiles. If we would like to show the variability between datasets, we can show using the ratio of positive vs negative changes. This is actually really interesting because now we can show that all-day quantiles are really sensitive to frequency changes and the interpolation method (this is apparently a known result and he will send me some papers), and wet-day quantiles are insensitive to frequency changes and interpolation method but very sensitive to the underlying station network! I think that is pretty profound.

Regarding point number 3, maybe it is worthwhile chatting with a hydrologist and a statistician regarding the validity of scaling calculations.

Things to do

  1. Change title of paper
  2. Check wording of paper, especially in parts where I say rainfall has intensified. Instead say magnitude of wet day rainfall has increased.
  3. Redo figure 2 as follows:
  4. Redo Figure 3 as follows:
  5. Find papers that say that scaling cannot be calculated with wet-day quantiles.
  6. Find paper that show that all day quantiles are sensitive to frequency changes and interpolation method and wet day quantiles are sensitive to underlying station distribution.

Wednesday, 10th Jan 2018

We will now recreate maps from Figure 2 but with stippling

#load the data from file that we saved in the previous code chunk
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/calcRelQuantChangewConfints_NATN.RData")
# load true all al-day quantile calculation
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/NATN_alldayquantSignif.RData")

NATN.quantChange <- calc.quantChange(NATN.quantSignif)
NATN.wetquantChange <- calc.quantChange(NATN.wetquantSignif)
NATN.allquantChange <- calc.quantChange(NATN.alldayquantSignif)

plot.quantChangePanelwStippling <- function(quantChange, quantChangeSignif, lon, lat, mask) {
  dim = dim(quantChange)
  xvals <- seq(0,1,length.out = dim[3])
  zrange <- range(quantChange, na.rm = T)
  brks <- seq(-0.6,1,0.2)
  zbrks <- brks; zbrks[1] <- min(brks[1],zrange[1]); zbrks[length(brks)] <- max(brks[length(brks)],zrange[2]) # since 0.3 < zrange[1]
  col <- create.col.list(brks,0)
  col[1] <- "indianred1"
  quant.index <- c(seq(31,by = 30,length.out = 9),seq(286,301,by = 3)) # c(seq(1,by = 15,length.out = 20),295,seq(298,301))
  layout(mat = matrix(c(seq(1,15),rep(16,5)), nrow = 4,ncol = 5,byrow = T), heights = c(1,1,1,0.5))
  par(mar = c(0,0,0,0),oma = c(1,1.5,1,1.5))
  for (i in 1:15) {
    image(lon,lat,quantChange[,,quant.index[i]], 
          breaks = zbrks, xaxt = 'n', yaxt = 'n', xlab = "",col = col)
    stippling.data <- !is.na(quantChangeSignif[,,quant.index[i]])
    stippling.xy <- expand.grid(lon, lat)[stippling.data,]
    points(stippling.xy[,1], stippling.xy[,2], pch = 20, cex = 0.25, col = "gray")
    image(lon, lat, mask, col="gray", add = T)
    map("world", "Australia", add = T)
    if (i < 21) {
      mtext(sprintf("%.2f",xvals[quant.index[i]]),side = 3,line = -1.8,adj = 0.05, cex = 1.3)
    } else {
      mtext(sprintf("%.3f",xvals[quant.index[i]]),side = 3,line = -1.8,adj = 0.05, cex = 1.3)
    }
  }
  colkey(col = col, side = 3, clim = range(brks), width = 5, at = brks, labels = paste(round(zbrks,1)), cex.axis = 1.6)
}

plot.quantChangePanelwnonSigniffade <- function(quantChange, quantChangeSignif, lon, lat, mask, fade) {
  dim = dim(quantChange)
  xvals <- seq(0,1,length.out = dim[3])
  zrange <- range(quantChange, na.rm = T)
  brks <- seq(-0.12,0.2, length.out = 9)
  zbrks <- brks; zbrks[1] <- min(brks[1],zrange[1]); zbrks[length(brks)] <- max(brks[length(brks)],zrange[2]) # since 0.3 < zrange[1]
  col <- create.col.list(brks,0)
  col[1] <- "indianred1"
  fadedcol <- apply(col2rgb(col),2,FUN = function(x) {return(rgb(red = x[1]/255, green = x[2]/255, blue = x[3]/255, alpha = fade))})
  quant.index <- c(seq(31,by = 30,length.out = 9),seq(286,301,by = 3)) # c(seq(1,by = 15,length.out = 20),295,seq(298,301))
  layout(mat = matrix(c(seq(1,15),rep(16,5)), nrow = 4,ncol = 5,byrow = T), heights = c(1,1,1,0.5))
  par(mar = c(0,0,0,0),oma = c(1,1.5,1,1.5))
  for (i in 1:15) {
    image(lon,lat,quantChange[,,quant.index[i]], 
          breaks = zbrks, xaxt = 'n', yaxt = 'n', xlab = "",col = fadedcol)
    image(lon, lat, quantChangeSignif[,,quant.index[i]],
          breaks = zbrks,col = col, add = T)
    image(lon, lat, mask, col="gray", add = T)
    map("world", "Australia", add = T)
    if (i < 21) {
      mtext(sprintf("%.2f",xvals[quant.index[i]]),side = 3,line = -1.8,adj = 0.05, cex = 1.3)
    } else {
      mtext(sprintf("%.3f",xvals[quant.index[i]]),side = 3,line = -1.8,adj = 0.05, cex = 1.3)
    }
  }
  colkey(col = col, side = 3, clim = range(brks), width = 5, at = brks, labels = paste(round(zbrks,2)), cex.axis = 1.6)
}

# plot.quantChangePanel(NATN.quantChange[[1]], lon, lat, AWAP.missing)
# plot.quantChangePanel(NATN.wetquantChange[[1]], lon, lat, AWAP.missing)
# Now all grids, not just significant ones
# plot.quantChangePanelwStippling(NATN.quantSignif$relQuantChange, NATN.quantChange$quantChange.signif, lon, lat, AWAP.missing)
# plot.quantChangePanel(NATN.wetquantSignif$relQuantChange, lon, lat, AWAP.missing)
plot.quantChangePanelwnonSigniffade(NATN.alldayquantSignif$relQuantChange, NATN.allquantChange$quantChange.signif, lon, lat, AWAP.missing, fade = 0.3)
#AWAP
# load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/calcRelQuantChangewConfints_AWAP.RData")
# 
# AWAP.quantChange <- calc.quantChange(AWAP.quantSignif)
# plot.quantChangePanelwnonSigniffade(AWAP.quantSignif$relQuantChange, AWAP.quantChange$quantChange.signif, lon, lat, AWAP.missing, fade = 0.3)

# load true all al-day quantile calculation
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/AWAP_allquantSignif.RData")

AWAP.allquantChange <- calc.quantChange(AWAP.allquantSignif)
plot.quantChangePanelwnonSigniffade(AWAP.allquantSignif$relQuantChange, NATN.allquantChange$quantChange.signif, lon, lat, AWAP.missing, fade = 0.3)
# #BOA
# load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/calcRelQuantChangewConfints_BRNS.RData")
# 
# BRNS.quantChange <- calc.quantChange(BRNS.quantSignif)
# plot.quantChangePanelwnonSigniffade(BRNS.quantSignif$relQuantChange, BRNS.quantChange$quantChange.signif, lon, lat, BRNS.missing, fade = 0.3)

# load true all al-day quantile calculation
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/BRNS_allquantsignif.RData")

BRNS.allquantChange <- calc.quantChange(BRNS.allquantSignif)
plot.quantChangePanelwnonSigniffade(BRNS.allquantSignif$relQuantChange, NATN.allquantChange$quantChange.signif, lon, lat, BRNS.missing, fade = 0.3)

NNI AWAP BOA

# Create Figure 2a (wet day quantile change maps) but with a modified 2mm threshold
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/calcRelQuantChangewConfints_NNI_wetDays_2mmthresh.RData")

NATN.wet2quantchange <- calc.quantChange(NATN.quantSignif)
plot.quantChangePanelwnonSigniffade(NATN.quantSignif$relQuantChange, NATN.wet2quantchange$quantChange.signif, lon, lat, AWAP.missing, fade = 0.3)

Difference in number of grids with positive changes and negative changes for each map

Now we create plots of difference between the number of grid cells that show positive changes in quantiles and those that show negative changes in quantiles. We create two such plots, one using the wet day distribution and the other using all day distribution, with the samples where both periods are 0mm are discarded (max frequency).

First run the code chunk loading the lat, lon, dates and masks and the required packages.

calc.postonegchanges <- function(relquantchange) {
  return(apply(relquantchange, 3, FUN = function(a) {
    #2770 is the number of land grids for AUSTRALIA
    return((length(which(a > 0)) - length(which(a < 0)))/2770)
  }))
}
# create a combined AWAP and BRNS mask
mask <- AWAP.missing
mask[BRNS.missing] <- T
# wet days first
# load the data for wet day quantiles first
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/calcRelQuantChangewConfints_NNI_wetDays.RData")
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/calcRelQuantChangewConfints_AWAP_wetDays.RData")
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/calcRelQuantChangewConfints_BRNS_wetDays.RData")

# create an array of masks with the same dimension as NATN.QR.signif.wet
mask_array <- array(dim = dim(NATN.quantSignif$relQuantChange), data = mask)

numQuantiles = 301
NATN.QR.signif.wet <- NATN.quantSignif$relQuantChange
# apply mask
NATN.quantSignif$relQuantChange[mask_array] <- NA
NATN.QR.signif.wet[mask_array] <- NA
for (i in 1:length(lon)) {
  for (j in 1:length(lat)) {
    for (q in 1:numQuantiles) {
      if (!is.na(NATN.QR.signif.wet[i,j,q])) {
        if (NATN.QR.signif.wet[i,j,q] <= NATN.quantSignif$confints.QR[i,j,q,2] && NATN.QR.signif.wet[i,j,q] >= NATN.quantSignif$confints.QR[i,j,q,1]) {
          NATN.QR.signif.wet[i,j,q] <- NA
        }
      }
    }
  }
}
# AWAP
numQuantiles = 301
AWAP.QR.signif.wet <- AWAP.quantSignif$relQuantChange
# apply mask
AWAP.quantSignif$relQuantChange[mask_array] <- NA
AWAP.QR.signif.wet[mask_array] <- NA
for (i in 1:length(lon)) {
  for (j in 1:length(lat)) {
    for (q in 1:numQuantiles) {
      if (!is.na(AWAP.QR.signif.wet[i,j,q])) {
        if (AWAP.QR.signif.wet[i,j,q] <= AWAP.quantSignif$confints.QR[i,j,q,2] && AWAP.QR.signif.wet[i,j,q] >= AWAP.quantSignif$confints.QR[i,j,q,1]) {
          AWAP.QR.signif.wet[i,j,q] <- NA
        }
      }
    }
  }
}
# BRNS
numQuantiles = 301
BRNS.QR.signif.wet <- BRNS.quantSignif$relQuantChange
# apply mask
BRNS.quantSignif$relQuantChange[mask_array] <- NA
BRNS.QR.signif.wet[mask_array] <- NA
for (i in 1:length(lon)) {
  for (j in 1:length(lat)) {
    for (q in 1:numQuantiles) {
      if (!is.na(BRNS.QR.signif.wet[i,j,q])) {
        if (BRNS.QR.signif.wet[i,j,q] <= BRNS.quantSignif$confints.QR[i,j,q,2] && BRNS.QR.signif.wet[i,j,q] >= BRNS.quantSignif$confints.QR[i,j,q,1]) {
          BRNS.QR.signif.wet[i,j,q] <- NA
        }
      }
    }
  }
}

NNI.postonegchanges.wet <- calc.postonegchanges(NATN.QR.signif.wet)
AWAP.postonegchanges.wet <- calc.postonegchanges(AWAP.QR.signif.wet)
BOA.postonegchanges.wet <- calc.postonegchanges(BRNS.QR.signif.wet)
# Average changes in quantiles
NNI.AvgQuantChange.wet <- apply(NATN.quantSignif$relQuantChange, 3, mean, na.rm = T)
AWAP.AvgQuantChange.wet <- apply(AWAP.quantSignif$relQuantChange, 3, mean, na.rm = T)
BOA.AvgQuantChange.wet <- apply(BRNS.quantSignif$relQuantChange, 3, mean, na.rm = T)
# Now Hybrid All day quantiles
# load the data for all day quantiles first
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/calcRelQuantChangewConfints_NATN.RData")
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/calcRelQuantChangewConfints_AWAP.RData")
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/calcRelQuantChangewConfints_BRNS.RData")

numQuantiles = 301
NATN.QR.signif <- NATN.quantSignif$relQuantChange
for (i in 1:length(lon)) {
  for (j in 1:length(lat)) {
    for (q in 1:numQuantiles) {
      if (!is.na(NATN.QR.signif[i,j,q])) {
        if (NATN.QR.signif[i,j,q] <= NATN.quantSignif$confints.QR[i,j,q,2] && NATN.QR.signif[i,j,q] >= NATN.quantSignif$confints.QR[i,j,q,1]) {
          NATN.QR.signif[i,j,q] <- NA
        }
      }
    }
  }
}
# AWAP
numQuantiles = 301
AWAP.QR.signif <- AWAP.quantSignif$relQuantChange
for (i in 1:length(lon)) {
  for (j in 1:length(lat)) {
    for (q in 1:numQuantiles) {
      if (!is.na(AWAP.QR.signif[i,j,q])) {
        if (AWAP.QR.signif[i,j,q] <= AWAP.quantSignif$confints.QR[i,j,q,2] && AWAP.QR.signif[i,j,q] >= AWAP.quantSignif$confints.QR[i,j,q,1]) {
          AWAP.QR.signif[i,j,q] <- NA
        }
      }
    }
  }
}
# BRNS
numQuantiles = 301
BRNS.QR.signif <- BRNS.quantSignif$relQuantChange
for (i in 1:length(lon)) {
  for (j in 1:length(lat)) {
    for (q in 1:numQuantiles) {
      if (!is.na(BRNS.QR.signif[i,j,q])) {
        if (BRNS.QR.signif[i,j,q] <= BRNS.quantSignif$confints.QR[i,j,q,2] && BRNS.QR.signif[i,j,q] >= BRNS.quantSignif$confints.QR[i,j,q,1]) {
          BRNS.QR.signif[i,j,q] <- NA
        }
      }
    }
  }
}

# using only the significant grids
NNI.postonegchanges <- calc.postonegchanges(NATN.QR.signif)
AWAP.postonegchanges <- calc.postonegchanges(AWAP.QR.signif)
BOA.postonegchanges <- calc.postonegchanges(BRNS.QR.signif)
# Average changes in quantiles
NNI.AvgQuantChange <- apply(NATN.QR.signif, 3, mean, na.rm = T)
AWAP.AvgQuantChange <- apply(AWAP.QR.signif, 3, mean, na.rm = T)
BOA.AvgQuantChange <- apply(BRNS.QR.signif, 3, mean, na.rm = T)

# find out which quantiles have atleast 5% significant grids
# All day quantiles
NNI_percentSignif <- apply(NATN.QR.signif, 3, FUN = function(x) {sum(!is.na(x))/2502*100})
AWAP_percentSignif <- apply(AWAP.QR.signif, 3, FUN = function(x) {sum(!is.na(x))/2502*100})
BOA_percentSignif <- apply(BRNS.QR.signif, 3, FUN = function(x) {sum(!is.na(x))/2062*100})
# Wet day quantiles
NNI_percentSignif.wet <- apply(NATN.QR.signif.wet, 3, FUN = function(x) {sum(!is.na(x))/2502*100})
AWAP_percentSignif.wet <- apply(AWAP.QR.signif.wet, 3, FUN = function(x) {sum(!is.na(x))/2502*100})
BOA_percentSignif.wet <- apply(BRNS.QR.signif.wet, 3, FUN = function(x) {sum(!is.na(x))/2062*100})
# Now TRUE All day quantiles
# load the data for all day quantiles first
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/NATN_alldayquantSignif.RData")
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/AWAP_allquantSignif.RData")
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/BRNS_allquantsignif.RData")

numQuantiles = 301
# apply mask
NATN.alldayquantSignif$relQuantChange[mask_array] <- NA
NATN.QR.signif <- NATN.alldayquantSignif$relQuantChange
for (i in 1:length(lon)) {
  for (j in 1:length(lat)) {
    for (q in 1:numQuantiles) {
      if (!is.na(NATN.QR.signif[i,j,q])) {
        if (NATN.QR.signif[i,j,q] <= NATN.alldayquantSignif$confints.QR[i,j,q,2] && NATN.QR.signif[i,j,q] >= NATN.alldayquantSignif$confints.QR[i,j,q,1]) {
          NATN.QR.signif[i,j,q] <- NA
        }
      }
    }
  }
}
# AWAP
numQuantiles = 301
# apply mask
AWAP.allquantSignif$relQuantChange[mask_array] <- NA
AWAP.QR.signif <- AWAP.allquantSignif$relQuantChange
for (i in 1:length(lon)) {
  for (j in 1:length(lat)) {
    for (q in 1:numQuantiles) {
      if (!is.na(AWAP.QR.signif[i,j,q])) {
        if (AWAP.QR.signif[i,j,q] <= AWAP.allquantSignif$confints.QR[i,j,q,2] && AWAP.QR.signif[i,j,q] >= AWAP.allquantSignif$confints.QR[i,j,q,1]) {
          AWAP.QR.signif[i,j,q] <- NA
        }
      }
    }
  }
}
# BRNS
numQuantiles = 301
# apply mask
BRNS.allquantSignif$relQuantChange[mask_array] <- NA
BRNS.QR.signif <- BRNS.allquantSignif$relQuantChange
for (i in 1:length(lon)) {
  for (j in 1:length(lat)) {
    for (q in 1:numQuantiles) {
      if (!is.na(BRNS.QR.signif[i,j,q])) {
        if (BRNS.QR.signif[i,j,q] <= BRNS.allquantSignif$confints.QR[i,j,q,2] && BRNS.QR.signif[i,j,q] >= BRNS.allquantSignif$confints.QR[i,j,q,1]) {
          BRNS.QR.signif[i,j,q] <- NA
        }
      }
    }
  }
}

# using only the significant grids
NNI.postonegchanges <- calc.postonegchanges(NATN.QR.signif)
AWAP.postonegchanges <- calc.postonegchanges(AWAP.QR.signif)
BOA.postonegchanges <- calc.postonegchanges(BRNS.QR.signif)
# Average changes in quantiles
NNI.AvgQuantChange <- apply(NATN.alldayquantSignif$relQuantChange, 3, mean, na.rm = T)
AWAP.AvgQuantChange <- apply(AWAP.allquantSignif$relQuantChange, 3, mean, na.rm = T)
BOA.AvgQuantChange <- apply(BRNS.allquantSignif$relQuantChange, 3, mean, na.rm = T)

# find out which quantiles have atleast 5% significant grids
# All day quantiles
# there are 2502 non-missing grids in AWAP
# there are 2062 non-missing grids in BRNS and combined masks
NNI_percentSignif <- apply(NATN.QR.signif, 3, FUN = function(x) {sum(!is.na(x))/2062*100})
AWAP_percentSignif <- apply(AWAP.QR.signif, 3, FUN = function(x) {sum(!is.na(x))/2062*100})
BOA_percentSignif <- apply(BRNS.QR.signif, 3, FUN = function(x) {sum(!is.na(x))/2062*100})
# Wet day quantiles
NNI_percentSignif.wet <- apply(NATN.QR.signif.wet, 3, FUN = function(x) {sum(!is.na(x))/2062*100})
AWAP_percentSignif.wet <- apply(AWAP.QR.signif.wet, 3, FUN = function(x) {sum(!is.na(x))/2062*100})
BOA_percentSignif.wet <- apply(BRNS.QR.signif.wet, 3, FUN = function(x) {sum(!is.na(x))/2062*100})
# Now the plotting
cols <- c('#1b9e77','#d95f02','#7570b3')
layout(mat = matrix(c(1:5,5), ncol = 2, byrow = T), heights = c(1,1,0.1))
par(mar = c(4,5,0,0) + 0.1, oma = c(0,0,0.5,0) + 0.1, xpd = F, bg = NA)
xvals <- seq(0,1,length.out = 301)
yrange <- range(NNI.postonegchanges.wet[which(NNI_percentSignif.wet > 5)], 
                AWAP.postonegchanges.wet[which(AWAP_percentSignif.wet > 5)], 
                BOA.postonegchanges.wet[which(BOA_percentSignif.wet > 5)])
plot(xvals[which(NNI_percentSignif.wet > 5)], NNI.postonegchanges.wet[which(NNI_percentSignif.wet > 5)], 
     ylim = yrange, xlim = range(xvals),
     pch = 19, cex = 0.5, cex.axis = 1.5, cex.lab = 1.5, col = cols[1], xlab = "Quantiles", 
     ylab = "Relative difference in number of grids")
points(xvals[which(AWAP_percentSignif.wet > 5)], AWAP.postonegchanges.wet[which(AWAP_percentSignif.wet > 5)], 
       pch = 19, cex = 0.5, col = cols[2])
points(xvals[which(BOA_percentSignif.wet > 5)], BOA.postonegchanges.wet[which(BOA_percentSignif.wet > 5)], 
       pch = 19, cex = 0.5, col = cols[3])
abline(h = 0.0)
mtext("(a)", cex = 2, side = 3, line = -2, adj = 0)
#par(mar = c(4,6,2,2) + 0.1)
xvals <- seq(0,1,length.out = 301)
yrange = range(NNI.postonegchanges[which(NNI_percentSignif > 5)], 
               AWAP.postonegchanges[which(AWAP_percentSignif > 5)], 
               BOA.postonegchanges[which(BOA_percentSignif > 5)], na.rm = T)
plot(xvals[which(NNI_percentSignif > 5)], NNI.postonegchanges[which(NNI_percentSignif > 5)], 
     ylim = yrange, xlim = range(xvals), 
     pch = 19, cex = 0.5, cex.axis = 1.5, cex.lab = 1.5, col = cols[1], xlab = "Quantiles", 
     ylab = "Relative Difference in number of grids")
points(xvals[which(AWAP_percentSignif > 5)], AWAP.postonegchanges[which(AWAP_percentSignif > 5)], 
       pch = 19, cex = 0.5, col = cols[2])
points(xvals[which(BOA_percentSignif > 5)], BOA.postonegchanges[which(BOA_percentSignif > 5)], 
       pch = 19, cex = 0.5, col = cols[3])
abline(h = 0.0)
mtext("(b)", cex = 2, side = 3, line = -2, adj = 0)
# Average quantile change
# wet days
yrange <- range(NNI.AvgQuantChange.wet, 
                AWAP.AvgQuantChange.wet, 
                BOA.AvgQuantChange.wet)
plot(xvals, NNI.AvgQuantChange.wet, 
     ylim = yrange, xlim = range(xvals), 
     pch = 19, cex = 0.5, cex.axis = 1.5, cex.lab = 1.5, col = cols[1], xlab = "Quantiles", 
     ylab = "Average Quantile Change")
points(xvals, AWAP.AvgQuantChange.wet, 
       pch = 19, cex = 0.5, col = cols[2])
points(xvals, BOA.AvgQuantChange.wet, pch = 19, cex = 0.5, col = cols[3])
abline(h = 0.0)
mtext("(c)", cex = 2, side = 3, line = -2, adj = 0)
#par(mar = c(4,6,2,2) + 0.1)
# all days
yrange <- range(NNI.AvgQuantChange, 
               AWAP.AvgQuantChange, 
               BOA.AvgQuantChange, na.rm = T)
xvals <- seq(0,1,length.out = 301)
plot(xvals, NNI.AvgQuantChange, 
     ylim = yrange, xlim = range(xvals), 
     pch = 19, cex = 0.5, cex.axis = 1.5, cex.lab = 1.5, col = cols[1], xlab = "Quantiles", 
     ylab = "Average Quantile Change")
points(xvals, AWAP.AvgQuantChange, 
       pch = 19, cex = 0.5, col = cols[2])
points(xvals, BOA.AvgQuantChange, 
       pch = 19, cex = 0.5, col = cols[3])
abline(h = 0.0)
mtext("(d)", cex = 2, side = 3, line = -2, adj = 0)
par(fig = c(0,1,0,1), oma = c(0,0,0,0), mar = c(0,0,0,0), new = T)
plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n")
legend("bottom", legend = c("NNI", "AWAP", "BOA"), col = cols, 
      pch = 19, bty = "n", horiz = T, cex = 1.7)

library(sp)
library(raster)

poly.xy <- data.frame(x = c(142, 142, 120, 112, 130, 142), y = c(-10, -31, -31, -21, -10, -10))
poly <- Polygon(poly.xy)
polys <- Polygons(list(poly), 1)
sppoly <- SpatialPolygons(list(polys))
proj4string(sppoly) <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

# run above code chunks to get the NATN.alldayquantSignig object
# run previous code chunk to load data and create R objects

# turn NATN.alldayquantSignif$relQuantChange into raster
NATN.quantChange.rbrick <- flip(t(brick(NATN.alldayquantSignif$relQuantChange, xmn = min(lat), xmx = max(lat), 
                                          ymn = min(lon), ymx = max(lon), 
                                          crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84     +towgs84=0,0,0")), direction = 'y')
BRNS.quantChange.rbrick <- flip(t(brick(BRNS.allquantSignif$relQuantChange, xmn = min(lat), xmx = max(lat), 
                                          ymn = min(lon), ymx = max(lon), 
                                          crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84     +towgs84=0,0,0")), direction = 'y')
AWAP.quantChange.rbrick <- flip(t(brick(AWAP.allquantSignif$relQuantChange, xmn = min(lat), xmx = max(lat), 
                                          ymn = min(lon), ymx = max(lon), 
                                          crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84     +towgs84=0,0,0")), direction = 'y')

# mask NATN.quantChange object 
# mask one way
NATN.quantChange.freqinc <- mask(NATN.quantChange.rbrick, sppoly)
# mask the other way
NATN.quantChange.freqdec <- mask(NATN.quantChange.rbrick, sppoly, inverse = T)


# Average spatially
NATN.quantChange.freqinc.avg <- cellStats(NATN.quantChange.freqinc, 'mean', na.rm = T)
NATN.quantChange.freqdec.avg <- cellStats(NATN.quantChange.freqdec, 'mean', na.rm = T)

# mask AWAP.quantChange object 
# mask one way
AWAP.quantChange.freqinc <- mask(AWAP.quantChange.rbrick, sppoly)
# mask the other way
AWAP.quantChange.freqdec <- mask(AWAP.quantChange.rbrick, sppoly, inverse = T)


# Average spatially
AWAP.quantChange.freqinc.avg <- cellStats(AWAP.quantChange.freqinc, 'mean', na.rm = T)
AWAP.quantChange.freqdec.avg <- cellStats(AWAP.quantChange.freqdec, 'mean', na.rm = T)

# mask BRNS.quantChange object 
# mask one way
BRNS.quantChange.freqinc <- mask(BRNS.quantChange.rbrick, sppoly)
# mask the other way
BRNS.quantChange.freqdec <- mask(BRNS.quantChange.rbrick, sppoly, inverse = T)


# Average spatially
BRNS.quantChange.freqinc.avg <- cellStats(BRNS.quantChange.freqinc, 'mean', na.rm = T)
BRNS.quantChange.freqdec.avg <- cellStats(BRNS.quantChange.freqdec, 'mean', na.rm = T)

# Find how many grids inside and outside sppoly

# Convert AWAP.missing and BRNS.missing into a raster
AWAP.nonmissing <- array(data = 1, dim = dim(AWAP.missing))
AWAP.nonmissing[which(AWAP.missing)] <- NA
AWAP.nonmissing.ras <- flip(t(raster(AWAP.nonmissing, xmn = min(lat), xmx = max(lat), 
                                          ymn = min(lon), ymx = max(lon), 
                                          crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84     +towgs84=0,0,0")), direction = 'y')
# Now count the non missing grids in inc freq region
nonmissing.freqinc <- cellStats(mask(AWAP.nonmissing.ras, sppoly), 'sum', na.rm = T) #1171
nonmissing.freqdec <- cellStats(mask(AWAP.nonmissing.ras, sppoly, inverse = T), 'sum', na.rm = T) #1331

# do the same for BRNS
BRNS.nonmissing <- array(data = 1, dim = dim(BRNS.missing))
BRNS.nonmissing[which(BRNS.missing)] <- NA
BRNS.nonmissing.ras <- flip(t(raster(BRNS.nonmissing, xmn = min(lat), xmx = max(lat), 
                                          ymn = min(lon), ymx = max(lon), 
                                          crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84     +towgs84=0,0,0")), direction = 'y')
# Now count the non missing grids in inc freq region
nonmissing.freqinc <- cellStats(mask(BRNS.nonmissing.ras, sppoly), 'sum', na.rm = T) #906
nonmissing.freqdec <- cellStats(mask(BRNS.nonmissing.ras, sppoly, inverse = T), 'sum', na.rm = T) #1156

# Relative difference between number of grids showing positive changes and 
# those showing negative changes for all quantiles 

# Now create a sp object with only the significant grids (not all grids)
# turn NATN.alldayquantSignif$relQuantChange into raster
NATN.quantChange.Signif.rbrick <- flip(t(brick(NATN.QR.signif.wet, xmn = min(lat), xmx = max(lat), 
                                          ymn = min(lon), ymx = max(lon), 
                                          crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84     +towgs84=0,0,0")), direction = 'y')
BRNS.quantChange.Signif.rbrick <- flip(t(brick(BRNS.QR.signif.wet, xmn = min(lat), xmx = max(lat), 
                                          ymn = min(lon), ymx = max(lon), 
                                          crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84     +towgs84=0,0,0")), direction = 'y')
AWAP.quantChange.Signif.rbrick <- flip(t(brick(AWAP.QR.signif.wet, xmn = min(lat), xmx = max(lat), 
                                          ymn = min(lon), ymx = max(lon), 
                                          crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84     +towgs84=0,0,0")), direction = 'y')

# mask NATN.quantChange.Signif object 
# mask one way
NATN.quantChange.Signif.freqinc <- mask(NATN.quantChange.Signif.rbrick, sppoly)
# mask the other way
NATN.quantChange.Signif.freqdec <- mask(NATN.quantChange.Signif.rbrick, sppoly, inverse = T)

# mask AWAP.quantChange.Signif object 
# mask one way
AWAP.quantChange.Signif.freqinc <- mask(AWAP.quantChange.Signif.rbrick, sppoly)
# mask the other way
AWAP.quantChange.Signif.freqdec <- mask(AWAP.quantChange.Signif.rbrick, sppoly, inverse = T)

# mask BRNS.quantChange.Signif object 
# mask one way
BRNS.quantChange.Signif.freqinc <- mask(BRNS.quantChange.Signif.rbrick, sppoly)
# mask the other way
BRNS.quantChange.Signif.freqdec <- mask(BRNS.quantChange.Signif.rbrick, sppoly, inverse = T)

# Percentage of significant grids
# NATN
numSignif_freqdec_NATN <- apply(!is.na(getValues(NATN.quantChange.Signif.freqdec)), 2, FUN = function(n) {sum(n)/1156*100})
numSignif_freqinc_NATN <- apply(!is.na(getValues(NATN.quantChange.Signif.freqinc)), 2, FUN = function(n) {sum(n)/906*100})
# AWAP
numSignif_freqdec_AWAP <- apply(!is.na(getValues(AWAP.quantChange.Signif.freqdec)), 2, FUN = function(n) {sum(n)/1156*100})
numSignif_freqinc_AWAP <- apply(!is.na(getValues(AWAP.quantChange.Signif.freqinc)), 2, FUN = function(n) {sum(n)/906*100})
# BRNS
numSignif_freqdec_BRNS <- apply(!is.na(getValues(BRNS.quantChange.Signif.freqdec)), 2, FUN = function(n) {sum(n)/1156*100})
numSignif_freqinc_BRNS <- apply(!is.na(getValues(BRNS.quantChange.Signif.freqinc)), 2, FUN = function(n) {sum(n)/906*100})

# areas where wet day freq has decreased
# NATN
gtzero_freqdec_NATN <- apply(getValues(NATN.quantChange.freqdec), 2, FUN = function(n) {length(which(n > 0))})
ltzero_freqdec_NATN <- apply(getValues(NATN.quantChange.freqdec), 2, FUN = function(n) {length(which(n < 0))})
relChange_freqdec_NATN <- (gtzero_freqdec_NATN - ltzero_freqdec_NATN)/1156
# AWAP
gtzero_freqdec_AWAP <- apply(getValues(AWAP.quantChange.freqdec), 2, FUN = function(n) {length(which(n > 0))})
ltzero_freqdec_AWAP <- apply(getValues(AWAP.quantChange.freqdec), 2, FUN = function(n) {length(which(n < 0))})
relChange_freqdec_AWAP <- (gtzero_freqdec_AWAP - ltzero_freqdec_AWAP)/1156
# BRNS
gtzero_freqdec_BRNS <- apply(getValues(BRNS.quantChange.freqdec), 2, FUN = function(n) {length(which(n > 0))})
ltzero_freqdec_BRNS <- apply(getValues(BRNS.quantChange.freqdec), 2, FUN = function(n) {length(which(n < 0))})
relChange_freqdec_BRNS <- (gtzero_freqdec_BRNS - ltzero_freqdec_BRNS)/1156

# areas where wet day freq has increased
# NATN
gtzero_freqinc_NATN <- apply(getValues(NATN.quantChange.freqinc), 2, FUN = function(n) {length(which(n > 0))})
ltzero_freqinc_NATN <- apply(getValues(NATN.quantChange.freqinc), 2, FUN = function(n) {length(which(n < 0))})
relChange_freqinc_NATN <- (gtzero_freqinc_NATN - ltzero_freqinc_NATN)/906
# AWAP
gtzero_freqinc_AWAP <- apply(getValues(AWAP.quantChange.freqinc), 2, FUN = function(n) {length(which(n > 0))})
ltzero_freqinc_AWAP <- apply(getValues(AWAP.quantChange.freqinc), 2, FUN = function(n) {length(which(n < 0))})
relChange_freqinc_AWAP <- (gtzero_freqinc_AWAP - ltzero_freqinc_AWAP)/906
# BRNS
gtzero_freqinc_BRNS <- apply(getValues(BRNS.quantChange.freqinc), 2, FUN = function(n) {length(which(n > 0))})
ltzero_freqinc_BRNS <- apply(getValues(BRNS.quantChange.freqinc), 2, FUN = function(n) {length(which(n < 0))})
relChange_freqinc_BRNS <- (gtzero_freqinc_BRNS - ltzero_freqinc_BRNS)/906
# quantile(c(NATN.QR.signif.freqinc.avg, AWAP.QR.signif.freqinc.avg, BRNS.QR.signif.freqinc.avg),
#         probs = seq(0, 1, by = 0.01), na.rm = T)

quants <- seq(0, 1, length.out = 301)
cols <- c('#1b9e77','#d95f02','#7570b3')
par(mar = c(5,5,0,0) + 0.1, oma = c(0,0,0,0) + 0.1, bg = NA)
layout(mat = matrix(c(1:5,5), ncol = 2, byrow = T), heights = c(1,1,0.1))

# relChange
# areas where wet day freq has inc
yrange = range(c(relChange_freqinc_NATN[which(numSignif_freqinc_NATN > 5)], 
                 relChange_freqinc_BRNS[which(numSignif_freqinc_AWAP > 5)], 
                 relChange_freqinc_AWAP[which(numSignif_freqinc_BRNS > 5)]), na.rm = T)
plot(quants[which(numSignif_freqinc_NATN > 5)], relChange_freqinc_NATN[which(numSignif_freqinc_NATN > 5)], 
     ylim = yrange, xlim = range(quants), pch = 19, cex = 0.5, col = cols[1],
     xlab = "Quantiles", ylab = "Relative Difference in number of grids", cex.axis = 1.5, cex.lab = 1.5)
points(quants[which(numSignif_freqinc_AWAP > 5)], relChange_freqinc_AWAP[which(numSignif_freqinc_AWAP > 5)], 
       col = cols[2], pch = 19, cex = 0.5)
points(quants[which(numSignif_freqinc_BRNS > 5)], relChange_freqinc_BRNS[which(numSignif_freqinc_BRNS > 5)], 
       col = cols[3], pch = 19, cex = 0.5)
abline(h = 0)
text(0,yrange[2],"(a)", cex = 2.2)
# areas where wet day freq has dec
yrange = range(c(relChange_freqdec_NATN[which(numSignif_freqdec_NATN > 5)], 
                 relChange_freqdec_BRNS[which(numSignif_freqdec_AWAP > 5)], 
                 relChange_freqdec_AWAP[which(numSignif_freqdec_BRNS > 5)]), na.rm = T)
plot(quants[which(numSignif_freqdec_NATN > 5)], relChange_freqdec_NATN[which(numSignif_freqdec_NATN > 5)], 
     ylim = yrange, xlim = range(quants), pch = 19, cex = 0.5, col = cols[1],
     xlab = "Quantiles", ylab = "Relative Difference in number of grids", cex.axis = 1.5, cex.lab = 1.5)
points(quants[which(numSignif_freqdec_AWAP > 5)], relChange_freqdec_AWAP[which(numSignif_freqdec_AWAP > 5)], 
       col = cols[2], pch = 19, cex = 0.5)
points(quants[which(numSignif_freqdec_BRNS > 5)], relChange_freqdec_BRNS[which(numSignif_freqdec_BRNS > 5)], 
      col = cols[3], pch = 19, cex = 0.5)
abline(h = 0)
text(0,yrange[2],"(b)", cex = 2.2)

# Area average
# areas where we t day freq has inc
quants <- seq(0, 1, length.out = 301)
yrange <- range(c(NATN.quantChange.freqinc.avg,
                  AWAP.quantChange.freqinc.avg,
                  BRNS.quantChange.freqinc.avg), na.rm = T)
plot(quants, NATN.quantChange.freqinc.avg,
     ylim = yrange, xlim = range(quants), col = cols[1],
     xlab = "Quantiles", ylab = "Average Quantile Change", cex.axis = 1.5, cex.lab = 1.5, pch = 19, cex = 0.5)
points(quants, AWAP.quantChange.freqinc.avg,
       col = cols[2], pch = 19, cex = 0.5)
points(quants, BRNS.quantChange.freqinc.avg,
       col = cols[3], pch = 19, cex = 0.5)
abline(h = 0)
text(0,yrange[2],"(c)", cex = 2.2)
# areas where wet day freq has dec
yrange <- range(c(NATN.quantChange.freqdec.avg,
                  AWAP.quantChange.freqdec.avg,
                  BRNS.quantChange.freqdec.avg), na.rm = T)
plot(quants, NATN.quantChange.freqdec.avg,
     ylim = yrange, xlim = range(quants), col = cols[1],
     xlab = "Quantiles", ylab = "Average Quantile Change", cex.axis = 1.5, cex.lab = 1.5, pch = 19, cex = 0.5)
points(quants, AWAP.quantChange.freqdec.avg, col = cols[2], pch = 19, cex = 0.5)
points(quants, BRNS.quantChange.freqdec.avg, col = cols[3], pch = 19, cex = 0.5)
abline(h = 0)
text(0,yrange[2],"(d)", cex = 2.2)

# legend
par(fig = c(0,1,0,1), oma = c(0,0,0,0), mar = c(0,0,0,0), new = T)
plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n")
legend("bottom", legend = c("NNI", "AWAP", "BOA"), col = cols, 
      pch = 19, bty = "n", horiz = T, cex = 1.7)

# determine brks and col for plotting
zrange <- range(AWAP.wetratiodiff.Seas, BRNS.wetratiodiff.Seas, NATN.wetratiodiff.Seas, na.rm = T)
brks <- seq(-2, 2, length.out = 9)
        #seq(-2,2,0.2)
zbrks <- brks; zbrks[1] <- zrange[1]; zbrks[length(zbrks)] <- zrange[2]
col <- c('#8c510a','#bf812d','#dfc27d','#f6e8c3','#c7eae5','#80cdc1','#35978f','#01665e')
      #create.col.list(brks, 0)
# Map of region division on total frequency change
# load all relevant functions and data and objects first
par(mar = c(0,0,0,0) + 0.1, bg = NA)
layout(mat = matrix(c(1,2), ncol = 1, byrow = T), heights = c(1,0.15))
plot.wetratiodiff(NATN.wetratiodiff, brks, col, zrange, AWAP.missing)
plot(sppoly, add = T, lwd = 7)
par(mar = c(0,5,0,5) + 0.1, bg = NA)
colkey(col=col, side=1, clim=range(brks), width=8, at=brks, labels = paste(round(zbrks,1)))

# Convert back into arrays so we can use pre-written functions
# NATN
NATN.QR.signif.freqinc.array <- array(apply(as.array(NATN.QR.signif.freqinc), 3, t), dim = c(length(lon), length(lat), 301))
NATN.QR.signif.freqinc.array <- NATN.QR.signif.freqinc.array[,length(lat):1,]
# AWAP
AWAP.QR.signif.freqinc.array <- array(apply(as.array(AWAP.QR.signif.freqinc), 3, t), dim = c(length(lon), length(lat), 301))
AWAP.QR.signif.freqinc.array <- AWAP.QR.signif.freqinc.array[,length(lat):1,]
# BRNS
BRNS.QR.signif.freqinc.array <- array(apply(as.array(BRNS.QR.signif.freqinc), 3, t), dim = c(length(lon), length(lat), 301))
BRNS.QR.signif.freqinc.array <- BRNS.QR.signif.freqinc.array[,length(lat):1,]

# NATN
NATN.QR.signif.freqdec.array <- array(apply(as.array(NATN.QR.signif.freqdec), 3, t), dim = c(length(lon), length(lat), 301))
NATN.QR.signif.freqdec.array <- NATN.QR.signif.freqdec.array[,length(lat):1,]
# AWAP
AWAP.QR.signif.freqdec.array <- array(apply(as.array(AWAP.QR.signif.freqdec), 3, t), dim = c(length(lon), length(lat), 301))
AWAP.QR.signif.freqdec.array <- AWAP.QR.signif.freqdec.array[,length(lat):1,]
# BRNS
BRNS.QR.signif.freqdec.array <- array(apply(as.array(BRNS.QR.signif.freqdec), 3, t), dim = c(length(lon), length(lat), 301))
BRNS.QR.signif.freqdec.array <- BRNS.QR.signif.freqdec.array[,length(lat):1,]
calc.postonegchanges <- function(relquantchange, totGrids) {
  return(apply(relquantchange, 3, FUN = function(a) {
    #2770 is the number of land grids for AUSTRALIA
    return((length(which(a > 0)) - length(which(a < 0)))/totGrids)
  }))
}
# calculate rel diff in num of pos grids and neg grids
# NATN
NATN.freqinc.postonegchanges <- calc.postonegchanges(NATN.QR.signif.freqinc.array, 1171)
NATN.freqdec.postonegchanges <- calc.postonegchanges(NATN.QR.signif.freqdec.array, 1331)
# AWAP
AWAP.freqinc.postonegchanges <- calc.postonegchanges(AWAP.QR.signif.freqinc.array, 1171)
AWAP.freqdec.postonegchanges <- calc.postonegchanges(AWAP.QR.signif.freqdec.array, 1331)
# BRNS
BRNS.freqinc.postonegchanges <- calc.postonegchanges(BRNS.QR.signif.freqinc.array, 906)
BRNS.freqdec.postonegchanges <- calc.postonegchanges(BRNS.QR.signif.freqdec.array, 1156)
# Plotting
quantile(c(NATN.freqinc.postonegchanges, AWAP.freqinc.postonegchanges, BRNS.freqinc.postonegchanges),
         probs = seq(0, 1, by = 0.01), na.rm = T)
yrange = c(-0.03, 0.3)
quants <- seq(0, 1, length.out = 301)
plot(quants, NATN.freqinc.postonegchanges, ylim = yrange, type = "l")
lines(quants, AWAP.freqinc.postonegchanges, col = 'red')
lines(quants, BRNS.freqinc.postonegchanges, col = 'blue')
abline(h = 0)

# freqdec
quantile(c(NATN.freqdec.postonegchanges, AWAP.freqdec.postonegchanges, BRNS.freqdec.postonegchanges),
         probs = seq(0, 1, by = 0.01), na.rm = T)
yrange = c(-0.55, 0.7)
quants <- seq(0, 1, length.out = 301)
plot(quants, NATN.freqdec.postonegchanges, ylim = yrange, type = "l")
lines(quants, AWAP.freqdec.postonegchanges, col = 'red')
lines(quants, BRNS.freqdec.postonegchanges, col = 'blue')
abline(h = 0)

Even in regions with a decrease in intensity there are more grids with an increase in quantile for quantiles greater than around 0.9.

Mean SDII (intensity) plots for comparison with wet-day quantiles

Run the setup code chunk first!

# NATN
# Import SDII data
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/NatgridBarnes/Prcp_daily_natgrid_0.5x0.5deg_1958-2013_land_sdii.nc")
NNI.SDII <- var.get.nc(nc, "Prcp")
close.nc(nc)
# calculate the relative change in the intensity
NNI.SDII.1 <- apply(NNI.SDII[,,1:28], c(1,2), mean, na.rm = T)
NNI.SDII.2 <- apply(NNI.SDII[,,29:56], c(1,2), mean, na.rm = T)
NNI.SDII.relchange <- (NNI.SDII.2 - NNI.SDII.1)/apply(NNI.SDII[,,1:56], c(1,2), mean, na.rm = T)
#AWAP
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/AWAP/AWAP_PRCP_1958-2015_0.5deg_remapcon2_land_noneg_sdii_invertlat.nc")
AWAP.SDII <- var.get.nc(nc, "pre")
close.nc(nc)
# calculate the relative change in the intensity
AWAP.SDII.1 <- apply(AWAP.SDII[,,1:28], c(1,2), mean, na.rm = T)
AWAP.SDII.2 <- apply(AWAP.SDII[,,29:56], c(1,2), mean, na.rm = T)
AWAP.SDII.relchange <- (AWAP.SDII.2 - AWAP.SDII.1)/apply(AWAP.SDII[,,1:56], c(1,2), mean, na.rm = T)
#BRNS
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/NatgridBarnes/Prcp_daily_barnesgrid_0.5x0.5deg_1958-2013_Req1_land_sdii.nc")
BOA.SDII <- var.get.nc(nc, "Prcp")
close.nc(nc)
# calculate the relative change in the intensity
BOA.SDII.1 <- apply(BOA.SDII[,,1:28], c(1,2), mean, na.rm = T)
BOA.SDII.2 <- apply(BOA.SDII[,,29:56], c(1,2), mean, na.rm = T)
BOA.SDII.relchange <- (BOA.SDII.2 - BOA.SDII.1)/apply(BOA.SDII[,,1:56], c(1,2), mean, na.rm = T)

#plotting
plot.relChange <- function(relChange, lon, lat, mask, data.name = "NNI") {
  dim = dim(relChange)
  zrange <- range(relChange, na.rm = T)
  brks <- seq(-0.09,0.15, by = 0.03)
  zbrks <- brks; zbrks[1] <- min(brks[1],zrange[1]); zbrks[length(brks)] <- max(brks[length(brks)],zrange[2]) # since 0.3 < zrange[1]
  col <- create.col.list(brks,0)
  col[1] <- "indianred1"
  layout(mat = matrix(c(1,2), nrow = 2, byrow = T), heights = c(1,0.2))
  par(mar = c(5,4,2,2)+0.1)
    image(lon,lat,relChange, 
          breaks = zbrks, xaxt = 'n', yaxt = 'n', xlab = paste(data.name, "relative change in Intensity
  between 1958-85 and 1986-2013"), cex.lab = 1.4, col = col, ylab = "")
    image(lon, lat, mask, col="gray", add = T)
    map("world", "Australia", add = T)
    par(mar = c(0,4,0,2)+0.1)
  colkey(col = col, side = 3, clim = range(brks), width = 5, at = brks, labels = zbrks, cex.axis = 1.6)
}

plot.relChange(NNI.SDII.relchange, lon, lat, AWAP.missing, "NNI")
plot.relChange(AWAP.SDII.relchange, lon, lat, AWAP.missing, "AWAP")
plot.relChange(BOA.SDII.relchange, lon, lat, BRNS.missing, "BOA")

# NATN
# Import SDII data
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/NatgridBarnes/Prcp_daily_natgrid_0.5x0.5deg_1958-2013_land_dailyavgprcp.nc")
NNI.SDII <- var.get.nc(nc, "Prcp")
close.nc(nc)
# calculate the relative change in the intensity
NNI.SDII.1 <- apply(NNI.SDII[,,1:28], c(1,2), mean, na.rm = T)
NNI.SDII.2 <- apply(NNI.SDII[,,29:56], c(1,2), mean, na.rm = T)
NNI.SDII.relchange <- (NNI.SDII.2 - NNI.SDII.1)/apply(NNI.SDII[,,1:56], c(1,2), mean, na.rm = T)
#AWAP
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/AWAP/AWAP_PRCP_1958-2015_0.5deg_remapcon2_land_noneg_dailyavgprcp_invertlat.nc")
AWAP.SDII <- var.get.nc(nc, "pre")
close.nc(nc)
# calculate the relative change in the intensity
AWAP.SDII.1 <- apply(AWAP.SDII[,,1:28], c(1,2), mean, na.rm = T)
AWAP.SDII.2 <- apply(AWAP.SDII[,,29:56], c(1,2), mean, na.rm = T)
AWAP.SDII.relchange <- (AWAP.SDII.2 - AWAP.SDII.1)/apply(AWAP.SDII[,,1:56], c(1,2), mean, na.rm = T)
#BRNS
nc <- open.nc("/srv/ccrc/data45/z3289452/DailyPRCP_netCDFs/NatgridBarnes/Prcp_daily_barnesgrid_0.5x0.5deg_1958-2013_Req1_land_dailyavgprcp.nc")
BOA.SDII <- var.get.nc(nc, "Prcp")
close.nc(nc)
# calculate the relative change in the intensity
BOA.SDII.1 <- apply(BOA.SDII[,,1:28], c(1,2), mean, na.rm = T)
BOA.SDII.2 <- apply(BOA.SDII[,,29:56], c(1,2), mean, na.rm = T)
BOA.SDII.relchange <- (BOA.SDII.2 - BOA.SDII.1)/apply(BOA.SDII[,,1:56], c(1,2), mean, na.rm = T)

#plotting
plot.relChange <- function(relChange, lon, lat, mask, data.name = "NNI") {
  dim = dim(relChange)
  zrange <- range(relChange, na.rm = T)
  brks <- seq(-0.10,0.30, by = 0.05)
  zbrks <- brks; zbrks[1] <- min(brks[1],zrange[1]); zbrks[length(brks)] <- max(brks[length(brks)],zrange[2]) # since 0.3 < zrange[1]
  col <- create.col.list(brks,0)
  col[1] <- "indianred1"
  layout(mat = matrix(c(1,2), nrow = 2, byrow = T), heights = c(1,0.2))
  par(mar = c(5,4,2,2)+0.1)
    image(lon,lat,relChange, 
          breaks = zbrks, xaxt = 'n', yaxt = 'n', xlab = paste(data.name, "relative change in daily avg precipitation
  between 1958-85 and 1986-2013"), cex.lab = 1.4, col = col, ylab = "")
    image(lon, lat, mask, col="gray", add = T)
    map("world", "Australia", add = T)
    par(mar = c(0,4,0,2)+0.1)
  colkey(col = col, side = 3, clim = range(brks), width = 5, at = brks, labels = zbrks, cex.axis = 1.6)
}

plot.relChange(NNI.SDII.relchange, lon, lat, AWAP.missing, "NNI")
plot.relChange(AWAP.SDII.relchange, lon, lat, AWAP.missing, "AWAP")
plot.relChange(BOA.SDII.relchange, lon, lat, BRNS.missing, "BOA")

# load WET_NATN_1 and WET_NATN_2
# Create and compare PDFs at lon = 123 (index = 23) and lat = -18 (index = 52)
TS1 <- WET_AWAP_1[80, 40,] #WET_AWAP_1  #WET_NATN_1[23, 52,]
TS2 <- WET_AWAP_2[80, 40,] #WET_AWAP_2  #WET_NATN_2[23, 52,]
start <- min(which(sort(TS1) > 0), which(sort(TS2) > 0))
TS1 <- TS1[start:length(TS1)]
TS2 <- TS2[start:length(TS2)]

# Create PDF
density_TS1 <- density(TS1, na.rm = T)
density_TS2 <- density(TS2, na.rm = T)

# Using quantiles
quant1 <- quantile(TS1, probs = seq(0, 1, length.out = 100), na.rm = T)
quant2 <- quantile(TS2, probs = seq(0, 1, length.out = 100), na.rm = T)

cols <- c('#d95f02','#7570b3')
bgs <- paste0(cols, "77")
# density plot
yrange <- range(density_TS1$y, density_TS2$y)
plot(density_TS1, col = cols[1], main = "Precipitation density", xlim = c(0,20), ylim = yrange)
#polygon(density_TS1, col = bgs[1])
lines(density_TS2, col = cols[2])
polygon(density_TS2, col = bgs[2])
legend("topright", legend = c("Period 1", "Period 2"), col = cols, lwd = 4, bty = "n", cex = 2)

# probability plot
breaks <- seq(0, max(TS1, TS2, na.rm = T), length.out = 1000)
hist1 <- hist(TS1, breaks = breaks)
counts1 <- hist1$counts
counts2 <- hist(TS2, breaks = breaks)$counts
smoothed1 <- smooth.spline(counts1, df = 999)
smoothed2 <- smooth.spline(counts2, df = 999)


# plot(breaks[2:length(breaks)], counts1, col = cols[1], type = "l")
# lines(breaks[2:length(breaks)], counts2, col = cols[2])

plot(smoothed1, col = cols[1], type = "l")
lines(smoothed2, col = cols[2])

Check for field significance of positive and negative change in quantiles

#################################################################################################################
# Null distribution of number of positive grids and number of negative grids
#################################################################################################################

quantiles <- round(seq(0, 1, length.out = 301), 3)
plotted_quantiles <- c(seq(0.1, 0.9, by = 0.1), seq(0.95, 1.0, 0.01))
plotted_quantiles.index <- c(seq(31,271, by = 30), seq(286, 301, by = 3))

# load the data
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/NATN_wetdayquantSignif_fldSignifNullDist.RData")

numPos.confints <- t(apply(NATN.fldSignifNullDist$numpos.nullDist, 2, function(d) quantile(d, probs=c(0.05,0.95))))
numNeg.confints <- t(apply(NATN.fldSignifNullDist$numneg.nullDist, 2, function(d) quantile(d, probs=c(0.05,0.95))))

# Now count the number of positive and negative grids for each quantile
numposgrids <- sapply(1:301, FUN = function(n) {length(which(NATN.wetdayquantSignif$relQuantChange[,,n] > 0))})
numneggrids <- sapply(1:301, FUN = function(n) {length(which(NATN.wetdayquantSignif$relQuantChange[,,1] < 0))})

# Fld significant pos/neg change quantiles
# pos
# upper
plotted_quantiles[which(plotted_quantiles.index %in% which(numposgrids > numPos.confints[,2]))]
##  [1] 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 0.95 0.96 0.97 0.98 0.99
## [15] 1.00
# lower
plotted_quantiles[which(plotted_quantiles.index %in% which(numposgrids < numPos.confints[,1]))]
## numeric(0)
# neg
# upper
plotted_quantiles[which(plotted_quantiles.index %in% which(numneggrids > numNeg.confints[,2]))]
## numeric(0)
# lower
plotted_quantiles[which(plotted_quantiles.index %in% which(numneggrids < numNeg.confints[,1]))]
##  [1] 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 0.95 0.96 0.97 0.98 0.99
# Now do this for the all day quantiles
load("/srv/ccrc/data45/z3289452/GIT_REPOS/PhD_repository/RData/NATN_alldayquantSignif_fldSignifNullDist.RData")

numPos.confints <- t(apply(NATN.fldSignifNullDist$numpos.nullDist, 2, function(d) quantile(d, probs=c(0.05,0.95))))
numNeg.confints <- t(apply(NATN.fldSignifNullDist$numneg.nullDist, 2, function(d) quantile(d, probs=c(0.05,0.95))))

# Now count the number of positive and negative grids for each quantile
numposgrids <- sapply(1:301, FUN = function(n) {length(which(NATN.alldayquantSignif$relQuantChange[,,n] > 0))})
numneggrids <- sapply(1:301, FUN = function(n) {length(which(NATN.alldayquantSignif$relQuantChange[,,1] < 0))})

# Fld significant pos/neg change quantiles
# pos
# upper
plotted_quantiles[which(plotted_quantiles.index %in% which(numposgrids > numPos.confints[,2]))]
##  [1] 0.10 0.50 0.60 0.70 0.80 0.90 0.95 0.96 0.97 0.98 0.99 1.00
# lower
plotted_quantiles[which(plotted_quantiles.index %in% which(numposgrids < numPos.confints[,1]))]
## numeric(0)
# neg
# upper
plotted_quantiles[which(plotted_quantiles.index %in% which(numneggrids > numNeg.confints[,2]))]
## numeric(0)
# lower
plotted_quantiles[which(plotted_quantiles.index %in% which(numneggrids < numNeg.confints[,1]))]
##  [1] 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 0.95 0.96 0.97 0.98 0.99 1.00