S.Sénési for Météo-France - sept 2019 to march 2021
import os
do_test = True
figure_name = "Fig8-27"
version = "" # a suffix for figure filename. Use e.g. "_V1" for legibility
confidence_factor = 1.645 # For AR6 comprehensive scheme : Multiplicative factor applied to control run
# variability for deciding a change is significant (besides sqrt(2))
sign_threshold = 0.66 # For AR6 simple scheme : threshold on cross-model change sign agreeement fraction
same_models_for_var = False
# Figure title
title = "Effect on precipitation of first versus second 2 degrees \n"+\
"of global warming (vs 1850-1900)"
outdir = "./figures"
#See doc for data_versions in sibling directory data_versions
data_versions_tag = "20201228_derived"
excluded_models = [ ]
included_models = None
variability_excluded_models= []
variability_models = None
data_versions_dir = os.getenv("CAMMAC")+"/data"
variable = "pr"
table = "Amon" # Script was yet tested only for a monthly table
variable_transformation = "plain" # Could be 'iav', 'gini', 'welsh', 'dry'...
seasons = ["DJF","JJA"] # any CDO season, not tested for "ANN". Graph is tuned for showing 2 seasons
experiment = "ssp585"
first_delta = 2.0 # Temperature change for the first interval (usually 2°)
second_delta = 4.0 # Temperature change for the ssecond interval (usually 4°)
proj_period = "2015-2099" # period investigated for the warming
ref_experiment = "historical"
ref_period = "1850-1900"
window_half_size = 10 # For time filtering of atmospheric temperature before analyzing 2K and 4K warming (unit=year)
field_type = "rmean" # Type of change field : mean or rmean (for mean of relative changes) or rmeans (for relative change fo means)
threshold = 0.1/(24*3600) # A threshold on seasonal means for individual relative changes. Can be :None. Here: 0.1 mm/day converted to kg m2 s-1
#threshold = None
# Plot tuning below is for precipitation and rmean (relative mean)
plot_args = dict(color="AR6_Precip_12",
colors="-80. -40. -20. -10. -5. 0 5. 10. 20. 40. 80. ")
with_variability = True # Should we use variability for stippling and hatching
scheme = "AR6" # Which hatching scheme ? AR5 or AR6
#
#
# Other details
figure_details = {"page_width":2450,"page_height":3444, "insert_width":2000,"pt":60, "ybox":133,"y":40}
common_grid = "r360x180"
variability_sampling_args = {"house_keeping":True,"compute":True,"detrend":True,"shift":100,"nyears":20,"number":20}
if do_test :
version = "_test"
ref_period = "1850"
included_models = ["CNRM-CM6-1"]
variability_models = ["GFDL-CM4", "CNRM-CM6-1"]#,"HadGEM3-GC31-LL","HadGEM3-GC31-MM"]
seasons = ["DJF","DJF"]
field_type = "rmeans"
import sys, os
from climaf.api import *
climaf.cache.stamping=False
from CAMMAClib.changes import change_figure, global_change
from CAMMAClib.ancillary import create_labelbar, feed_dic
from CAMMAClib.variability import agreement_fraction_on_sign, variability_AR5, stippling_hatching_masks_AR5,\
agreement_fraction_on_lower,lowchange_conflict_masks_AR6
from CAMMAClib.mips_et_al import read_versions_dictionnary, institute_for_model, mip_for_experiment, \
models_for_experiments_multi_var, models_for_experiments, TSU_metadata
# Load some user settings, if available
settings_file=os.getenv("CAMMAC_USER_PYTHON_CODE_DIR",".")+'/cammac_user_settings.py'
if os.path.exists(settings_file) :
exec(compile(open(settings_file).read(),settings_file,'exec'))
metadata=""
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
# Read dictionnary of data versions (see sibling directory data_versions)
data_versions=read_versions_dictionnary(data_versions_tag,data_versions_dir)
# Identify models with data for relevant experiments : projection and reference, for 'tas' and
# variable of interest
models=models_for_experiments_multi_var(data_versions,[(variable,table),("tas","Amon")],
[ref_experiment,experiment],excluded_models,included_models)
# Also identify models with data for computing variability
control_models=models_for_experiments(data_versions,variable,table,
["piControl"],variability_excluded_models,variability_models)
if with_variability and len(control_models)==0 :
raise ValueError("No model has data for computing variability for %s and %s (%s,%s)"%\
(variable,table,variability_excluded_models,variability_models))
# compute ensemble of warming series along projection period of choosen experiment
GSAT=global_change("tas","Amon",experiment,proj_period,ref_experiment,ref_period,models,
data_versions,filter_length=2*window_half_size+1)
max_change=dict()
models_warming_enough=[]
models_not_warming_enough=[]
year2=dict()
year4=dict()
for model,variant in models :
max_change[model]=cvalue(ccdo_fast(GSAT[model],operator="timmax"))
if max_change[model]>= 4. :
models_warming_enough.append((model,max_change[model]))
metadata+=TSU_metadata([ref_experiment,experiment],[(model,variant)],"tas","Amon",data_versions)
# year=2025 ;
year=int(proj_period.split("-")[0]) + window_half_size
found2=False
for v in cMA(GSAT[model]).flatten().data :
if v >= first_delta and not found2 :
found2=True
year2[(model,variant)]=year
if v >= second_delta :
year4[(model,variant)]=year
break
year+=1
else :
models_not_warming_enough.append((model,max_change[model]))
if with_variability :
for model,realization in control_models :
metadata+=TSU_metadata("piControl",[(model,realization)],variable,table,data_versions)
print "\nThese models don't reach %s K warming"%second_delta,models_not_warming_enough
print "\nThese %d models DO reach %s K warming"%(len(models_warming_enough),second_delta), \
models_warming_enough
print "\nYears of %s warming"%first_delta,year2
print "\nYears of %s warming"%second_delta,year4
cases=["2","4_2","4_2_2"]
def compute_change_fields() :
global threshold
fields=dict() # Returned dict of fields for plot : means of changes, stippling masks, hatching masks
# Strucure is :
# fields[season][case][choice]
# with 'choice' varying among : mean,median,rmean,rmedian,stippling,hatching
# where prefix 'r' means 'relative')
diffs =dict() # diffs[season][case][case][model] where 'case' is either 'plain' or 'relative'
if threshold is not None :
if type(threshold) is str :
threshold=eval(threshold)
threshold_string="%g"%threshold
metadata=""
for season in seasons:
print season+" : ",
fields[season]=dict()
#
if with_variability :
variabilities=cens()
#print control_models
for model,realization in control_models :
# Store model internal variability
variabilities[model]=regridn(
variability_AR5(model,realization,variable,table,data_versions,season=season,
**variability_sampling_args),
cdogrid=common_grid)
# Compute median variability across models
variability= ccdo_ens(variabilities ,operator="enspctl,50")
else:
variability=None
#
for model,realization in year4 :
#print model,
# Compute reference field
grid,version,_= data_versions[ref_experiment][variable][table][model][realization]
dref = dict(project="CMIP6", experiment=ref_experiment,
model=model, institute=institute_for_model(model),
period=ref_period, variable=variable, table=table,
version=version, grid=grid,mip=mip_for_experiment(ref_experiment),
realization=realization)
metadata+=TSU_metadata(ref_experiment,[(model,realization)],variable,table,data_versions)
ref = ccdo(ds(**dref),operator="timmean -selseason,%s"%season)
feed_dic(diffs,regridn(ref,cdogrid=common_grid),season,"2","ref",model)
#
# Move to projection experiment
dic = dref.copy()
_,version,_=data_versions[experiment][variable][table][model][realization]
dic.update(experiment=experiment,mip=mip_for_experiment(experiment), version=version)
metadata+=TSU_metadata(experiment,[(model,realization)],variable,table,data_versions)
#
# Compute field and changes at first 2°C
#
period = "%d-%d"%(year2[(model,realization)]-window_half_size,
year2[(model,realization)]+window_half_size)
dic.update(period=period)
rr2 = ccdo(ds(**dic),operator="timmean -selseason,%s"%season)
feed_dic(diffs,regridn(rr2,cdogrid=common_grid),season,"2","proj",model)
#
# plain change
diff=ccdo2(rr2,ref,operator="sub")
regridded_diff=regridn(diff,cdogrid=common_grid)
feed_dic(diffs,regridded_diff,season,"2","change",model)
#
# normalized change
if model in variabilities :
nchange=ccdo2(regridded_diff,variabilities[model],operator="div")
feed_dic(diffs,nchange,season,"2","nchange",model)
#
# relative change
if threshold is not None :
thresholded_ref=ccdo_fast(ref,operator="setrtomiss,-1.e+10,"+threshold_string)
else :
thresholded_ref=ref
rr2_relative = ccdo_fast(ccdo2(diff,thresholded_ref,operator="div"),operator="mulc,100.")
feed_dic(diffs,regridn(rr2_relative,cdogrid=common_grid),season,"2","rchange",model)
#
# Store proj value for 2K as ref value for case 4K-2K
feed_dic(diffs,regridn(rr2,cdogrid=common_grid),season,"4_2","ref",model)
#
# Compute field at 4°C, and changes vs first 2°C
#
period = "%d-%d"%(year4[(model,realization)]-10,year4[(model,realization)]+10)
dic.update(period=period)
rr4 = ccdo(ds(**dic),operator="timmean -selseason,%s"%season)
feed_dic(diffs,regridn(rr4,cdogrid=common_grid),season,"4_2","proj",model)
#
# plain change
diff=ccdo2(rr4,rr2,operator="sub")
regridded_diff=regridn(diff,cdogrid=common_grid)
feed_dic(diffs,regridded_diff,season,"4_2","change",model)
#
# normalized change
if model in variabilities :
nchange=ccdo2(regridded_diff,variabilities[model],operator="div")
feed_dic(diffs,nchange,season,"4_2","nchange",model)
#
# relative change from 2K to 4K
if threshold is not None :
thresholded_rr2=ccdo_fast(rr2,operator="setrtomiss,-1,"+threshold_string)
else :
thresholded_rr2=rr2
rr4_relative = ccdo_fast(ccdo2(diff,thresholded_rr2,operator="div"),operator="mulc,100.")
feed_dic(diffs,regridn(rr4_relative,cdogrid=common_grid),season,"4_2","rchange",model)
#
# Compute diff between 4K-2K and 2K
#
for opt in ["change","rchange","nchange"] :
if opt in diffs[season]["4_2"] and model in diffs[season]["4_2"][opt] :
tmp=ccdo2(diffs[season]["4_2"][opt][model], diffs[season]["2"][opt][model],
operator="sub")
feed_dic(diffs,tmp,season,"4_2_2",opt,model)
print
#
# Choose field type for computing stippling/hatching
if "mean" in field_type :
choice="mean" # For cases mean and rmean
else:
choice="median" # For cases median and rmedian
#
# Compute ensemble statistics, and stippling
for case in cases :
fields[season][case]=dict()
fields[season][case]["mean"] = ccdo_ens(cens(diffs[season][case]["change"]) ,
operator="ensmean")
fields[season][case]["median"] = ccdo_ens(cens(diffs[season][case]["change"]) ,
operator="enspctl,50")
fields[season][case]["rmean"] = ccdo_ens(cens(diffs[season][case]["rchange"]),
operator="ensmean")
fields[season][case]["rmedian"]= ccdo_ens(cens(diffs[season][case]["rchange"]),
operator="enspctl,50")
#
# Compute relative change of ensemble means
if case != "4_2_2" :
meanr=ccdo_ens(cens(diffs[season][case]["ref"]) , operator="ensmean")
mean2=ccdo_ens(cens(diffs[season][case]["proj"]) , operator="ensmean")
fields[season][case]["rmeans"]=ccdo2(ccdo2(mean2,meanr,operator="sub"),
meanr,operator="mulc,100 -div")
else :
fields[season][case]["rmeans"]=ccdo2(fields[season]["4_2"]["rmeans"],
fields[season]["2"]["rmeans"],operator="sub")
#
agreef = agreement_fraction_on_sign(cens(diffs[season][case]["change"]))
if scheme=="AR5" :
fields[season][case]["stippling"],fields[season][case]["hatching"]=\
stippling_hatching_masks_AR5(
fields[season][case][choice],variability,agreef)
ceval(fields[season][case]["stippling"])
elif scheme == "AR6" :
if "nchange" in diffs[season][case]:
nchanges=regridn(cens(diffs[season][case]["nchange"]),cdogrid=common_grid)
agree_low=agreement_fraction_on_lower(nchanges,confidence_factor,sign_threshold)
fields[season][case]["lowchange"],fields[season][case]["conflict"] = \
lowchange_conflict_masks_AR6(agreef,agree_low)
elif scheme == "AR6S" :
changes=regridn(cens(diffs[season][case]["change"]),cdogrid=common_grid)
agreef=agreement_fraction_on_sign(changes)
fields[season][case]["sign_mask"]=ccdo_fast(agreef, operator="lec,%g"%sign_threshold)
#
return fields
#
fields=compute_change_fields()
if not os.path.exists(outdir):
os.makedirs(outdir)
with open("%s/%s%s_md.txt"%(outdir,figure_name,version),"w") as f:
f.write(metadata)
fig_for_label=change_figure(
variable, variable_transformation,fields[seasons[0]]["2"][field_type],
relative=("rme" in field_type),labelbar="True",
custom_plot=plot_args,number=0,title="some_dummy_title")
#
figfile_for_label="./tmp_fig_for_label.png"
#os.system("rm -f %s"%figfile_for_label)
cdrop(fig_for_label)
cfile(fig_for_label,figfile_for_label,ln=True)
create_labelbar(figfile_for_label, "./insert.png",scheme=scheme)#,missing=False,y_offset=570)
#os.system("rm -f %s"%figfile_for_label)
plots=dict()
#
titles={ seasons[0] : { "2" : "a) %s, first 2~S~o~"%seasons[0],
"4_2" : "c) %s, second 2~S~o~"%seasons[0],
"4_2_2" : "e) %s, second 2~S~o~N~ - first 2~S~o~N~ ( c)-a) )"%seasons[0]},
seasons[1] : { "2" : "b) %s, first 2~S~o~"%seasons[1],
"4_2" : "d) %s, second 2~S~o~"%seasons[1],
"4_2_2" : "f) %s, second 2~S~o~N~ - first 2~S~o~N~ ( d)-b) )"%seasons[1]}}
#
for season in seasons :
plots[season]=dict()
for case in cases :
mask1="" ; mask2=""
pattern1="" ; pattern2=""
if scheme == "AR5" and "hatching" in fields[season][case]:
mask1=fields[season][case]["hatching"]
pattern1="hatching"
mask2=fields[season][case]["stippling"]
pattern2="stippling"
elif scheme == "AR6" and "conflict" in fields[season][case]:
mask1=fields[season][case]["conflict"]
pattern1="crosses"
mask2=fields[season][case]["lowchange"]
pattern2="backslashes"
elif scheme == "AR6S":
mask2=fields[season][case]["sign_mask"]
pattern2="slashes"
#
plots[season][case]=change_figure(
variable, variable_transformation,
fields[season][case][field_type],
mask1=mask1, pattern1=pattern1,
mask2=mask2, pattern2=pattern2,
relative=("rme" in field_type),
title=titles[season][case], number="%d"%len(year4), labelbar="False",
custom_plot=plot_args)
# Write fields
fn="%s=%s=%s=%s=%s=%s=%s=%s.nc"%\
(season,variable,experiment,field_type,variable_transformation,\
ref_period,case,figure_name+version)
cfile(fields[season][case][field_type],fn,ln=True)
#
fn="%s=%s=%s=%s=%s=%s=%s=%s.nc"%\
(season,variable,experiment,pattern1,variable_transformation,\
ref_period,case,figure_name+version)
if mask1 != "" : cfile(mask1,fn,ln=True)
#
fn="%s=%s=%s=%s=%s=%s=%s=%s.nc"%\
(season,variable,experiment,pattern2,variable_transformation,\
ref_period,case,figure_name+version)
if mask2 != "" : cfile(mask2,fn,ln=True)
#
plot_lines= [
[plots[seasons[0]]["2"] ,plots[seasons[1]]["2"] ],
[plots[seasons[0]]["4_2"] ,plots[seasons[1]]["4_2"] ],
[plots[seasons[0]]["4_2_2"],plots[seasons[1]]["4_2_2"]]]
#
fig=cpage(plot_lines,title=title, insert="./insert.png", **figure_details)
outfile="%s_change_2K_4K_2seasons_%s%s.png"%(variable,data_versions_tag,version)
cdrop(fig)
#os.system("rm -f %s"%(outdir+"/"+outfile))
cfile(fig,outdir+"/"+outfile)
os.system("cd %s ; ln -sf %s %s%s.png"%(outdir,outfile,figure_name,version))
#
small=outfile.replace(".png",".small.png")
os.system("cd %s ; convert -geometry 50%% %s %s"%(outdir,outfile,small))
os.system("cd %s ; ln -sf %s %s%s_small.png"%(outdir,small,figure_name,version))
#iplot(fig)
#iplot(plots["DJF"]["4_2"])
#iplot(plot(fields["DJF"]["2"]["mean"],scale=3600.*24.,min=-2,max=2,delta=0.4))