S.Sénési for Météo-France - sept 2019 to march 2021
from __future__ import print_function
import os
version = ""
input_dir = "./changes"
outdir = "./figures"
figure_name = "fig_SOD_BoxTS.X_f3_h"
panel = "h"
title = ""
yaxis_title = "%% change in space averages "
#
data_versions_tag = "20210201_derived"
excluded_models = []
variables = [("pr","mean"),("mrro","mean"), ("pr","std"), ("mrro","std")]
option = "mean"
combined_seasons = [ "tropics_JJA", "extra_summer" ]
scenarios = ["ssp585"]
#
colors = {"tropics_JJA" : {"pr":"indianred1","mrro":"cornflowerblue"},
"extra_summer" : {"pr":"firebrick3","mrro":"blue3"},
"globe_annual" : { "pr":"firebrick3","mrro":"blue3","prw":"gray40","tas":"gray40"},
"extra_annual" : { "pr":"firebrick3","mrro":"blue3","prw":"gray40","tas":"gray40"},
"tropics_annual" : { "pr":"firebrick3","mrro":"blue3","prw":"gray40","tas":"gray40"}}
thicknesses = {"tropics_JJA" :10.0,"extra_summer" : 4.0 , "globe_annual":10.0, "tropics_annual":10.0, "extra_annual":10.0}
xyMarkLineModes = {"mean":"MarkLines","std":"MarkLines"}
xyDashPatterns = {"mean":0,"std":1 }
xyMarkerSizes = {"pr_mean" : 0.015,"mrro_mean" : 0.015,"pr_std" : 0.01,"mrro_std": 0.01, "prw_mean": 0.01,"tas_mean":0.015, "tas_std":0.01}
nice = {"pr":"precipitation","prw":"precipitable water","mrro":"runoff",
"tropics_JJA":"Tropical land JJA",
"extra_summer" : "Extra-tropical summer ~C~(NH:JJA, SH:DJF)",
"globe_annual" : "Globe annual", "tropics_annual" : "Tropics annual",
"extra_annual" : "Extra-tropics annual",
"mean":"annual mean", "std" : "inter-annual variability",
"ssp126":"SSP1-1.9","ssp126":"SSP1-2.6","ssp245":"SSP2-4.5","ssp585":"SSP5-8.5",
}
only_warmer_CI = True # Should we show Confidence Interval bars only at warm end
show_variability_CI = True # Should we show Confidence Interval bars for variability
show_mean_CI = True # Should we show Confidence Interval bars for means
show_tas_CI = False # Should we show Confidence Interval bars for GSAT
plot_intermediate_CIs = False # Should we show Confidence Interval bars for all warming levels
xy_ranges = (0.,6.5,-15.,60.)
#
max_warming_level = 5.
min_models_nb = 7
do_test = True
if do_test :
input_dir = "/data/ssenesi/prod_20210201/FigBoxTS.X_f3_h/changes"
version = ""
variables = [("mrro","mean")]
variables = [("pr","mean"),("mrro","mean"), ("pr","std"), ("mrro","std"),("prw","mean")]
combined_seasons = [ "tropics_annual" ]
title = "Hydrological variables change over tropical land"
xy_ranges = (1.,6.5,-10.,50.)
yaxis_title = "% change "
do_test_extra = False
if do_test_extra :
input_dir="/data/ssenesi/prod_20210201/FigBoxTS.X_f3_h_extra/changes"
version="extra"
variables = [("mrro","mean")]
variables = [("pr","mean"),("mrro","mean"), ("pr","std"), ("mrro","std"),("prw","mean")]
combined_seasons = [ "extra_annual" ]
title = "Hydrological variables change over extra-tropical land"
xy_ranges = (1.,6.5,0.,60.)
yaxis_title = "% change "
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
import numpy as np
import numpy.ma as ma
import Ngl
import os
import json
import copy
from IPython.display import Image
def ensemble_stat_series(stats,variable,stat,season,option) :
""" Compute mean or quantile on an ensemble """
#
global excluded_models, scenarios
series=stats[variable][stat][season]["ens"]
periods=series.keys()
periods.sort()
ret=[]
for p in periods :
if compute_option != "parametric" and float(p)> max_warming_level :
#print("skipping %g",float(p))
ret.append(missing_value)
continue
ens=series[p]
#if len(ens) < min_models_nb:
# print("skipping %g",float(p))
# ret.append(missing_value)
# continue
if option == "x":
ret.append(np.float64(p)) ############################################################### a revoir pour le cas parametric
continue
####
ens=ens.values()
if option == "mean":
ret.append(np.mean(ens))
elif option == "median":
ret.append(np.median(ens))
#
elif option in ["min","max"] :
l=[ value for value in ens ]
l.sort()
if option=="min" : ret.append(np.float64(l[0]))
if option=="max" : ret.append(np.float64(l[-1]))
#
elif option in ["nq5","nq95","nq17","nq83" ] : #percentiles with gaussian hypothesis
mean=np.mean(ens)
std1=np.std(ens,ddof=1)
if option == "nq5" : ret.append(mean - 1.645*std1)
if option == "nq95" : ret.append(mean + 1.645*std1)
if option == "nq17" : ret.append(mean - 0.954*std1)
if option == "nq83" : ret.append(mean + 0.954*std1)
#
else:
raise ValueError("Time to use some library for stats (%s) !"%option)
#print( "\ntime series for %s %s %s %s"%(variable,stat,season,option), ret)
return(ret)
def extract_last_defined_values(all_error_bars) :
for x,center,up,low in all_error_bars :
if x != missing_value and center != missing_value :
rep=(x,center,up,low)
return rep
def season_plot(wks,stats,season,variables,option,errors,only_warmer=only_warmer_CI,
show_variability_CI=show_variability_CI,show_tas_CI=show_tas_CI,
show_mean_CI=show_mean_CI, plot_all_CIs=plot_intermediate_CIs,
rank=1,key="",xy_ranges=(0.,6.5,-15.,60.),title="",models=None):
"""
rank : allows to decide for a shift of the legend
key : text for the legend
"""
#GTAS_shift=0.89 # Shift for GTAS between 1995-2014 and 1850-1900, based on AR6 WL Box table
#GTAS_shift=0. # No shift if GTAS computed for 1850-1900
#
# Gathering time series of ensemble mean of time statistics for required variables
means=[] # their mean value
x=[] # their x-axis, which is a warming level if compute_option != parametric
for variable,stat in variables :
#d=stats[variable][stat][season][option]
#periods=sorted(d)
#means.append([ d[p] for p in periods ])
print(variable,stat)
means.append(ensemble_stat_series(stats,variable,stat,season,option))
x.append (ensemble_stat_series(stats,variable,stat,season,"x" ))
#
if len(means)==0 :
raise ValueError("No matching variable in stats "+repr(stats.keys()))
#print()"means=",means)
if len(means)>0 :
y=np.array(means,np.float32)
y=ma.masked_values(y,missing_value)
if compute_option == "parametric" :
periods = stats["tas"]["mean"]["globe_year"]["mean"].keys()
periods.sort()
GTAS = [ stats["tas"]["mean"]["globe_year"]["mean"][p] for p in periods ]
else:
GTAS=np.array(x)
#-- set resources
res = Ngl.Resources() #-- generate an res object for plot
res.nglDraw = False # don't draw
res.nglFrame = False # don't advance frame
res.gsnMaximize = False
res.gsnPaperOrientation="landscape"
res.tiMainString = title
res.tiMainFontHeightF = 0.02
#
res.trXMinF,res.trXMaxF,res.trYMinF,res.trYMaxF = xy_ranges
#
res.tmYMajorGrid = True
dpat = Ngl.new_dash_pattern(wks,"_____$_____$_____$_____$_____$_____$_____$_____$_____$_____")
res.tmYMajorGridLineDashPattern = dpat #15
#res.tmYLMode = "Explicit"
res.tmYLValues = [-10.,0.,10.,20.,30.,40.,50.]
#res.tmYLMode = "Automatic"
res.tmYLMode = "Manual"
res.tmYLTickStartF = -20.
res.tmYLDataBottomF = -20.
res.tmYLDataTopF = 40.
res.tmYLTickSpacingF= 20
#
res.tmXBMode = "Explicit"
res.tmXBValues = [ 1.5,2,2.5,3,3.5,4,4.5,5]
res.tmXBLabels = [ "","2","","3","","4","","5"]
res.tmXBTickStartF = 0.
res.tmXBTickSpacingF= 1.
res.tmXBDataLeftF = 1.4
res.tmXBDataRightF = 5.4
#
scenario_string=""
if scenarios is not None :
scenario_string=" -"
for s in scenarios :
scenario_string += " "+nice[s]+","
scenario_string=scenario_string[0:-1]
res.tiYAxisString = yaxis_title + " - multi-model mean" + scenario_string
if models is not None :
res.tiYAxisString = yaxis_title + " - %d models mean "%len(models) + scenario_string
res.tiYAxisFontHeightF = 0.016
res.tiXAxisString = "GSAT change since pre-industrial (degrees)" #-- x-axis title
res.tiXAxisFontHeightF = 0.016
#-- xy-plot resources
res.xyLineColors = [colors[season][v] for v,s in variables] #-- set 3 different colors for lines
res.xyLineThicknessF = [thicknesses[season]] #-- line thickness for all
res.xyDashPatterns = [ xyDashPatterns [s] for v,s in variables]
#
res.xyMarkLineModes = [ xyMarkLineModes[s] for v,s in variables]
res.xyMarker = 16
res.xyMarkerSizes = [ xyMarkerSizes.get("%s_%s"%(v,s),0.01) for v,s in variables]
res.xyMarkerColors = [colors[season][v] for v,s in variables]
#-- legend resources
res.lgAutoManage = False
res.lgLabelJust = "CenterLeft"
res.lgBoxBackground = 0
res.pmLegendDisplayMode = "Always"
res.lgLabelFontAspectF = 1.5
#res.lgLabelConstantSpacingF =0.0
if rank == 1 :
res.xyExplicitLegendLabels = [" "+ nice[v]+" "+nice[s] for v,s in variables ] #-- set explicit legend labels
res.pmLegendParallelPosF = -0.45
else :
res.xyExplicitLegendLabels = [ "" for v,s in variables ] #-- no legend labels
res.pmLegendParallelPosF = 0.15
#
res.pmLegendOrthogonalPosF = 0.4 #-- move the legend upwards
res.pmLegendZone = 0 #-- legend zone: 0 = topLeft; 6 = topRight
res.lgJustification = "TopLeft" #-- legend justification
res.pmLegendWidthF = 0.15 #-- change width
res.pmLegendHeightF = 0.1
res.lgLabelFontHeightF = 0.015
res.pmLegendSide = "Top" #-- Change location
res.lgPerimOn = False #-- turn off the perimeter
res.xyLineLabelFontColors = [colors[season][v] for v,s in variables]
#
if key is None :
res.pmLegendOrthogonalPosF = 0.47 #1.6 #-- move the legend upwards
#-- creates the plot
plot=Ngl.xy(wks,GTAS,y,res)
#
# Plotting model numbers for each (or the first) variable
#
if models is None :
count=0
deltay=3
processed=[]
#print("variables=",variables)
for var,stat in variables :
if var in processed : continue
processed.append(var)
if compute_option != "parametric" :
resn= Ngl.Resources()
resn.txFontHeightF = 0.017
resn.txJust ="TopRight"
resn.txFontColor = colors[season][var]
y_model_numbers = res.trYMinF + 6 + count*deltay
#print "y_model_numbers=",y_model_numbers
Ngl.add_text(wks, plot, "#", 2. - 0.3, y_model_numbers, resn)
resn.txJust ="TopCenter"
#
ensembles=stats[var][stat][season]["ens"]
periods=list(ensembles.keys())
periods.sort()
for p in periods :
mn=len(set(ensembles[p]))
#if mn >= min_models_nb :
#print(var,stat,p,mn,min_models_nb)
Ngl.add_text(wks, plot, "%d"%mn, p, y_model_numbers, resn)
count += 1
#
# Key label
############
#
if key is not None :
rest= Ngl.Resources()
rest.txFontHeightF = 0.02
rest.txJust ="TopCenter"
rest.txFuncCode ="~"
if rank==1:
xt=1.2
elif rank== 2 :
xt=5.2
Ngl.add_text(wks, plot, key, xt, res.trYMaxF - 3., rest)
#
#
# Handling error bars
#####################
#
#
CI_vars=variables
if show_mean_CI is False :
CI_vars=[ (var,stat) for (var,stat) in CI_vars if stat != "mean" ]
if show_variability_CI is False :
CI_vars=[ (var,stat) for (var,stat) in CI_vars if stat != "std" ]
# CI_vars=[ (var,stat) for (var,stat) in CI_vars if ( var != "prw" or stat != "std" ) ]
if show_tas_CI and compute_option == "parametric" :
CI_vars = CI_vars + [("tas","mean")]
plot_CI(stats,season,wks,plot,rank,xy_ranges,CI_vars,errors,plot_intermediate_CIs,only_warmer)
return plot
def plot_CI(stats,season,wks,plot,rank,xy_ranges,CI_vars,errors,plot_intermediate_CIs,only_warmer):
def dic2list(stat):
l=ensemble_stat_series(stats,v,s,season,stat)
return np.array(l,np.float32)
#
delta_scale=2.
for v,s in CI_vars :
#index+=1
if errors in ["second","butlast"] :
ups =dic2list("butlast")
lows =dic2list("second")
elif errors in ["nq5", "nq95"] :
ups =dic2list("nq95")
lows =dic2list("nq5")
elif errors in ["max", "min"] :
ups =dic2list("max")
lows =dic2list("min")
elif errors == "means" and "wmean" in stats[v][season] :
ups =dic2list("wmean")
lows =dic2list("iwmean")
else : raise ValueError("errors option '%s' not (yet?) available"%errors)
#
centers=dic2list("mean")
#
polyres = Ngl.Resources()
if s=="std" :
polyres.gsMarkerIndex = 1
else :
#polyres.gsMarkerIndex = 2
polyres.gsMarkerIndex = 1
#polyres.gsMarkerSizeF = .06
polyres.gsMarkerSizeF = .10
polyres.gsMarkerColor = colors[season][v]
polyres.gsLineColor = colors[season][v]
polyres.gsLineThicknessF = 20. #2.
patt=xyDashPatterns[s]
if patt==1 : patt=16
polyres.gsLineDashPattern = patt
#polyres.gsLineDashPattern = xyDashPatterns["mean"]
#
# Set x-wise shift, dependent on variable+stat
if v=="pr" :
if s != "std":
delta=0.
else :
delta=0.25
elif v=="mrro" :
if s!="std" :
delta=0.07
else :
delta=0.32
elif v=="prw" :
delta=-0.07
else : delta=0.
#
x=ensemble_stat_series(stats,v,s,season,"x" )
all_error_bars=zip(x,centers,ups,lows) ####################################### etait GTAS
if plot_intermediate_CIs :
toplot=all_error_bars
first=None
else :
first_CI=all_error_bars[0]
last_CI=extract_last_defined_values(all_error_bars)
toplot= [first_CI,last_CI]
first=True
if only_warmer :
toplot=[ last_CI ]
first=False
for t,center,up,low in toplot :
if t > 1.e+10 : continue
if toplot != all_error_bars :
if first : t0=xy_ranges[0]+0.1
elif first is False: t0=xy_ranges[1]-1.5
else :
t0=t
if rank==1 : t0=t0+0.5
if v != "tas" :
Ngl.add_polymarker(wks,plot,t0+delta*delta_scale,up,polyres)
Ngl.add_polymarker(wks,plot,t0+delta*delta_scale,center,polyres)
Ngl.add_polymarker(wks,plot,t0+delta*delta_scale,low,polyres)
Ngl.add_polyline (wks,plot,[t0+delta*delta_scale,t0+delta*delta_scale],[up,low],polyres)
if v=="prw" :
x_prw=t0+delta*delta_scale
y_prw=center
elif "prw" in [ var for var,stat in CI_vars ] :
Ngl.add_polymarker(wks,plot,x_prw+low-center,y_prw,polyres)
Ngl.add_polymarker(wks,plot,x_prw, y_prw,polyres)
Ngl.add_polymarker(wks,plot,x_prw+up-center, y_prw,polyres)
Ngl.add_polyline (wks,plot,[x_prw+up-center,x_prw+low-center],[y_prw,y_prw],polyres)
if first : first=False
#
#if index==0 and rank==1 and first != None :
if rank==1 and first != None :
rest= Ngl.Resources()
rest.txFontHeightF = 0.015
rest.txJust ="Right"
Ngl.add_text(wks, plot, "[5%-95%]~C~range" , xy_ranges[1]-1.4, xy_ranges[2]+5., rest)
def combined_figure(stats,outname,variables,option,errors,models):
#
wks_type = "png" #-- output type of workstation
wks = Ngl.open_wks(wks_type,outname)
#
# If actually combining seasons, the key must show each season name
# otherwise, we assume that the title includes it
if len(combined_seasons)==1 :
key=None
else :
key=nice[combined_seasons[0]]
plot1=season_plot(wks,stats,combined_seasons[0],variables,option,
errors,only_warmer_CI,show_variability_CI,rank=1,key=key,
xy_ranges=xy_ranges,title=title,models=models)
#return "debug"
if len(combined_seasons) == 2 :
key=nice[combined_seasons[1]]
plot2=season_plot(wks,stats,combined_seasons[1], variables,option,
errors,only_warmer_CI,show_variability_CI,rank=2, key=key,
xy_ranges=xy_ranges,models=models)
Ngl.overlay(plot1,plot2)
Ngl.draw(plot1)
#
Ngl.frame(wks) # flush graphics to the wks
Ngl.delete_wks(wks) #-- this need to be done to cope with opened wks max number
return outname+"."+wks_type
def filter_on_scenarios_models_warming(stats,scenarios,excluded_models,warming):
"""
In dict stats, keep only those results for a list of scenarios, not produced
by some models, and reaching a given warming level
Also returns the list of models reaching the level for all variables
"""
# First idenitify for each variable which models reach the warming level
# ans store it as a set in next dict
models_by_variable=dict()
for variable in stats :
if variable in [ 'option','ref_period']:continue
if 'mean' not in stats[variable] :
#print("no mean for %s"%variable)
continue
for season in stats[variable]['mean'] :
series=stats[variable]['mean'][season]["ens"]
#
# Get list of models_scenarios reaching the warming level
periods=list(series.keys())
periods.sort()
models=None
for p in periods :
if float(p) >= warming and models is None :
all_mod_scens=list(series[p].keys())
models=set([ mod_scen.split("_")[0] for mod_scen in all_mod_scens ])
if models is None : models=set()
models_by_variable[variable]=models
#print( "\n",variable,stat,season,"%d models reaching %g : "%(len(models),warming) , sorted(list(models)))
#
# Intersects those sets across variables
models=None
for variable in list(models_by_variable.keys()):
if models is None :
models=models_by_variable[variable]
else :
models=models.intersection(models_by_variable[variable])
models=list(models)
# Withdraw excluded wished models
models = [ m for m in models if m not in excluded_models ]
models.sort()
#print( "\n","all_variables",stat,season,"%d models reaching %g : "%(len(models),warming) , models)
# filter the list of models in each 'leaf' entry of dict stats
# i.e. each period in stats[variable][stat][season]["ens"]
for variable in stats :
if variable in [ 'option','ref_period']:continue
for stat in stats[variable] :
if stat not in ['mean',"std"]: continue
for season in stats[variable][stat] :
print("filtering on",variable,stat,season)
series=stats[variable][stat][season]["ens"]
#
# Get list of models_scenarios reaching the warming level
periods=list(series.keys())
periods.sort()
for p in periods :
values_dict=series[p]
filtered_dict=dict()
keys=list(values_dict.keys())
for mod_scen in keys :
mod=mod_scen.split("_")[0]
if mod in models :
# Keep only ensemble entries which are relevant for the scenarios list
if scenarios is None :
filtered_dict[mod_scen]=values_dict[mod_scen]
else:
for scenario in scenarios :
if "_"+scenario in mod_scen :
filtered_dict[mod_scen]=values_dict[mod_scen]
series[p]=filtered_dict
return models
fn="%s/stats_CMIP6_%s%s.json"%(input_dir,data_versions_tag,version)
print(fn)
with open (fn,"r") as f :
stats=json.load(f)
#print(stats.keys())
ref_period=stats["ref_period"]
compute_option=stats.get("option","parametric").encode('ascii')
missing_value=1.e+20
models=None
#filter_on_scenario_and_models(stats,scenarios,excluded_models)
#print stats
models=filter_on_scenarios_models_warming(stats,scenarios,excluded_models,max_warming_level)
print(len(models),models)
#stats=stats['ssp585']
with open("t.json","w") as f :
json.dump(stats,f,separators=(',', ': '),indent=1,ensure_ascii=True)
filen="%s/%s_%s_%s%s"%(outdir,figure_name,compute_option,data_versions_tag,version)
! mkdir -p {outdir}
figfile=combined_figure(stats,filen,variables,option, errors="nq95",models=models)
print("Image available as %s"%figfile)
Image(figfile)
all_metadata=""
for variable in stats.keys():
if variable not in [ var for var,stat in variables ] : continue
metadata=stats[variable]["metadata"]
for scenario in metadata :
for model in metadata[scenario] :
if model in models and model not in excluded_models :
string = metadata[scenario][model]
string.replace("\n"," %s\n"%panel)
all_metadata += string
with open("%s_md.txt"%(filen),"w") as f: f.write(all_metadata)