CAMMAC https://cammac.readthedocs.io

Build a figure showing long-term change for 1 variable and one SSP

S.Sénési for Météo-France - sept 2019 to march 2021

A few commands below are specific to the Notebook environment, and can be safely commented out

If using a notebook, use wide display

In [ ]:
from IPython.core.display import display, HTML, Image
display(HTML("<style>.container { width:100% !important; }</style>"))

Default settings (some may be overriden by Papermill - this would show in next cell in the execution output notebook)

In [ ]:
import os
figure_filename      = None # if set, will supersede the automatic one
# Version number will be a suffix for automatic figure filename. Use e.g. "_V1" for legibility
version              = ""
manual_title         = None  # If set, will replace automatic title 
metadata_file_prefix = "Fig_simple"

#See doc for data_versions in sibling directory data_versions
data_versions_tag   = "20201210_derived"
data_versions_dir   = os.getenv("CAMMAC")+"/data"
excluded_models     = []
included_models     = None   # If this is a list, only listed models will be used
#included_models     = ["CNRM-CM6-1","IPSL-CM6A-LR"]
variability_excluded_models= []
variability_models        = None


#
variable            = "pr"
#threshold           = 0.1/(24*3600) # 0.1 mm/day, in SI
threshold           = None
table               = "Amon"
derivation          = "plain"
variable_label      = "Precipitation"
#
experiment          = "ssp585"
proj_period         = "2081-2100"
#
ref_experiment      = "historical"
ref_period          = "1995-2014" 
#
season              = "DJF"        
field_type          = "means_rchange" 
print_statistics    = True   

# A set of predefined intervals adequate for % change
col100 = {"colors":"-100 -75 -50 -25 -10   0   10 25 50 75 100 "}
col40  = {"colors":" -40 -30 -20 -10 -5    0    5 10 20 30  40 "}
col20  = {"colors":" -20 -15 -10  -5 -2.5  0  2.5  5 10 15  20 "}
col10  = {"colors":" -10  -8  -6  -4 -2    0    2  4  6  8  10"}

do_plot             = True
scheme              = "AR6"
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

plot_args           = {"min":-50.,"max":50., "delta": 10.}
labelbar            = True

# If pre-computed fields for this SSP, season and projection_period are available, should we use it ?
# Set it to False for recomputing and printing fields and field changes statistics
# No problem if set to True and pre-computed fields does not exist
use_cached_proj_fields = False
drop_old_figures       = False
# Next can be {} to deactivate variability computation
variab_sampling_args= {"shift":100,"nyears":20,"number":20,"detrend":True,"house_keeping":False,"compute":True,}
#variab_sampling_args= {}
#
#
outdir              = "./figures"
cache_dir           = "./cache"
figure_details      = {"page_width":2450,"page_height":3444,"pt":48, "ybox":133,"y":52,"insert_width":2400}
common_grid         = "r360x180"


# If some basic fields are to be plotted for each model (one figure per field_type, one panel per model):
#   - should we restrict the plotted models to a given list (None means : plot all)
plot_only              = None
#   - which fields should be actually plotted
field_types_to_plot_for_all_models    = [ "reference" , "projection", "change", "rchange", "schange", "variability" ]
#   - 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. } ,
#   "change"     : { "min" :1000.,"max":-1000. , "delta":200.} , 
   "rchange"    : { "min" : 0., "max" : 1000. , "delta" : 100. ,"units" : "%" }#, "color":"AR6_Precip_12"} ,
#   "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 } ,
    }

# We can also plot for some 'selected' models, a single-model figure with one panel per field type 
models_to_plot = []
#models_to_plot = [ "CNRM-CM6-1"]
field_types_to_plot_for_selected_models = []  #e.g. ["rchange"]

# And a selection of field types in one panel per model
field_types_for_multi_types_plot = [ ]
models_for_multi_types_plot      = [ ]


do_test=True

