S.Sénési for Météo-France - sept 2019 to march 2021
import os
import os #
do_test = False
#
title = "Multi-model zonal mean long-term changes in P, E and P-E"
variables = ["pr", "evspsbl", "P-E"]
unit = "mm/day"
factor = 24*3600
table = "Amon"
columns = "Precipitation|Evaporation|Precip. minus evaporation"
figure_name = "Fig8-14" # Used for a symbolic link to an explicit filename,and for metadata filename
version = "" # Suffix added to figure name
outdir = "./figures"
# Data used
data_versions_tag = "20200719"
data_versions_dir = os.getenv("CAMMAC")+"/data"
excluded_models = {}
included_models = None
variability_excluded_models= { "P-E" : "ACCESS-ESM1-5"} # No common data period between pr and evspsbl as of 20200913
variability_models = None
default_fixed_fields_dir = os.getenv("CAMMAC")+"/data/fixed_fields"
#
# If variables list or next series of parameters is changed, the computation
# will still be OK in intermediate files but the companion plot script may
# have to be modified
#
ref_experiment = "historical"
ref_period = "1995-2014"
# The plot script will set correct colors only if experiment order is ssp126/ssp245/ssp585
experiments = ["ssp126", "ssp245", "ssp585"]
projection_period = "2081-2100"
#
# Ancillary parameters
#
print_statistics = True
common_grid = "r360x180"
variability_sampling_args = {"house_keeping":False,"compute":True,"detrend":True,"shift":100,"nyears":20,"number":20}
check_fixed_fields = True
#
#
# Toggle for getting profiling information
do_profiling = False
if do_test :
variables = ["pr"]
version = "_test"
#variables = ["pr", "pr", "pr"]
#columns = "Precipitation|Precipitation|Precipitation"
#included_models = ["IPSL-CM6A-LR","CESM2-WACCM","CanESM5"]
included_models = ["IPSL-CM6A-LR"]
variability_models = included_models
#ref_period = "2014"
# The plot script will set correct colors only if experiment order is ssp126/ssp245/ssp585
#experiments = ["ssp126", "ssp126", "ssp126"]
#projection_period = "2099-2100"
#variability_sampling_args = {"house_keeping":False,"compute":True,"detrend":True,"shift":100,"nyears":3,"number":3}
from IPython.core.display import display, HTML, Image
display(HTML("<style>.container { width:100% !important; }</style>"))
if do_profiling:
import cProfile
import io
import pstats
import sys
import os.path
from climaf.api import *
climaf.cache.stamping=False
from CAMMAClib.mips_et_al import read_versions_dictionnary, table_for_var_and_experiment, \
institute_for_model, models_for_experiments, TSU_metadata
from CAMMAClib.ancillary import prettier_label, feed_dic
from CAMMAClib.variability import variability_AR5
# 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'))
data_versions=read_versions_dictionnary(data_versions_tag,data_versions_dir)
if check_fixed_fields :
issue=False
for variable in variables :
for experiment in experiments :
for model,realization in models_for_experiments(data_versions,variable,table,[experiment],
excluded_models.get(variable,[]),included_models) :
grid,version,data_period=data_versions[experiment][variable][table][model][realization]
base_dict=dict(project="CMIP6", experiment=experiment,
model=model, institute=institute_for_model(model),
version=version, mip="ScenarioMIP",realization=realization,
table="fx",period="fx",variable="sftlf",grid=grid)
try :
cfile(ds(**base_dict))
except :
print "Issue for fixed field sftlf and model %s"%model
issue=True
if issue :
raise ValueError("Some sftlf fields are missing")
if do_profiling:
pr = cProfile.Profile()
pr.enable()
# Put intermediate results in dictionnaries
all_anomalies=dict() # anomalies of piControl periods, for all models
references=dict()
projections=dict()
changes=dict()
changes_land=dict()
count=dict()
files=""
metadata=""
panel={ "pr" : { "ssp126": "a","ssp245": "d","ssp585": "g", },
"evspsbl": { "ssp126": "b","ssp245": "e","ssp585": "h", },
"P-E" : { "ssp126": "c","ssp245": "f","ssp585": "i", }}
# Some operations have to be done on all dicts above
dicts=[all_anomalies,references,projections,changes,changes_land]
#
if print_statistics :
print "Values below are field medians"
print "exp. variable model reference projection change change_land"
print 120*"-"
#
for experiment in experiments :
for d in dicts : d[experiment]=dict()
count[experiment]=dict()
#
for variable in variables :
table=table_for_var_and_experiment(variable,experiment)
#
for d in dicts : d[experiment][variable]=cens()
count[experiment][variable]=0
#
# Populate an ensemble merging all piControl variability samples among models
#
control_pairs=models_for_experiments(data_versions,variable,table,["piControl"],
variability_excluded_models.get(variable,[]),variability_models)
if len(control_pairs) == 0 :
raise ValueError("No model provides %s for piControl according to the versions dictionnary, among"%variable+`variability_models`)
for model,realization in control_pairs :
print "%s %4s %-20s"%(experiment,variable, model),
metadata+=TSU_metadata("piControl",[(model,realization)],variable,table,data_versions,panel[variable][experiment])
sample=variability_AR5(model,realization,variable,table,data_versions, season="anm", project="CMIP6",
variability=False, **variability_sampling_args)
sample_mean=ccdo_ens(sample,operator="ensmean")
departures=ccdo2(sample,sample_mean,operator="sub")
departuresr=regridn(departures,cdogrid=common_grid)
for cperiod in sample :
all_anomalies[experiment][variable][model+"_"+cperiod]=departuresr[cperiod]
#
pairs=models_for_experiments(data_versions,variable,table,["historical",experiment],
excluded_models.get(variable,[]),included_models)
if len(pairs) == 0 :
raise ValueError("No model provides %s according to the versions dictionnary, among"%variable+`included_models`)
for model,realization in pairs :
print "%s %4s %-20s"%(experiment,variable, model),
#
# Build a dictionnary of facets for reference dataset
#
grid,version,_ = data_versions[ref_experiment][variable][table][model][realization]
base_dict=dict(project="CMIP6", experiment=ref_experiment,
model=model, institute=institute_for_model(model),
period=ref_period, variable=variable, table=table,
version=version, mip="CMIP",realization=realization,grid=grid)
metadata+=TSU_metadata(ref_experiment,[(model,realization)],variable,table,data_versions,panel[variable][experiment])
#
# Compute reference time mean over requested season
#
reference_dict=base_dict.copy()
reference=clim_average_fast(ds(**reference_dict),"anm")
references[experiment][variable][model]=reference
#
# Compute projection time mean over requested season
#
projection_dict=reference_dict.copy()
_,version,_ = data_versions[experiment][variable][table][model][realization]
projection_dict.update(mip="ScenarioMIP",experiment=experiment,
period=projection_period, realization=realization,version=version)
metadata+=TSU_metadata(experiment,[(model,realization)],variable,table,data_versions,panel[variable][experiment])
projection=clim_average_fast(ds(**projection_dict),'anm')
projections[experiment][variable][model]=projection
#
# Compute change and regrid to common grid
#
change=ccdo2(projection,reference,operator="sub")
changes[experiment][variable][model]=regridn(change,cdogrid=common_grid)
#
# Same for land only
#
# Compute a land mask with 1 if land, else missing
sftlf=base_dict.copy()
sftlf.update(table="fx",period="fx",variable="sftlf",experiment="piControl")
land=ccdo_fast(ccdo_fast(ds(**sftlf),operator="setrtomiss,-1,99.9"),operator="divc,100")
# Apply land mask
change_land=ccdo2(change,land,operator="mul")
changes_land[experiment][variable][model]=regridn(change_land,cdogrid=common_grid)
#
#
if print_statistics :
print " %7.2g %7.2g %7.2g %7.2g "%(
cvalue(ccdo_fast(reference,operator="fldpctl,50")),\
cvalue(ccdo_fast(projection,operator="fldpctl,50")),\
cvalue(ccdo_fast(change,operator="fldpctl,50")),\
cvalue(ccdo_fast(change_land,operator="fldpctl,50"))),
count[experiment][variable] = count[experiment][variable] + 1
print
#if not do_test : csync()
#
# Compute ensemble statistics over models
#
#print "keys=",changes[experiment][variable].keys()
#print "values=",changes[experiment][variable].values()
ensavg = ccdo_ens(changes[experiment][variable],operator='ensmean')
ensstd1 = ccdo_ens(changes[experiment][variable],operator='ensstd1')
ensavg5 = ccdo2(ensavg,ccdo_fast(ensstd1,operator="mulc,1.645"),operator="sub")
ensavg95= ccdo2(ensavg,ccdo_fast(ensstd1,operator="mulc,1.645"),operator="add")
ensavg_land=ccdo_ens(changes_land[experiment][variable],operator='ensmean')
#
# Same on ensemble of anomalies for control run (ensemble over time slices and models)
#
variavg = ccdo_ens(all_anomalies[experiment][variable],operator='ensmean')
varistd1 = ccdo_ens(all_anomalies[experiment][variable],operator='ensstd1')
variab5 = ccdo2(variavg,ccdo_fast(varistd1,operator="mulc,1.645"),operator="sub")
variab95 = ccdo2(variavg,ccdo_fast(varistd1,operator="mulc,1.645"),operator="add")
#
# Compute zonal means, group it as a CliMAF ensemble of 1d fields
#
lines=cens()
lines["pctl5"] = ccdo_fast(ensavg5,operator="zonmean")
lines["variab5"] = ccdo_fast(variab5,operator="zonmean")
lines["mean"] = ccdo_fast(ensavg,operator="zonmean")
lines["variab95"] = ccdo_fast(variab95,operator="zonmean")
lines["pctl95"] = ccdo_fast(ensavg95,operator="zonmean")
lines["land_mean"] = ccdo_fast(ensavg_land,operator="zonmean")
#
# Create data file for this experiment and variable
#
fn="out_"+experiment+"_"+variable+".nc"
efile(lines,fn,force=True)
files=files+" "+fn
print
# Make sure to write CliMAF cache index, accounting for all these new fields
#if not do_test : csync()
print "Done"
if do_profiling:
pr.disable()
s = io.BytesIO()
ps = pstats.Stats(pr, stream=s).sort_stats('cumulative')
ps.print_stats(40)
print s.getvalue()
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)
#count={'ssp126': {'pr': 19, 'P-E': 18}, 'ssp585': {'pr': 19, 'P-E': 19}, 'ssp245': {'pr': 21, 'P-E': 20}}
files="out_ssp126_pr.nc out_ssp126_evspsbl.nc out_ssp126_P-E.nc out_ssp245_pr.nc out_ssp245_evspsbl.nc out_ssp245_P-E.nc out_ssp585_pr.nc out_ssp585_evspsbl.nc out_ssp585_P-E.nc" columns="Precipitation|Evaporation|Precipitation minus evaporation" lines="SSP1-1.9 (mm/day)|SSP2-4.5 (mm/day)|SSP5-8.5 (mm/day)" labels="3 models| 3 models| 3 models| 3 models| 3 models| 3 models| 3 models| 3 models| 3 models" title="title" nx=3 ny=3 command="ncl -Q "+CAMMAC+"/notebooks/change_zonal_mean.ncl "+\ "'title=\"%s\"' 'files=\"%s\"' 'nx=%d' 'ny=%d' "%(title,files,nx,ny) +\ "'columns=\"%s\"' 'lines=\"%s\"' 'graphs=\"%s\"' 'format=\"1200x1800\"' "%(columns,lines,labels) print "executing : ",command ! {command}
figure_filename="zonal_mean_%s_%s_%s_%s%s.png"%(variables[0],variables[1],variables[2],data_versions_tag,version)
#
labels=""
lines=""
for experiment in experiments :
lines=lines+"%s (%s)|"%(prettier_label.get(experiment,experiment),unit)
for variable in variables :
labels+=" %s models|"%count[experiment][variable]
labels=labels[:-1]
lines=lines[:-1]
#
nx=len(variables)
ny=len(experiments)
#
#print "files=",files
#print "columns=",columns
#print "lines=",lines
#print "labels=",labels
#
command="ncl -Q "+CAMMAC+"/notebooks/change_zonal_mean.ncl "+\
"'title=\"%s\"' 'files=\"%s\"' 'nx=%d' 'ny=%d' "%(title,files,nx,ny) +\
"'columns=\"%s\"' 'lines=\"%s\"' 'graphs=\"%s\"' 'format=\"1200x1800\"' "%(columns,lines,labels)+\
"'fact=%g' "%factor
print "executing : ",command
! rm -f {figure_filename}
! {command}
! mv zonal_means.png {outdir}/{figure_filename}
! ln -sf {figure_filename} {outdir}/{figure_name}.png
#
Image(figure_filename)
#pb=changes_land["ssp585"]["pr"]
#for model in pb :
# zm=ccdo_fast(pb[model],operator="zonmean")