##From the LAT-LON file sent by Nick, Using ERA5 data, compute SPEI for each lat-lon coordinate ##NOTES: ##To compute we need "monthly climatic water balance (precipitation minus potential evapotranspiration)....in mm" #(from https://cran.r-project.org/web/packages/SPEI/SPEI.pdf) #This requires computing PET and here is done below using the Thornthwaite method that uses Tave and latitude. #Temperature is required in degrees Celcius and Precip as accumulated monthly total ##The script does both and converts degK to degC and converts to mm and acumulates to monthly ##1. SPEI ##According to the classifications of McKee et al. (1993) and Paulo et al. (2012) for the SPEI drought classes, we have 5 classes: ## > -0.5 Non-Drought (Class 1) ## -0.5 and -1 Mild drought (Class 2) ## -1.5 and -1 Moderate drought (Class 3) ## -2 and -1.5 Severe drought (Class 4) <-- WE TEST FOR SEVERE AND EXTREME ## < -2 Extreme drought (Class 5) #SPEI and SPI have teh same index threshold classification ##WMO <-- USING THIS ** for thresholds ##2.0+ extremely wet ##1.5 to 1.99 very wet** ##1.0 to 1.49 moderately wet ##-.99 to .99 near normal ##-1.0 to -1.49 moderately dry ##-1.5 to -1.99 severely dry** ##-2 and less extremely dry ##2. Hydrological parameters (From ERA5 documentation on monthly data) ##The hydrological parameters are in units of "m of water per day" and so they should be multiplied by 1000 to convert to kgm-2day-1 or mmday-1. ##Which means they do still need to be accumulated to monthly accumulations to compute SPEI ##Where are we working setwd("./Scripts") ##Packages to install #install.packages("tidyverse") #install.packages("tidync") library(tidync) library(RNetCDF) library(ncdf4) #install.packages("SPEI") #library(spei) library(SPEI) library("stringr") library(ncdf.tools) library(tools) library(lattice) library(lubridate) #install.packages("boot") library(boot) options(scipen = 999) #So numbers are not in scientific notation #####CONFIGURE SECTION##### ##Set up resolution and time period res <- "9km" timeperiod <- "1988-2018" ##Change to 2008-2018 for this analysis syear <- 1988 ##For if you want to plot a time series later; Change to 2008 for 2008 analysis eyear <- 2018 ##For if you want to plot a time series later filenamevarppt <- "PPT" filenamevart2m <- "T2M" #####END CONFIGURE SECTION##### #0.Where to write the output data and constructing the filename outdir <- "../OutData" outfilename1 <- "AllRespn.spei_thorn" outfilename2 <- res outfilename3 <- timeperiod outfilename4 <- "txt" outfile1 <-paste(outfilename1,outfilename2,outfilename3,outfilename4,sep=".") print(outfile1) ##AfroBarometer files corlatlondir <- "../../Data/AfroBarom/Data" corlatlonfile <- "Afro.FinalVersion.CT.csv" ##Lat-lon coordinates of all responses file <- paste(corlatlondir,corlatlonfile,sep = "/") print(file) ##Read in the station data print(paste("Reading in the station data for",file,sep=" ")) stn <-read.csv(file, header = FALSE, sep=",") head(stn) ##Convert to a data frame stn.df <- as.data.frame(stn, header = TRUE) head(stn.df) ##Read in the netcdf file #1. ERA5 Data...compute PET using Thornthwaite method ERAindirroot <- "../../Data/ERA5" ERAindir <- paste(ERAindirroot,res,sep="/") infilename1 <- "Afr.ERA5.Mon" infilename2ppt <- filenamevarppt infilename2t2m <- filenamevart2m infilename3 <- timeperiod infilename4 <- res infilename5 <- "nc" fnameppt <- paste(infilename1,infilename2ppt,infilename3,infilename4,infilename5,sep=".") fnamet2m <- paste(infilename1,infilename2t2m,infilename3,infilename4,infilename5,sep=".") print(fnameppt) print(fnamet2m) ncinppt <- nc_open(paste(ERAindir, fnameppt,sep="/")) ncint2m <- nc_open(paste(ERAindir, fnamet2m,sep="/")) #str(ncin) eravarppt <- ncinppt$var[[1]]$name eravart2m <- ncint2m$var[[1]]$name print(eravarppt) print(eravart2m) #2. Time t <- ncvar_get(ncinppt, "time") tunits <- ncatt_get(ncinppt, "time", "units") tunits1 <- ncatt_get(ncinppt, "time", "units") tustr <- strsplit(tunits1$value, " ") tunits2 <- paste(unlist(tustr)[3],unlist(tustr)[4], sep =" ") # convert time units from "hours since..." to sensible notation & create dataframe eradatadates <- as.POSIXct(t*3600,origin=tunits2) #origin='1900-01-01 00:00' #3. Extract ERA5 ppt data at location ##LOOPING THRU ROWS print("STARTING LOOP THROUGH ROWS ie LAT-LONS") for(row in 1:nrow(stn.df)) { print("STARTING LOOP THROUGH ROWS, ID and LAT-LON are:") ##Tests # row <- 1 # row <- 701 # row <- 29592 id <- as.character(stn.df[row, "V2"]) #Do as.character else levels are included and then loops through all levels. Levels are every row's ID lat <- as.numeric(as.character(stn.df[row, "V4"])) #Need to convert character to numeric lon <- as.numeric(as.character(stn.df[row, "V3"])) print(paste(id,lat,lon,sep=",")) ##To deal with Rodrigues Island NAs... skipping.... if(is.na(lat)) { next } # get values at location; use ncvar_get() where "start(xdim,ydim,zdim,tdim)" extracts the lat-lon coordinate across all time steps print("Figuring out lat lon coordinate to extract for and extracting") ##PPT era.latlon.ppt <- ncvar_get(ncinppt, eravarppt, start= c(which.min(abs(ncinppt$dim$longitude$vals - lon)), # look for closest long - X dim which.min(abs(ncinppt$dim$latitude$vals - lat)), # look for closest lat - Y dim 1), # Z dim count = c(1,1,-1) # T dim: count '-1' means 'all values along that dimension'that dimension' ) ##TAVE era.latlon.t2m <- ncvar_get(ncint2m, eravart2m, start= c(which.min(abs(ncint2m$dim$longitude$vals - lon)), # look for closest long - X dim which.min(abs(ncint2m$dim$latitude$vals - lat)), # look for closest lat - Y dim 1), # Z dim count = c(1,1,-1) # T dim: count '-1' means 'all values along that dimension'that dimension' ) ##Compute PET #Construct a dataframe of ppt, t2m, then add PET and cwb #First for rainfall Convert from M to mm and accumulate monthly totals, because #computing PET needs total monthly rainfall (and temperature in degrees Celcius) era.latlon.mm.a <- era.latlon.ppt*1000 ##ACCUMULATE TO MONTHLY TOTAL #Data are in mm/day/month so need to be accumulated up to monthly total ## library(lubridate) <-- DONE AT THE TOP mm.day.mon <- data.frame(dates= eradatadates, obs = era.latlon.mm.a) #str(mm.day.mon) daysinmon <- as.numeric(days_in_month(mm.day.mon$dates)) mm.day.mon$obs <- mm.day.mon$obs*daysinmon era.latlon.mm <- mm.day.mon$obs era.latlon.mm ##Checking the number of timesteps are still the same numrows.t2m <- length(era.latlon.t2m) numrows.ppt <- length(era.latlon.mm) numrows.t2m/12 ## Making the dataframe and computing PET and CWB df1 <- as.data.frame(era.latlon.mm) df.t2m <- as.data.frame(era.latlon.t2m-273) #Thornthwaite needs Celcius df1$t2m <- df.t2m head(df1) #Compute and insert PET df1$pet <- thornthwaite(df.t2m, lat) #Check head(df1) #Compute and insert CWB df1$cwb <- df1$era.latlon.mm-df1$pet #Check head(df1) ##Make a seperate cwb dataframe for the SPEI trend analysis df.cwb <- as.data.frame(df1$cwb) ##Convert to a time series object to sanity check #cwb.ts1 <- ts(df1$cwb, end=c(2018,12), frequency=12) #plot(cwb.ts) #cwb.ts2 <- ts(df.cwb, end=c(2018,12), frequency=12) #plot(cwb.ts2) ## Compute 3-months SPEI spei_3 <- spei(df.cwb, 3) spei3_coeff <- as.numeric(spei_3$fitted) ##DROUGHT spei3_drought <- spei3_coeff spei3_drought[spei3_drought > -1.50 ] <- NA ##Set up the final dataframe with month and extracted value print("Making final dataframe from which to compute percentiles and slopes") datafinal.drought.spei3 <- data.frame(dates= eradatadates, obs = spei3_drought) ##Needed when you count the number of years in the regression ymddates <- as.Date(datafinal.drought.spei3$dates) #Count the number of "drought" events spei3_droughtcount <- datafinal.drought.spei3 spei3_droughtcount$obs[spei3_droughtcount$obs <= 10 ] <- 1 #10 covers wet and dry ##Annual count of spei years <- as.Date(ymddates,"%Y-%m-%d") annualcount.spei3_drought <- rowsum(spei3_droughtcount,na.rm=TRUE,format(ymddates,'%Y')) ##BOOTSTRAPPING to test trend significance print("BOOTSTRAPPING") print("USING 3-MONTH SPEI") #A function that will generate slope values for each bootstrap run yrs <- c(1:length(annualcount.spei3_drought$obs)) annualcount.spei3_dobs <- annualcount.spei3_drought$obs slope.func <- function(annualcount.spei3_drought,yrs) { slope.test <- lm(annualcount.spei3_dobs ~yrs) slope.test$coefficients[[2]] } ##Run the bootstrap - "R" is the number of bootstrapping runs testboot <- boot(annualcount.spei3_dobs,slope.func,R=1000,sim="ordinary",stype="i") ##Finding 95 percentile value of the bootstrap regressions - if your real regression value falls above the 95 percentile of all the regression values the trend is significant ex95 <- quantile(testboot$t,0.95) #$t are the boostapped values, $t0 is real data value ex5 <- quantile(testboot$t,0.05) #$t are the boostapped values, $t0 is real data value ##Print the 95 percentile value of the bootstrap regressions and the real value of the regression to compare. If your real regression value falls above (below) the 95 (5th) percentile of all the regression values the trend is significant. print(c(testboot$t0, ex5[[1]], ex95[[1]])) ##Assigning observed change where 0=no change (slope is 0); +/-1 is non-significant change; +/-2 is significant change print("Allocating out change slope") #slope is 0, doing a range between -0.1 to 0.1 to account for long decimals if((testboot$t0 <= 0.05)&(testboot$t0 >= -0.05)) { change = "0" } #Positive slope if((testboot$t0 > 0.05)&(testboot$t0 < ex95[[1]])) { change = "1" } if((testboot$t0 >= ex95[[1]])) { change = "2" } #Negative slope if((testboot$t0 < -0.05)&(testboot$t0 > ex5[[1]])) { change = "-1" } if((testboot$t0 <= ex5[[1]])) { change = "-2" } #print(change) print(paste(id,lat,lon,change,sep=",")) ################## #Write change signal to the table print("Writing to table") write.table((paste(id,lat,lon,change,testboot$t0,sep=",")),file=paste(outdir,outfile1,sep="/"),col.names=FALSE,quote=FALSE,sep="",row.names=FALSE,append=TRUE) print("MOVING TO NEXT ROW") print("") } print("DONE!")