# Expert CliMAF tuning  :
own_cache            = False
one_cache_per_process= False
cache_root           = None
raz                  = False
dig                  = False
In [ ]:
if do_test:
    #included_models     = ["CNRM-CM6-1", "IPSL-CM6-1","BCC-ESM1","CESM2-WACCM","CanESM5"]
    #included_models     = ["CNRM-CM6-1"]
    #included_models     = ["EC-Earth3" ]
    #variability_models        = included_models
    #variability_excluded_models= []
    #proj_period         = "2081"
    #ref_period          = "2014" 
    #plot_args           = {"min":-50.,"max":50., ""delta"": 10., "options": "gsnRightStringFontHeightF=2."}
    experiment           = "ssp370"
    scheme               = "AR6"
    season               = "ANN"
    #use_cached_proj_fields = False
    #variab_sampling_args= {}
    variab_sampling_args= {"house_keeping":False,"compute":True,"detrend":True,"shift":0,"nyears":20,"number":10}
    #field_types_to_plot_for_all_models = []
    own_cache            = False
    one_cache_per_process= False
    raz                  = False
    dig                  = False
    print_statistics=True
    #models_for_multi_types_plot      = [ "CNRM-CM6-1"]
    #field_types_for_multi_types_plot = [ "reference", "projection", "change_orig" , "rchange_orig"]
    ranges = {
      "reference"   : { "scale" : 24*3600 , "units" : "mm/d", "min" : 0 , "max" : 10 , "delta" : 0.5 },
      "projection"  : { "scale" : 24*3600 , "units" : "mm/d", "min" : 0 , "max" : 10 , "delta" : 0.5 },
      "change"      : { "scale" : 24*3600 , "units" : "mm/d", "min" : -10 , "max" : 10 , "delta" : 1 },
      "change_orig" : { "scale" : 24*3600 , "units" : "mm/d", "min" : -10 , "max" : 10 , "delta" : 1 },
      "rchange"     : { "min" : -100. , "max" : 100., "delta" : 10. },
      "rchange_orig": { "min" : -100. , "max" : 100., "delta" : 10. },
      "variability" : { "scale" : 24*3600 , "units" : "mm/d", "min" : 0 , "max" : 1 , "delta" : 0.05 },
    }
In [ ]:
if (len(field_types_to_plot_for_all_models) > 0 or len(field_types_to_plot_for_selected_models) > 0) and use_cached_proj_fields :
    raise ValueError("Must set use_cached_proj_fields to False for plotting individual models ")

Loading libraries

In [ ]:
import sys, os, os.path, time, resource

# Climaf setup (version >= 1.2.13 - see https://climaf.readthedocs.io)
if own_cache or one_cache_per_process:
    if cache_root is None :
        cache_root=os.getcwd()
    cache_dir=cache_root+"/tmp_climaf_cache"
    if one_cache_per_process :
        cache_dir += "_%d"%os.getpid()
    os.environ["CLIMAF_CACHE"]=cache_dir
from climaf.api import *
if own_cache or one_cache_per_process:
    climaf.cache.currentCache=cache_dir

# Climaf settings
if raz and (own_cache or one_cache_per_process): 
    craz()
climaf.cache.stamping             = False
climaf.driver.dig_hard_into_cache = dig

from CAMMAClib.changes     import change_figure_with_caching
from CAMMAClib.ancillary   import prettier_label, create_labelbar
from CAMMAClib.mips_et_al  import TSU_metadata, read_versions_dictionnary
# 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'))
In [ ]:
if not os.path.exists(outdir):    os.makedirs(outdir)
In [ ]:
do_check=False
if do_check :
    base_dict=dict(project="CMIP6", experiment="historical",
                        model="CNRM-CM6-1", institute="CNRM-CERFACS",
                        period="1995-2014", variable="pr", table="Amon", 
                        version="latest", grid="*",realization="r1i1p1f2")
    r=ccdo_fast(ds(**base_dict),operator="timmean -selseason,%s -seasmean"%season)
    #Image(cfile(plot(r)))
    #
    thresh=0.0001/(24*3600)
    s=ccdo_fast(r,operator="setrtomiss,-%g,%g"%(thresh,thresh))
    Image(cfile(plot(s)))

Next function combines global variables and its own arguments for calling change_figure_with_caching

In [ ]:
def afigure(variable, title, panel,deep=None,same=True,scheme=scheme):

    global metadata
    print(included_models)
    fig_file,fig,dic,variab_models,models = change_figure_with_caching(
        ref_period = ref_period, proj_period = proj_period, 
        variable = variable, table = table, ref_experiment = ref_experiment,
        experiment = experiment, season = season,
        derivation_label = derivation, 
        field_type = field_type,
        title = title, 
        custom_plot = plot_args, labelbar = labelbar, 
        data_versions_tag = data_versions_tag, data_versions_dir = data_versions_dir,
        excluded_models = excluded_models,models=included_models,
        outdir = outdir, drop=drop_old_figures,
        #
        common_grid = common_grid, 
        variab_sampling_args = variab_sampling_args,
        variability_models=variability_models, variability_excluded_models=variability_excluded_models,
        cache_dir = cache_dir, read = use_cached_proj_fields, write = True, 
        print_statistics =  print_statistics, deep = deep,  
        same_models_for_variability_and_changes=same,threshold=threshold,
        scheme=scheme,
        change_sign_agree_threshold = sign_threshold
        )
    metadata += TSU_metadata([ref_experiment,experiment],models,       variable,table,data_versions,panel)
    metadata += TSU_metadata(["piControl"],              variab_models,variable,table,data_versions,panel)
    return fig,dic
