S.Sénési for Météo-France - sept 2019 to march 2021
import os
figure_name = "FigP-E@+3K"
version = "" # a suffix for figure filename. Use e.g. "_V1" for legibility
# Figure title
#manual_title = "Effect of a 3 degrees warming on P-E (vs 1850-1900)"
manual_title = None
outdir = "./figures"
#See doc for data_versions in sibling directory data_versions
data_versions_tag = "20200918"
excluded_models = [ ]
included_models = None
data_versions_dir = os.getenv("CAMMAC")+"/data"
variable = "P-E"
table = "Amon" # Script was yet tested only for a monthly table
variable_transformation = "plain" # Could be 'iav', 'gini', 'welsh', 'dry'...
season = "ANN" #, "DJF","JJA" # any CDO season. Graph is tuned for showing 1 season
experiment = "ssp585"
#
ref_experiment = "historical"
ref_period = "1850-1900"
proj_period = "2015-2099" # period investigated for the warming
warming = 3.0 # Temperature change (degrees)
window_half_size = 10 # For time filtering of atmospheric temperature before analyzing warming (unit=year)
#
clim_experiment = "piControl"
clim_period = "1-100" # period for the climatology. For piControl, provide years relative to begin (starting with 1)
#
field_type = "rmeans" # Type of change field plotted : 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. ")
clim_contours = [ -2, 0 , 2 ] # mm/day
#
# Other details
figure_details = {"page_width":2450,"page_height":3444, "insert_width":2000,"pt":60, "ybox":133,"y":40}
common_grid = "r360x180"
#
# If some basic fields are to be plotted for ~ each model :
# - should we restrict the plotted models to given list (None means : plot all)
plot_only = None
# - which fields should be actually plotted
field_types_to_plot_for_all_models = [ "clim" , "rchange"]
#field_types_to_plot_for_all_models = [ "clim" , "rmeans"]
# - with which common plot_parameters
custom_plot_all_fields = { "proj" : "Robinson", "mpCenterLonF" : 0., "options" : "lbBoxEndCapStyle=TriangleBothEnds" }#, "focus":"land"}
# - should we use specific settings for page layout
figure_details_all_models = None
# - which range should be used
ranges = {} # The baseline value !
ranges={
# "reference" : { "min" : 0., "max" : 3000. , "delta" : 200. } ,
# "projection" : { "min" : 0., "max" : 3000. , "delta" : 200. } ,
"clim" : { "contours":"-3 -2 -1 0 1 2 3", "scale" : 24*3600, "units" : "mm/d", "min":-4. , "max":4., "delta":1} ,
"rchange" : { "color" : "AR6_Precip_12", "colors":"-80. -40. -20. -10. -5. 0 5. 10. 20. 40. 80. " } ,
# "schange" : { "colors": "-5 -2 -1 -0.5 -0.25 0. 0.25 0.5 1 2 5" , "units":"-", "color":"AR6_Precip_12" } ,
# "variability": { "min" : 0., "max" : 1. , "delta" : 0.1 } ,
}
contour_mask = CAMMAC+"/data/fixed_fields/sftlf_fx_CNRM-CM6-1_historical_r1i1p1f2_gr.nc"
do_test = True
if do_test :
version = "_test"
ref_period = "1850"
clim_period = "1-2"
included_models = ["CNRM-ESM2-1"]
import sys, os
from climaf.api import *
from climaf.cache import stamping
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
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,clim_experiment],excluded_models,included_models)
#print models
if False and variable=="P-E" :
# Check that all models+variant have a 'latest' version of evspsbl
issues=[]
for model,variant in models :
for exp in [ref_experiment,experiment,clim_experiment] :
check=ds(project="CMIP6", experiment=exp,
model=model, institute=institute_for_model(model),
period="*", variable="evspsbl", table="Amon",
version="latest", grid=grid,mip=mip_for_experiment(exp),
realization=variant)
if check.baseFiles() is None or len(check.baseFiles())==0 :
issues.append((model,variant,exp))
if len(issues) > 0 :
print "Issues with evspsbl for : ",issues
raise ValueError("")
# compute ensemble of warming time 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=[]
years=dict()
for model,variant in models :
max_change[model]=cvalue(ccdo_fast(GSAT[model],operator="timmax"))
if max_change[model]>= warming :
models_warming_enough.append((model,max_change[model]))
metadata+=TSU_metadata([ref_experiment,experiment],[(model,variant)],"tas","Amon",data_versions)
year=int(proj_period.split("-")[0]) + window_half_size
found=False
for v in cMA(GSAT[model]).flatten().data :
if v >= warming and not found :
years[(model,variant)]=year
found=True
break
year+=1
else :
models_not_warming_enough.append((model,max_change[model]))
print "\nThese models don't reach %s K warming"%warming,models_not_warming_enough
print "\nThese (%d) models DO reach %s K warming"%(len(models_warming_enough),warming), models_warming_enough
print "\nYears when %s warming reached : "%warming,years
def compute_change_fields() :
#
global metadata
#
if threshold is not None :
threshold_string="%f"%threshold
#
# diffs[model][case] stores the 'ref' field, 'proj' field, 'change' and 'rchange' (relative change) fields for each model
diffs =dict()
#
for model,realization in years :
#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),"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 prescribed warming level
#
period = "%d-%d"%(years[(model,realization)]-window_half_size,years[(model,realization)]+window_half_size)
dic.update(period=period)
rr = ccdo(ds(**dic),operator="timmean -selseason,%s"%season)
feed_dic(diffs,regridn(rr,cdogrid=common_grid),"proj",model)
#
# plain change
diff=ccdo2(rr,ref,operator="sub")
feed_dic(diffs,regridn(diff,cdogrid=common_grid),"change",model)
#
# relative change
if threshold is not None :
thresholded_ref=ccdo_fast(ref,operator="setrtomiss,-1.e+10,"+threshold_string)
else :
thresholded_ref=ref
rr_relative = ccdo_fast(ccdo2(diff,thresholded_ref,operator="div"),operator="mulc,100.")
feed_dic(diffs,regridn(rr_relative,cdogrid=common_grid),"rchange",model)
print
#
# Compute climatology on prescribed experiment
#
if len(clim_contours) > 0 :
clims = models_for_experiments(data_versions,variable,table,[clim_experiment],included_models=[model])
if len(clims) != 1 :
raise ValueError("Issue with data for the climatology for model %s : %s"%(model,clims))
else :
_,clim_realization=clims[0]
_,version,clim_data_period = data_versions[clim_experiment][variable][table][model][clim_realization]
dclim = dic.copy()
dclim.update(experiment = clim_experiment, mip = mip_for_experiment(clim_experiment),
version = version, realization = clim_realization)
metadata += TSU_metadata(clim_experiment,[(model,clim_realization)],variable,table,data_versions)
if clim_experiment != "piControl" :
dclim.update(period=clim_period)
else :
shifts=clim_period.split("-")
clim_data_begin=int(clim_data_period.split("-")[0])
clim_begin = clim_data_begin+int(shifts[0])-1
clim_end = clim_data_begin+int(shifts[1])-1
dclim.update(period="%d-%d"%(clim_begin,clim_end))
clim = ccdo(ds(**dclim),operator="timmean -selseason,%s"%season)
feed_dic(diffs,regridn(clim,cdogrid=common_grid),"clim",model)
#
#
# Compute ensemble statistics, and stippling
# A dict of ensemble 'stat' fields, with 'stat' varying among :
# mean,median,rmean,rmedian,stippling,hatching (where prefix 'r' means 'relative')
fields=dict()
fields["mean"] = ccdo_ens(cens(diffs["change"]) , operator="ensmean")
fields["median"] = ccdo_ens(cens(diffs["change"]) , operator="enspctl,50")
fields["rmean"] = ccdo_ens(cens(diffs["rchange"]), operator="ensmean")
fields["rmedian"]= ccdo_ens(cens(diffs["rchange"]), operator="enspctl,50")
if len(clim_contours) > 0 :
fields["clim"] = ccdo_ens(cens(diffs["clim"]) , operator="ensmean")
#
# Compute relative change of ensemble means
meanr=ccdo_ens(cens(diffs["ref"]) , operator="ensmean")
mean2=ccdo_ens(cens(diffs["proj"]) , operator="ensmean")
fields["rmeans"]=ccdo2(ccdo2(mean2,meanr,operator="sub"),meanr,operator="mulc,100 -div")
#
return diffs,fields
#
diffs,fields=compute_change_fields()
#ncview(fields["clim"])
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)
if len(clim_contours) > 0 :
contours = fields["clim"]
mm_per_day = 1./(24*3600) # One millimeter per day in S.I.
plot_args["contours"] = ""
for c in clim_contours :
plot_args["contours"] += "%g "%(c*mm_per_day)
plot_args["aux_options"] = "cnLineThicknessF=6.|cnLineColor=black|gsnContourZeroLineThicknessF=12."#+plot_args.get("aux_options","")
contour_suffix = "_cont"
# Apply mask on contour fields
if contour_mask is not None :
mask = fds(contour_mask)
invert_mask = ccdo(mask, operator="expr,'invert=100-sftlf'")
invert_mask = ccdo(invert_mask, operator="gec,100")
actual_mask = regridn(invert_mask,cdogrid=common_grid)
contours = ccdo2_flip(contours,actual_mask,operator="ifthen")
else :
contours = ""
contour_suffix = ""
#print plot_args ["aux_options"]
the_plot=change_figure(variable, variable_transformation, fields[field_type],
hatching=contours, shade=False,
relative=("rme" in field_type),title="", number="%d"%len(years), labelbar="True",
custom_plot=plot_args)
#
fn="%s=%s=%s=%s=%s=%s=%s=%s.nc"%\
(season,variable,experiment,field_type,variable_transformation,\
ref_period,proj_period,figure_name+version),
cfile( fields[field_type],fn,ln=True)
fn="%s=%s=%s=%s=%s=%s=%s=%s.nc"%\
(season,variable,clim_experiment,field_type,variable_transformation,\
clim_period,"NA",figure_name+version),
cfile( contours,fn,ln=True)
#cdrop(the_plot)
#iplot(the_plot)
if manual_title is None :
title = "Effect on %s %s of a %d degrees warming (vs %s)"%(season,variable,int(warming),ref_period)
else :
title=manual_title
fig=cpage([[the_plot]],title=title, **figure_details)
outfile="%s_change_%dK_%s_%s%s_%s%s.png"%(variable,int(warming),experiment,season,contour_suffix,data_versions_tag,version)
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+"_"+season,version))
outfile="%s_%s_%dK_%s_%s_%s%s.nc"%(variable,field_type,int(warming),experiment,season,data_versions_tag,version)
cfile(fields[field_type],outfile,ln=True)
#
outfile="%s_clim_%dK_%s_%s_%s%s.nc"%(variable,int(warming),clim_experiment,season,data_versions_tag,version)
cfile(fields["clim"],outfile,ln=True)
if figure_details_all_models is None :
figure_details_all_models=figure_details
for field_type in field_types_to_plot_for_all_models :
#
if field_type in diffs:
plotargs=custom_plot_all_fields.copy()
plotargs.update(ranges.get(field_type,{}))
plotargs.update(options="cnFillMode=CellFill")
#
if plot_only is not None :
ens=cens()
for model in plot_only :
if model in diffs[field_type] :
ens[model]=diffs[field_type][model]
else:
ens=cens(diffs[field_type])
allplots=plot(ens,**plotargs)
page=cpage(allplots,title="%s_%s_%s_%s"%(variable,field_type,experiment,season),**figure_details_all_models)
pagename="%s/all_models_%s_%s_%s_%s_%s.png"%(outdir,variable,field_type,experiment,season,data_versions_tag)
cfile(page,pagename)