#several major imporvements on all fronts
Valentin Heimhuber
7 years ago
parent
7302cf843c
commit
9eda19058e
@ -0,0 +1,227 @@
|
|||||||
|
# -*- coding: utf-8 -*-
|
||||||
|
"""
|
||||||
|
Created on Thu Jun 28 12:15:17 2018
|
||||||
|
|
||||||
|
@author: z5025317
|
||||||
|
"""
|
||||||
|
|
||||||
|
# -*- coding: utf-8 -*-
|
||||||
|
#==========================================================#
|
||||||
|
#Last Updated - March 2018
|
||||||
|
#@author: z5025317 Valentin Heimhuber
|
||||||
|
#code for creating future climate variability deviation plots for NARCLIM variables.
|
||||||
|
#This code is derived from P1_NARCliM_First_Pass_variab_deviation_plots.py to deal with variables for which we don't have NARCLIM deltas.
|
||||||
|
#These variables are Sea Level, Acidity and...
|
||||||
|
# it uses CSIRO CC in Australia CMIP-5 Deltas instead
|
||||||
|
|
||||||
|
|
||||||
|
#==========================================================#
|
||||||
|
#Load packages
|
||||||
|
#==========================================================#
|
||||||
|
import numpy as np
|
||||||
|
import os
|
||||||
|
import pandas as pd
|
||||||
|
import glob
|
||||||
|
import matplotlib
|
||||||
|
import matplotlib.pyplot as plt
|
||||||
|
from datetime import datetime
|
||||||
|
from datetime import timedelta
|
||||||
|
from matplotlib.backends.backend_pdf import PdfPages
|
||||||
|
from ggplot import *
|
||||||
|
import csv
|
||||||
|
matplotlib.style.use('ggplot')
|
||||||
|
# Load my own functions
|
||||||
|
os.chdir('C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/Paper#1/Analysis/Code')
|
||||||
|
import climdata_fcts as fct
|
||||||
|
import silo as sil
|
||||||
|
import re
|
||||||
|
#==========================================================#
|
||||||
|
|
||||||
|
#==========================================================#
|
||||||
|
# Set working direcotry (where postprocessed NARClIM data is located)
|
||||||
|
os.chdir('C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/Paper#1/')
|
||||||
|
#==========================================================#
|
||||||
|
|
||||||
|
|
||||||
|
#==========================================================#
|
||||||
|
#set input parameters
|
||||||
|
#==========================================================#
|
||||||
|
Case_Study_Name = 'CASESTUDY2'
|
||||||
|
Casestudy2_csv_path = "Data/NARCLIM_Site_CSVs/CASESTUDY2/CASESTDUY2_NARCLIM_Point_Sites.csv"
|
||||||
|
|
||||||
|
Estuary = 'HUNTER' # 'Belongil'
|
||||||
|
Base_period_start = '1970-01-01' #Start of interval for base period of climate variability
|
||||||
|
Base_period_DELTA_start = '1990-01-01' #Start of interval used for calculating mean and SD for base period (Delta base period)
|
||||||
|
Base_period_end = '2009-01-01' #use last day that's not included in period as < is used for subsetting
|
||||||
|
Clim_var_type = 'Acidity' #Name of climate variable in NARCLIM models '*' will create pdf for all variables in folder
|
||||||
|
[minplotDelta, maxplotDelta]=[0,1]
|
||||||
|
|
||||||
|
#Source for present day climate data (historic time series) can be either: 'Station' or 'SILO'
|
||||||
|
Location = 'Estuary' # pick locaiton for extracting the SILO data can be: Estuary, Catchment, or Ocean
|
||||||
|
|
||||||
|
presentdaybar = False #include a bar for present day variability in the plots?
|
||||||
|
present_day_plot = 'yes' #save a time series of present day
|
||||||
|
Version = "V1"
|
||||||
|
Stats = 'mindaily' #'maxdaily' #maximum takes the annual max Precipitation instead of the sum
|
||||||
|
ALPHA_figs = 1 #Set alpha of figure background (0 makes everything around the plot panel transparent)
|
||||||
|
#==========================================================#
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#==========================================================#
|
||||||
|
#set directory path for output files
|
||||||
|
output_directory = 'Output/' + Case_Study_Name + '/' + Estuary + '/' + '/Clim_Deviation_Plots/'
|
||||||
|
#output_directory = 'J:/Project wrl2016032/NARCLIM_Raw_Data/Extracted'
|
||||||
|
if not os.path.exists(output_directory):
|
||||||
|
os.makedirs(output_directory)
|
||||||
|
print('-------------------------------------------')
|
||||||
|
print("output directory folder didn't exist and was generated")
|
||||||
|
print('-------------------------------------------')
|
||||||
|
print('-------------------')
|
||||||
|
#==========================================================#
|
||||||
|
|
||||||
|
|
||||||
|
#==========================================================#
|
||||||
|
#read ensemble change CSV file
|
||||||
|
#==========================================================#
|
||||||
|
|
||||||
|
clim_var_csv_path = './Data/Non_NARCLIM/AUZ_CC_CSIRO_Deltas/Brisbane_Ocean_Variables.csv'
|
||||||
|
df = pd.read_csv(clim_var_csv_path, index_col=0, parse_dates = True)
|
||||||
|
df = df.filter(regex= 'RCP 8.5|Variable')
|
||||||
|
df = df[df.Variable=='Ocean pH']
|
||||||
|
Ensemble_Delta_df = pd.DataFrame(index=['Median', '10th', '90th'], columns =['near', 'far'])
|
||||||
|
Ensemble_Delta_df['near'] = re.split('\(|\)|\ to |', df.iloc[0,1])[0:3]
|
||||||
|
Ensemble_Delta_df['far'] = re.split('\(|\)|\ to |', df.iloc[0,2])[0:3]
|
||||||
|
|
||||||
|
Ensemble_Delta_df = Ensemble_Delta_df.astype(float).fillna(0.0)
|
||||||
|
#==========================================================#
|
||||||
|
|
||||||
|
|
||||||
|
#==========================================================#
|
||||||
|
#load present day climate variable time series
|
||||||
|
#==========================================================#
|
||||||
|
#Present_day_df,minplotDelta, maxplotDelta = fct.import_present_day_climdata_csv(Estuary, Clim_var_type)
|
||||||
|
Present_day_Var_CSV = glob.glob('./Data/Non_NARCLIM/Hunter_Lower_Acidity2.csv')
|
||||||
|
Present_day_df = pd.read_csv(Present_day_Var_CSV[0],parse_dates=True, index_col=0)
|
||||||
|
Present_day_df.columns = ['PH']
|
||||||
|
df = Present_day_df
|
||||||
|
df.index = df.index.to_period('D')
|
||||||
|
Present_day_df = df
|
||||||
|
Present_day_df = df.groupby(df.index).mean().reset_index()
|
||||||
|
Present_day_df.columns = ['Date','PH']
|
||||||
|
Present_day_df.index = Present_day_df['Date']
|
||||||
|
Present_day_df= Present_day_df['PH']
|
||||||
|
#Present_day_df.index = pd.to_datetime(Present_day_df.index,format='%Y%m%d')
|
||||||
|
|
||||||
|
if Stats == 'mindaily':
|
||||||
|
Present_day_df_annual = Present_day_df.resample('A').min()
|
||||||
|
Present_day_df_annual = Present_day_df_annual.replace(0, np.nan)
|
||||||
|
if Stats == 'mean':
|
||||||
|
Present_day_df_annual = Present_day_df.resample('A').mean()
|
||||||
|
Present_day_df_annual = Present_day_df_annual.replace(0, np.nan)
|
||||||
|
|
||||||
|
Present_Day_ref_df = Present_day_df_annual
|
||||||
|
|
||||||
|
#Subset to present day variability period
|
||||||
|
Present_Day_ref_df = pd.DataFrame(Present_Day_ref_df.loc[(Present_Day_ref_df.index >= Base_period_start) & (Present_Day_ref_df.index <= Base_period_end)])
|
||||||
|
#Subset to present day variability delta base period (this is the statistical baseline used for present day conditions)
|
||||||
|
Present_Day_Delta_ref_df = pd.DataFrame(Present_Day_ref_df.loc[(Present_Day_ref_df.index >= Base_period_DELTA_start) & (Present_Day_ref_df.index <= Base_period_end)])
|
||||||
|
Present_Day_Mean = np.percentile(Present_Day_Delta_ref_df, 50)
|
||||||
|
Present_Day_SD = np.std(Present_Day_Delta_ref_df)
|
||||||
|
#create data frame for floating stacked barplots
|
||||||
|
index=['-2std', '-1std', 'Med', '1std', '2std']
|
||||||
|
columns = ['present','near future', 'far future']
|
||||||
|
Plot_in_df = pd.DataFrame(index=index, columns =columns)
|
||||||
|
#
|
||||||
|
Plot_in_df['present'] = [float(Present_Day_Mean-2*Present_Day_SD),float(Present_Day_Mean-Present_Day_SD), float(Present_Day_Mean),
|
||||||
|
float(Present_Day_Mean+Present_Day_SD), float(Present_Day_Mean+2*Present_Day_SD)]
|
||||||
|
Plot_in_df['near future'] = [float(Present_Day_Mean + Ensemble_Delta_df.near['10th']),np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.near['Median']),
|
||||||
|
np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.near['90th'])]
|
||||||
|
Plot_in_df['far future'] = [float(Present_Day_Mean + Ensemble_Delta_df.far['10th']),np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.far['Median']),
|
||||||
|
np.NaN, float(Present_Day_Mean + Ensemble_Delta_df.far['90th'])]
|
||||||
|
#Create a second data frame that has the values in a way that they can be stacked up in bars with the correct absolute values
|
||||||
|
Plot_in_df2 = pd.DataFrame(index=index, columns =columns )
|
||||||
|
#
|
||||||
|
index=['-2std', '-1std', 'Med', '1std', '2std']
|
||||||
|
columns = ['near future', 'far future']
|
||||||
|
Plot_in_df2 = pd.DataFrame(index=index, columns =columns )
|
||||||
|
|
||||||
|
Plot_in_df2['near future'] = [float(Present_Day_Mean + Ensemble_Delta_df.near['10th']),np.NaN, float(Ensemble_Delta_df.near['Median']-Ensemble_Delta_df.near['10th']),
|
||||||
|
np.NaN, float(Ensemble_Delta_df.near['90th']-Ensemble_Delta_df.near['Median'])]
|
||||||
|
Plot_in_df2['far future'] = [float(Present_Day_Mean + Ensemble_Delta_df.far['10th']),np.NaN, float(Ensemble_Delta_df.far['Median']-Ensemble_Delta_df.far['10th']),
|
||||||
|
np.NaN, float(Ensemble_Delta_df.far['90th']-Ensemble_Delta_df.far['Median'])]
|
||||||
|
#transpose the data frame
|
||||||
|
Plot_in_df_tp = Plot_in_df2.transpose()
|
||||||
|
#do the individual plots
|
||||||
|
xmin = int(min(Plot_in_df.min(axis=1))-minplotDelta)
|
||||||
|
xmax = int(max(Plot_in_df.max(axis=1))+maxplotDelta)
|
||||||
|
#define colour scheme
|
||||||
|
#likert_colors = ['none', 'firebrick','firebrick','lightcoral','lightcoral']
|
||||||
|
likert_colors = ['none', 'darkblue', 'darkblue','cornflowerblue','cornflowerblue']
|
||||||
|
uni_colors = ['none', 'cornflowerblue', 'cornflowerblue','cornflowerblue','cornflowerblue']
|
||||||
|
|
||||||
|
#plot the stacked barplot
|
||||||
|
fig = plt.figure(figsize=(14,8))
|
||||||
|
ax=plt.subplot(2,4,3)
|
||||||
|
Plot_in_df_tp.plot.bar(stacked=True, color=uni_colors, edgecolor='none', legend=False, ax=ax, width=0.5)
|
||||||
|
df = pd.DataFrame(Plot_in_df.iloc[2,[1,2]])
|
||||||
|
index2=['Med']
|
||||||
|
columns2 = ['near future', 'far future']
|
||||||
|
Plot_in_df3 = pd.DataFrame(index=index2, columns =columns2)
|
||||||
|
Plot_in_df3['near future'] = df.iloc[1,0]
|
||||||
|
Plot_in_df3['far future'] = df.iloc[0,0]
|
||||||
|
Plot_in_df3 = Plot_in_df3.transpose()
|
||||||
|
plt.plot(Plot_in_df3['Med'], linestyle="", markersize=52,
|
||||||
|
marker="_", color='darkblue', label="Median")
|
||||||
|
z = plt.axhline(float(Present_Day_Mean-2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
|
||||||
|
z.set_zorder(-1)
|
||||||
|
z = plt.axhline(float(Present_Day_Mean+2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
|
||||||
|
z.set_zorder(-1)
|
||||||
|
z = plt.axhline(float(Present_Day_Mean-Present_Day_SD), linestyle='--', color='black', alpha=.5)
|
||||||
|
z.set_zorder(-1)
|
||||||
|
z = plt.axhline(float(Present_Day_Mean+Present_Day_SD), linestyle='--', color='black', alpha=.5)
|
||||||
|
z.set_zorder(-1)
|
||||||
|
z = plt.axhline(float(Present_Day_Mean), linestyle='--', color='red', alpha=.5)
|
||||||
|
z.set_zorder(-1)
|
||||||
|
plt.ylim(xmin, xmax)
|
||||||
|
plt.title(Clim_var_type)
|
||||||
|
ax.grid(False)
|
||||||
|
for tick in ax.get_xticklabels():
|
||||||
|
tick.set_rotation(0)
|
||||||
|
fig.tight_layout()
|
||||||
|
fig.patch.set_alpha(ALPHA_figs)
|
||||||
|
|
||||||
|
#plot the present day time series
|
||||||
|
ax=plt.subplot(2,2,1)
|
||||||
|
Present_Day_ref_df.plot(legend=False, ax=ax)
|
||||||
|
z = plt.axhline(float(Present_Day_Mean-2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
|
||||||
|
z.set_zorder(-1)
|
||||||
|
z = plt.axhline(float(Present_Day_Mean+2*Present_Day_SD), linestyle='-', color='black', alpha=.5)
|
||||||
|
z.set_zorder(-1)
|
||||||
|
z = plt.axhline(float(Present_Day_Mean-Present_Day_SD), linestyle='--', color='black', alpha=.5)
|
||||||
|
z.set_zorder(-1)
|
||||||
|
z = plt.axhline(float(Present_Day_Mean+Present_Day_SD), linestyle='--', color='black', alpha=.5)
|
||||||
|
z.set_zorder(-1)
|
||||||
|
z = plt.axhline(float(Present_Day_Mean), linestyle='--', color='red', alpha=.5)
|
||||||
|
z.set_zorder(-1)
|
||||||
|
#fig.patch.set_facecolor('deepskyblue')
|
||||||
|
fig.tight_layout()
|
||||||
|
fig.patch.set_alpha(ALPHA_figs)
|
||||||
|
plt.title(Clim_var_type + ' ' + Stats + ' annual present day')
|
||||||
|
plt.ylim(xmin, xmax)
|
||||||
|
ax.grid(False)
|
||||||
|
|
||||||
|
if presentdaybar == False:
|
||||||
|
out_file_name = Estuary + '_' + Clim_var_type + '_' + Stats + '_' + Base_period_start + '_' + Base_period_end + Version + '_' + '_NPDB.png'
|
||||||
|
out_path = output_directory + '/' + out_file_name
|
||||||
|
fig.savefig(out_path)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
@ -1,90 +0,0 @@
|
|||||||
#R code for creating ggplots of time series with smooth (GAM) and linear term
|
|
||||||
|
|
||||||
|
|
||||||
######################
|
|
||||||
#Import Libraries and set working directory
|
|
||||||
######################
|
|
||||||
library(zoo)
|
|
||||||
library(hydroTSM) #you need to install these packages first before you can load them here
|
|
||||||
library(lubridate)
|
|
||||||
library(mgcv)
|
|
||||||
library(ggplot2)
|
|
||||||
library(gridExtra)
|
|
||||||
library(scales)
|
|
||||||
options(scipen=999)
|
|
||||||
setwd("C:/Users/z5025317/OneDrive - UNSW/WRL_Postdoc_Manual_Backup/WRL_Postdoc/Projects/Paper#1/")
|
|
||||||
######################
|
|
||||||
|
|
||||||
######################
|
|
||||||
#Set inputs
|
|
||||||
######################
|
|
||||||
Case.Study <- "CASESTUDY2"
|
|
||||||
Estuary <- "HUNTER"
|
|
||||||
Climvar <- 'tasmean'
|
|
||||||
ggplotGAM.k <- 7
|
|
||||||
######################
|
|
||||||
|
|
||||||
|
|
||||||
######################
|
|
||||||
#Set input file paths
|
|
||||||
######################
|
|
||||||
AirT_CSV_Path <- "./Data/Ocean_Data/BOM_monthly_SL_Hunter_Newcastle.txt"
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
dat = readLines(AirT_CSV_Path)
|
|
||||||
dat = as.data.frame(do.call(rbind, strsplit(dat, split=" {2,10}")), stringsAsFactors=FALSE)
|
|
||||||
colnames(dat) <-dat[3,]
|
|
||||||
|
|
||||||
dat2 = dat[-c(1:3), ]
|
|
||||||
dat2 = dat2[-(720:726),]
|
|
||||||
dat2 = dat2[-(1:108),]
|
|
||||||
dat2 = dat2[,-1]
|
|
||||||
|
|
||||||
dat2$Date <- as.yearmon(dat2[,1], "%m %Y")
|
|
||||||
|
|
||||||
SeaLev.df <- dat2
|
|
||||||
head(SeaLev.df)
|
|
||||||
SeaLev.df$MSL <- as.numeric(SeaLev.df$Mean)
|
|
||||||
SeaLev.df$Julday1 <- seq(1,length(SeaLev.df[,1]),1)
|
|
||||||
linear.trend.model_EC_all <- lm(MSL ~ Julday1, SeaLev.df)
|
|
||||||
SeaLev.pvalNCV_ECall <- summary(linear.trend.model_EC_all )$coefficients[2,4]
|
|
||||||
SeaLev.lintrend <- summary(linear.trend.model_EC_all )$coefficients[2,1] * 12
|
|
||||||
|
|
||||||
######################
|
|
||||||
#Plot
|
|
||||||
######################
|
|
||||||
|
|
||||||
##################################### Full Time Period
|
|
||||||
p1air <- ggplot(SeaLev.df, aes(y=MSL, x=Date)) + geom_line(alpha=0.5) +
|
|
||||||
ggtitle(paste(Estuary, " - Linear and smooth trend in monthly mean sea level (BOM Gauge) | lin trend was ",
|
|
||||||
round(100*SeaLev.lintrend,3), ' cm/year with p=', round(SeaLev.pvalNCV_ECall,10), sep=" ")) +
|
|
||||||
theme(plot.title=element_text(face="bold", size=9)) +
|
|
||||||
geom_smooth(method='lm',fill="green", formula=y~x, colour="darkgreen", size = 0.5) +
|
|
||||||
stat_smooth(method=gam, formula=y~s(x, k=4), se=T, size=0.5, col="red") +
|
|
||||||
ylab("Monthly Mean Sea Level [m]") + xlab("Time")
|
|
||||||
|
|
||||||
#export to png
|
|
||||||
png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_MonthlyMeanSeaLevel_full_period_', Sys.Date(),".png", sep="")
|
|
||||||
png(file = png.name, width = 10.5, height = 7, units='in', res=500)
|
|
||||||
grid.arrange(p1air,ncol=1)
|
|
||||||
dev.off()
|
|
||||||
|
|
||||||
#multiple smooths
|
|
||||||
p1air <- ggplot(SeaLev.df, aes(y=MSL, x=Date)) + geom_line(alpha=0.5) +
|
|
||||||
ggtitle(paste(Estuary, " - Linear and smooth trend in monthly mean sea level (BOM Gauge) | lin trend was ",
|
|
||||||
round(100*SeaLev.lintrend,3), ' cm/year with p=', round(SeaLev.pvalNCV_ECall,10), sep=" ")) +
|
|
||||||
theme(plot.title=element_text(face="bold", size=9)) +
|
|
||||||
stat_smooth(method=gam, formula=y~s(x, k=13), se=T, size=0.5, col="red") +
|
|
||||||
#stat_smooth(method=gam, formula=y~s(x, k=8), se=T, size=0.5, cor="blue") +
|
|
||||||
stat_smooth(method=gam, formula=y~s(x, k=5), se=T, size=0.5, col="green") +
|
|
||||||
ylab("Monthly Mean Sea Level [m]") + xlab("Time")
|
|
||||||
|
|
||||||
#export to png
|
|
||||||
png.name <- paste('./Output/', Case.Study, '/', Estuary, '/Trends_MonthlyMeanSeaLevel_MultiGAM_full_period_', Sys.Date(),".png", sep="")
|
|
||||||
png(file = png.name, width = 10.5, height = 7, units='in', res=500)
|
|
||||||
grid.arrange(p1air,ncol=1)
|
|
||||||
dev.off()
|
|
||||||
Loading…
Reference in New Issue