In [ ]:
def mem() :
    return "%d Mo"%int(resource.getrusage(resource.RUSAGE_SELF).ru_maxrss/1024)

def timer(step="init"):
    global previous
    now=time.time()
    if step=='init' :
        print "At %s : %s"%(step,mem())
    else :
        print "At %s : %d s, %s"%(step,now-previous,mem())
    previous=now 

Compute figure

In [ ]:
data_versions=read_versions_dictionnary(data_versions_tag, data_versions_dir)
metadata=""

if manual_title is None :
    title="Change for %s %s and %s (%s vs %s)"%(variable,prettier_label.get(experiment,experiment),season,proj_period,ref_period)
else:
    title=manual_title
clog('info')
timer('init')
clog('error')
fig,dic=afigure(variable, title,"",same=False,scheme=scheme)
#print dic
timer("fig")
In [ ]:
#ncview(dic['ssp585']['DJF']['projection']['plain']['EC-Earth3'])

Write figure and metadata

In [ ]:
with open("%s/%s_md.txt"%(outdir,metadata_file_prefix),"w") as f:
    f.write(metadata)

if figure_filename is None :
    outfile="%s/change_map_%s_%s_%s_%s_%s_%s%s.png"%(outdir,variable,derivation,field_type,experiment,proj_period,data_versions_tag,version)
else :
    outfile="%s_%s.png"%(figure_filename,data_versions_tag)

#craz()
if do_plot :
    cdrop(fig)
    cfile(fig,outfile)
    timer("apres")
    print ("Figure available : %s"%outfile)
    os.system("convert -geometry 50%% %s %s"%(outfile,outfile.replace(".png",".small.png")))
#print cfile(cpage([[fig]],title="try"))
In [ ]:
if do_plot :
    # Call utility function for extracting labelbar and adding fill pattern for shadings
    create_labelbar(outfile, "./insert.png")
    figure_title="Just for testing new legend"
    # Create multi-panel figure
    page=cpage([[fig]],title=figure_title, insert="./insert.png", **figure_details)
    cdrop(page)
    #
    outfile2=outfile.replace(".png",".page.png")
    cfile(page,outdir+"/"+outfile2)

If using a notebook , display result

In [ ]:
#Image(outfile,width=800)

Ploting a field type for all models - iterate on some field types

In [ ]:
if figure_details_all_models is None :
    figure_details_all_models=figure_details
In [ ]:
for field_type in field_types_to_plot_for_all_models :
    #
    if field_type in dic[experiment][season]:
        plotargs=custom_plot_all_fields.copy()
        plotargs.update(ranges.get(field_type,{}))
        #if "_orig" in field_type :
        plotargs.update(options="cnFillMode=CellFill")
        #
        if plot_only is not None :
            ens=cens()
            for model in plot_only :
                if model in dic[experiment][season][field_type][derivation] :
                    ens[model]=dic[experiment][season][field_type][derivation][model]
        else:
            ens=cens(dic[experiment][season][field_type][derivation])
        allplots=plot(ens,**plotargs)
        page=cpage(allplots,title="%s_%s_%s_%s_%s"%(variable,derivation,field_type,experiment,season),**figure_details_all_models)
        pagename="%s/all_models_%s_%s_%s_%s_%s_%s.png"%(outdir,variable,derivation,field_type,experiment,season,data_versions_tag)
        cfile(page,pagename)
        

Ploting all field types for a model - iterate on some models

In [ ]:
#clog("debug")
for model in models_for_multi_types_plot :
    #
    ens=cens()
    for field_type in field_types_for_multi_types_plot :
        if field_type in dic[experiment][season] :
            plotargs=custom_plot_all_fields.copy()
            plotargs.update(ranges.get(field_type,{}))
            
            plotargs.update(options="cnFillMode=CellFill")
            #
            if model in dic[experiment][season][field_type][derivation] :
                field=dic[experiment][season][field_type][derivation][model]
                ens[field_type]=plot(field,title=field_type,**plotargs)
                print "%s plotted for %s"%(field_type,model)
    #
    page=cpage(ens,title="%s_%s_%s_%s_%s"%(model,variable,derivation,experiment,season),**figure_details_all_models)
    pagename="%s/all_fields_%s_%s_%s_%s_%s_%s.png"%(outdir,model,variable,derivation,experiment,season,data_versions_tag)
    cfile(page,pagename)
#Image(cfile(ens[model]))
        
In [ ]: