Commit 0915c0a9 authored by Sebastian Milinski's avatar Sebastian Milinski
Browse files

Added plotting script (for GMST)

parent e38e0c02
;*************************************************
; plotting tsurf time series
;
; Reads time series data from file (no 2D input allowed)
;
; internal calculations
; - seasonal means
; - annual mean
; - subtracting 1960-1990 average from all timeseries
;
;*************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" ; These four libraries are automatically
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" ; loaded from NCL V6.4.0 onward.
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" ; No need for user to explicitly load.
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
begin
;*************************************************
;*** settings ***
plotseason = 4 ; ( DJF | MAM | JJA | SON | yearmean )
if (all(isdefined((/"varname","optionstring"/)))) then
print("using production mode.")
devmode = 0
plottype = "pdf"
else
print("using devmode")
devmode = 1
plottype = "X11"
end if
if (devmode .eq. 1) then
varname = "tsurf"
optionstring = "global.mean"
end if
experiments = (/"hist","onepct","rcp26","rcp45","rcp85"/)
seasons = (/"DJF","MAM","JJA","SON","yearmean"/)
basedir = "/work/mh0033/m300265/processed_thesis"
plotdir = "/home/mpim/m300265/thesis/plots/ch1"
plotfilename = varname + "." + optionstring
plotfilepath = plotdir + "/" + plotfilename
system("if ! test -d " + plotdir +" ; then mkdir " + plotdir + " ; fi")
ylim_override = 0
; ymin = 0.4
; ymax = 1.3
print("selected options:")
print("Variable: " + varname)
print("options: " + optionstring)
print(" ")
print("loading data...")
tStrt = get_cpu_time()
; creating:
; data_hist_seas [ season | year | ens ]
; data_onepct_seas
; data_rcp26_seas
; data_rcp45_seas
; data_rcp85_seas
do expno = 0, dimsizes(experiments)-1
experiment = experiments(expno)
if (experiment .eq. "hist") then
startyear = 1850
endyear = 2005
r0 = 1
rn = 100
ens_string = experiment + "_" + startyear + "-" + endyear + "_ens_" + r0 + "-" + rn
datdir = basedir + "/" + experiment
ifile = datdir + "/" + ens_string + "." + varname + "." + optionstring + ".nc"
fdat = addfile(ifile,"r")
data_raw = fdat->$varname$(:,:,0,0)
data_raw&ens = ispan(r0,rn,1)
data_hist_seas = new((/dimsizes(seasons),dimsizes(data_raw&time)/12,dimsizes(data_raw&ens)/),typeof(data_raw))
data_hist_seas!0 = "season"
data_hist_seas&season = seasons
data_hist_seas!1 = "year"
data_hist_seas&year = ispan(startyear,endyear,1)
data_hist_seas!2 = "ens"
data_hist_seas&ens = ispan(r0,rn,1)
do iens = 0, dimsizes(data_hist_seas&ens)-1
data_hist_seas(0:3,:,iens) = (/month_to_seasonN(data_raw(:,iens),seasons(0:3))/)
data_hist_seas(4,:,iens) = (/month_to_annual(data_raw(:,iens), 1)/)
end do
delete(fdat)
delete(data_raw)
else if (experiment .eq. "onepct") then
startyear = 1850
endyear = 1999
r0 = 1
rn = 100
ens_string = experiment + "_" + startyear + "-" + endyear + "_ens_" + r0 + "-" + rn
datdir = basedir + "/" + experiment
ifile = datdir + "/" + ens_string + "." + varname + "." + optionstring + ".nc"
fdat = addfile(ifile,"r")
data_raw = fdat->$varname$(:,:,0,0)
data_raw&ens = ispan(r0,rn,1)
data_onepct_seas = new((/dimsizes(seasons),dimsizes(data_raw&time)/12,dimsizes(data_raw&ens)/),typeof(data_raw))
data_onepct_seas!0 = "season"
data_onepct_seas&season = seasons
data_onepct_seas!1 = "year"
data_onepct_seas&year = ispan(startyear,endyear,1)
data_onepct_seas!2 = "ens"
data_onepct_seas&ens = ispan(r0,rn,1)
do iens = 0, dimsizes(data_onepct_seas&ens)-1
data_onepct_seas(0:3,:,iens) = (/month_to_seasonN(data_raw(:,iens),seasons(0:3))/)
data_onepct_seas(4,:,iens) = (/month_to_annual(data_raw(:,iens), 1)/)
end do
delete(fdat)
delete(data_raw)
else if (experiment .eq. "rcp26") then
startyear = 2006
endyear = 2099
r0 = 1
rn = 100
ens_string = experiment + "_" + startyear + "-" + endyear + "_ens_" + r0 + "-" + rn
datdir = basedir + "/" + experiment
ifile = datdir + "/" + ens_string + "." + varname + "." + optionstring + ".nc"
fdat = addfile(ifile,"r")
data_raw = fdat->$varname$(:,:,0,0)
data_raw&ens = ispan(r0,rn,1)
data_rcp26_seas = new((/dimsizes(seasons),dimsizes(data_raw&time)/12,dimsizes(data_raw&ens)/),typeof(data_raw))
data_rcp26_seas!0 = "season"
data_rcp26_seas&season = seasons
data_rcp26_seas!1 = "year"
data_rcp26_seas&year = ispan(startyear,endyear,1)
data_rcp26_seas!2 = "ens"
data_rcp26_seas&ens = ispan(r0,rn,1)
do iens = 0, dimsizes(data_rcp26_seas&ens)-1
data_rcp26_seas(0:3,:,iens) = (/month_to_seasonN(data_raw(:,iens),seasons(0:3))/)
data_rcp26_seas(4,:,iens) = (/month_to_annual(data_raw(:,iens), 1)/)
end do
delete(fdat)
delete(data_raw)
else if (experiment .eq. "rcp45") then
startyear = 2006
endyear = 2099
r0 = 1
rn = 100
ens_string = experiment + "_" + startyear + "-" + endyear + "_ens_" + r0 + "-" + rn
datdir = basedir + "/" + experiment
ifile = datdir + "/" + ens_string + "." + varname + "." + optionstring + ".nc"
fdat = addfile(ifile,"r")
data_raw = fdat->$varname$(:,:,0,0)
data_raw&ens = ispan(r0,rn,1)
data_rcp45_seas = new((/dimsizes(seasons),dimsizes(data_raw&time)/12,dimsizes(data_raw&ens)/),typeof(data_raw))
data_rcp45_seas!0 = "season"
data_rcp45_seas&season = seasons
data_rcp45_seas!1 = "year"
data_rcp45_seas&year = ispan(startyear,endyear,1)
data_rcp45_seas!2 = "ens"
data_rcp45_seas&ens = ispan(r0,rn,1)
do iens = 0, dimsizes(data_rcp45_seas&ens)-1
data_rcp45_seas(0:3,:,iens) = (/month_to_seasonN(data_raw(:,iens),seasons(0:3))/)
data_rcp45_seas(4,:,iens) = (/month_to_annual(data_raw(:,iens), 1)/)
end do
delete(fdat)
delete(data_raw)
else if (experiment .eq. "rcp85") then
startyear = 2006
endyear = 2099
r0 = 1
rn = 30
ens_string = experiment + "_" + startyear + "-" + endyear + "_ens_" + r0 + "-" + rn
datdir = basedir + "/" + experiment
ifile = datdir + "/" + ens_string + "." + varname + "." + optionstring + ".nc"
fdat = addfile(ifile,"r")
data_raw = fdat->$varname$(:,:,0,0)
data_raw&ens = ispan(r0,rn,1)
data_rcp85_seas = new((/dimsizes(seasons),dimsizes(data_raw&time)/12,dimsizes(data_raw&ens)/),typeof(data_raw))
data_rcp85_seas!0 = "season"
data_rcp85_seas&season = seasons
data_rcp85_seas!1 = "year"
data_rcp85_seas&year = ispan(startyear,endyear,1)
data_rcp85_seas!2 = "ens"
data_rcp85_seas&ens = ispan(r0,rn,1)
do iens = 0, dimsizes(data_rcp85_seas&ens)-1
data_rcp85_seas(0:3,:,iens) = (/month_to_seasonN(data_raw(:,iens),seasons(0:3))/)
data_rcp85_seas(4,:,iens) = (/month_to_annual(data_raw(:,iens), 1)/)
end do
delete(fdat)
delete(data_raw)
else
print( "Error: " + experiment + " not found.")
end if
end if
end if
end if
end if
end do
obsdir = "/work/mh0033/m300265/GMST/"
obsfile = "NCEP_air2m.GMST.nc"
ifile = obsdir + obsfile
fdat = addfile(ifile,"r")
data_raw = fdat->air(:,0,0)
startyear = 1948
endyear = 2017
; data_raw&time = ispan(1948, 2017, 1)
; data_raw = data_raw-273.15
; obs_reflevel = avg(data_raw({1961:1990}))
; data_raw = data_raw-obs_reflevel
data_obs_seas = new((/dimsizes(seasons),dimsizes(data_raw&time)/12/),typeof(data_raw))
data_obs_seas!0 = "season"
data_obs_seas&season = seasons
data_obs_seas!1 = "year"
data_obs_seas&year = ispan(startyear,endyear,1)
data_obs_seas(0:3,:) = (/month_to_seasonN(data_raw,seasons(0:3))/)
data_obs_seas(4,:) = (/month_to_annual(data_raw,1)/)
delete(fdat)
delete(data_raw)
obs_reflevel = avg(data_obs_seas(4,{1961:1990}))
hist_reflevel = dim_avg_n(dim_avg_n(data_hist_seas(4,{1961:1990},:),0),0)
data_obs_seas = data_obs_seas - obs_reflevel
data_hist_seas = data_hist_seas - hist_reflevel
data_onepct_seas = data_onepct_seas - hist_reflevel
data_rcp26_seas = data_rcp26_seas - hist_reflevel
data_rcp45_seas = data_rcp45_seas - hist_reflevel
data_rcp85_seas = data_rcp85_seas - hist_reflevel
printVarSummary(data_obs_seas)
printVarSummary(data_hist_seas)
print("loading data: " + (get_cpu_time() - tStrt) + " seconds")
;*************************************************
; print("Processing...")
; tStrt = get_cpu_time()
; print("processing: " + (get_cpu_time() - tStrt) + " seconds")
;*************************************************
;*************************************************
print("Plotting...")
;filename defined in the beginning
wks = gsn_open_wks(plottype,plotfilepath)
; timeseries
res2 = True ; plot mods desired
res2@gsnDraw = False ; draw all plots when drawing panel
res2@gsnFrame = False ; don't advance frame to show all plots on the same page
res2@tiMainString = " " ; add title
res2@tiXAxisString = " ";"year"
res2@tiYAxisString = " ";"GMST rel. to 1961-1990 [K]"
res2@xyDashPattern = 0 ; Make curves all solid
res2@xyLineColor = "blue"
res2@pmLegendDisplayMode = "Never"
res2@vpHeightF = 0.3
res2@vpWidthF = 0.9
res2@trYMinF = -0.4
res2@trYMaxF = 1.0
res2@tmXBMode = "Manual"
res2@tmYRBorderOn = False
res2@tmYROn = False
res2@tmXTBorderOn = False
res2@tmXTOn = False
linecol = new((/101,4/),float)
rgba = namedcolor2rgba("lightskyblue")
linecol(:,:) = conform(linecol,rgba(0,:),1)
linecol(0,:) = namedcolor2rgba("black")
res2@xyLineColors = linecol
linewidth = new((/101/),float)
linewidth(:) = 2.0 ;1
linewidth(0) = 4.0
res2@xyLineThicknesses = linewidth
;*************************************************
; panel plot res
;*************************************************
resP = True ; modify the panel plot
resP@txString = " "
resP@gsnFrame = False
resP@gsnPanelLabelBar = False ; add common colorbar
resP@gsnMaximize = True ; maximizes plot and rotates to landscape if necessary, needs to be set to force portrait or landscape
resP@gsnPaperOrientation = "landscape"
;************************************************
; Set resources for customizing a simple legend
;************************************************
; genres = True
; genres@XPosPercent = 10 ; move to the right
; genres@ItemSpacePercent = 6
; textres = True
; textres@lgLabelFontHeights = 0.015
; textres@lgLabels = (/"ensemble stddev",polyorder + "th-order polyfit", runave_length + "-year runmean"/)
; textres@lgPerimOn = False ; no perimeter
; textres@lgItemCount = 3 ; how many
; lineres = True
; lineres@lgLineThicknesses = 3 ; line thickness
; lineres@LineLengthPercent = 8 ; expressed as %, 0->100, length of line
; lineres@lgLineLabelFontHeights = 0.015 ; font height
; lineres@lgDashIndex = 0 ; line patterns
; lineres@lgLineColors = (/"red","black","blue"/)
;*************************************************
; plot timeseries
; Obs
res2@trXMinF = floor(min(data_obs_seas&year))
res2@trXMaxF = ceil(max(data_obs_seas&year))
res2@tmXBTickStartF = res2@trXMinF
obs = gsn_csm_xy(wks,data_obs_seas&year,data_obs_seas(plotseason,:),res2) ; create plot
obs_copy = gsn_csm_xy(wks,data_obs_seas&year,data_obs_seas(plotseason,:),res2) ; create same again to reuse for overlay
gsn_panel(wks,obs,(/1,1/),resP)
frame(wks)
; Add a few hist members
linecol(0,:) = namedcolor2rgba("lightskyblue")
res2@xyLineColors = linecol
ens_hist_1 = gsn_csm_xy(wks,data_hist_seas&year,data_hist_seas(season|plotseason,ens|1,year|:),res2) ; create plot
overlay(obs_copy,ens_hist_1)
gsn_panel(wks,obs_copy,(/1,1/),resP)
frame(wks)
ens_hist_2 = gsn_csm_xy(wks,data_hist_seas&year,data_hist_seas(season|plotseason,ens|2,year|:),res2) ; create plot
overlay(obs_copy,ens_hist_2)
gsn_panel(wks,obs_copy,(/1,1/),resP)
frame(wks)
ens_hist_3 = gsn_csm_xy(wks,data_hist_seas&year,data_hist_seas(season|plotseason,ens|3,year|:),res2) ; create plot
overlay(obs_copy,ens_hist_3)
gsn_panel(wks,obs_copy,(/1,1/),resP)
frame(wks)
; hist
linewidth(0) = 2.0
res2@xyLineThicknesses = linewidth
linecol(0,:) = namedcolor2rgba("lightskyblue")
res2@xyLineColors = linecol
ens_hist = gsn_csm_xy(wks,data_hist_seas&year,data_hist_seas(season|plotseason,ens|:,year|:),res2) ; create plot
linecol(0,:) = namedcolor2rgba("blue")
res2@xyLineColors = linecol
linewidth(0) = 4.0
res2@xyLineThicknesses = linewidth
ensmean_hist = gsn_csm_xy(wks,data_hist_seas&year,dim_avg_n_Wrap(data_hist_seas(season|plotseason,ens|:,year|:),0),res2) ; create plot
overlay(ens_hist,obs)
gsn_panel(wks,ens_hist,(/1,1/),resP)
frame(wks)
overlay(ens_hist,ensmean_hist)
overlay(ens_hist,obs)
gsn_panel(wks,ens_hist,(/1,1/),resP)
frame(wks)
;*************************************************
; whole time period covered by ensemble
;*************************************************
res2@trXMinF = floor(min(data_hist_seas&year))
res2@trXMaxF = ceil(max(data_rcp26_seas&year))
res2@tmXBTickStartF = res2@trXMinF
res2@trYMinF = -1.0
res2@trYMaxF = 4.5
linewidth(0) = 2.0
res2@xyLineThicknesses = linewidth
linecol(0,:) = namedcolor2rgba("lightskyblue")
res2@xyLineColors = linecol
ens_hist = gsn_csm_xy(wks,data_hist_seas&year,data_hist_seas(season|plotseason,ens|:,year|:),res2) ; create plot
linewidth(0) = 4.0
res2@xyLineThicknesses = linewidth
linecol(0,:) = namedcolor2rgba("blue")
res2@xyLineColors = linecol
ensmean_hist = gsn_csm_xy(wks,data_hist_seas&year,dim_avg_n_Wrap(data_hist_seas(season|plotseason,ens|:,year|:),0),res2) ; create plot
linewidth(0) = 2.0
res2@xyLineThicknesses = linewidth
rgba = namedcolor2rgba("lightpink")
linecol(:,:) = conform(linecol,rgba(0,:),1)
res2@xyLineColors = linecol
ens_onepct = gsn_csm_xy(wks,data_onepct_seas&year,data_onepct_seas(season|plotseason,ens|:,year|:),res2) ; create plot
linewidth(0) = 4.0
res2@xyLineThicknesses = linewidth
linecol(0,:) = namedcolor2rgba("red")
res2@xyLineColors = linecol
ensmean_onepct = gsn_csm_xy(wks,data_onepct_seas&year,dim_avg_n_Wrap(data_onepct_seas(season|plotseason,ens|:,year|:),0),res2) ; create plot
overlay(ens_hist,ens_onepct)
linewidth(0) = 2.0
res2@xyLineThicknesses = linewidth
rgba = namedcolor2rgba("gold")
linecol(:,:) = conform(linecol,rgba(0,:),1)
res2@xyLineColors = linecol
ens_rcp26 = gsn_csm_xy(wks,data_rcp26_seas&year,data_rcp26_seas(season|plotseason,ens|:,year|:),res2) ; create plot
linewidth(0) = 4.0
res2@xyLineThicknesses = linewidth
linecol(0,:) = namedcolor2rgba("goldenrod4")
res2@xyLineColors = linecol
ensmean_rcp26 = gsn_csm_xy(wks,data_rcp26_seas&year,dim_avg_n_Wrap(data_rcp26_seas(season|plotseason,ens|:,year|:),0),res2) ; create plot
overlay(ens_hist,ens_rcp26)
linewidth(0) = 2.0
res2@xyLineThicknesses = linewidth
rgba = namedcolor2rgba("palegreen")
linecol(:,:) = conform(linecol,rgba(0,:),1)
res2@xyLineColors = linecol
ens_rcp45 = gsn_csm_xy(wks,data_rcp45_seas&year,data_rcp45_seas(season|plotseason,ens|:,year|:),res2) ; create plot
linewidth(0) = 4.0
res2@xyLineThicknesses = linewidth
linecol(0,:) = namedcolor2rgba("green4")
res2@xyLineColors = linecol
ensmean_rcp45 = gsn_csm_xy(wks,data_rcp45_seas&year,dim_avg_n_Wrap(data_rcp45_seas(season|plotseason,ens|:,year|:),0),res2) ; create plot
overlay(ens_hist,ens_rcp45)
linewidth(0) = 2.0
res2@xyLineThicknesses = linewidth
rgba = namedcolor2rgba("magenta1")
linecol(:,:) = conform(linecol,rgba(0,:),1)
res2@xyLineColors = linecol
ens_rcp85 = gsn_csm_xy(wks,data_rcp85_seas&year,data_rcp85_seas(season|plotseason,ens|:,year|:),res2) ; create plot
linewidth(0) = 4.0
res2@xyLineThicknesses = linewidth
linecol(0,:) = namedcolor2rgba("magenta4")
res2@xyLineColors = linecol
ensmean_rcp85 = gsn_csm_xy(wks,data_rcp85_seas&year,dim_avg_n_Wrap(data_rcp85_seas(season|plotseason,ens|:,year|:),0),res2) ; create plot
overlay(ens_hist,ens_rcp85)
linewidth(0) = 4.0
res2@xyLineThicknesses = linewidth
linecol(0,:) = namedcolor2rgba("black")
res2@xyLineColors = linecol
obs = gsn_csm_xy(wks,data_obs_seas&year,data_obs_seas(plotseason,:),res2) ; create plot
overlay(ens_hist,ensmean_hist)
overlay(ens_hist,ensmean_onepct)
overlay(ens_hist,ensmean_rcp26)
overlay(ens_hist,ensmean_rcp45)
overlay(ens_hist,ensmean_rcp85)
overlay(ens_hist,obs)
gsn_panel(wks,ens_hist,(/1,1/),resP)
if (devmode .eq. 1) then
drawNDCGrid(wks) ; This is for debugging purposes
end if
frame(wks)
print(plottype)
print(plotfilename)
end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment