S.Sénési for Météo-France - sept 2019 to march 2021
fn="%s/statsCMIP6%s%s.json"%(input_dir,data_versions_tag,input_version)
filen="%s/%s_%s%s.json"%(outdir,"Hydro_vars_change",data_versions_tag,output_version)
from __future__ import print_function
import os
input_version = ""
input_dir = "./changes"
outdir = "./figures"
output_version = ""
#
data_versions_tag = "20210201_derived"
excluded_models = []
variables = [("pr","mean"),("mrro","mean"), ("pr","std"), ("mrro","std")]
max_warming = 5.
do_test = True
if do_test :
input_version = ""
output_version = "_tropic"
variables = [("pr","mean"),("mrro","mean"), ("pr","std"), ("mrro","std")]
#variables = [("mrro","mean")]
input_dir = "/home/ssenesi/CAMMAC/notebooks/prod_20210201/FigBoxTS.X_f3_h/changes"
outdir = input_dir
do_test_extra = False
if do_test_extra :
input_version = "extra"
output_version = "_extra"
variables = [("pr","mean"),("mrro","mean"), ("pr","std"), ("mrro","std")]
#variables = [("mrro","mean")]
input_dir = "/home/ssenesi/CAMMAC/notebooks/prod_20210201/FigBoxTS.X_f3_h_extra/changes"
outdir = input_dir
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 os
import json
import copy
def ensemble_stat_series(stats,variable,stat,season,option) :
""" Compute mean or quantile on an ensemble """
#
series=stats[variable][stat][season]["ens"]
periods=series.keys()
periods.sort()
ret=[]
for p in periods :
if float(p)> max_warming :
#print("skipping %g",float(p))
ret.append(missing_value)
continue
ens=series[p]
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 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,input_version)
print(fn)
with open (fn,"r") as f :
stats=json.load(f)
#print(stats.keys())
ref_period=stats["ref_period"]
missing_value=1.e+20
scenarios=["ssp585"]
#scenarios=None
models=filter_on_scenarios_models_warming(stats,scenarios,excluded_models,max_warming)
print(len(models),models)
# stat[variable][time_statistics][season_region][ensemble_statistics][warming_level]
for var in list(stats.keys()) :
if type(stats[var]) is not dict :
stats.pop(var)
else :
for time_stat in list(stats[var].keys()) :
if time_stat not in ["mean","std"] :
stats[var].pop(time_stat)
else :
for season in stats[var][time_stat] :
for estat in list(stats[var][time_stat][season].keys()):
if estat not in ["mean","nq5","nq95","ens"]:
stats[var][time_stat][season].pop(estat)
else:
if estat == "ens" : continue
# Recompute stat on filtered dic
values=ensemble_stat_series(stats,var,time_stat,season,estat)
thekeys=ensemble_stat_series(stats,var,time_stat,season,"x")
stats[var][time_stat][season][estat]={ k:v for (k,v) in zip(thekeys,values)}
for estat in [ "nq17","nq83"]:
values=ensemble_stat_series(stats,var,time_stat,season,estat)
thekeys=ensemble_stat_series(stats,var,time_stat,season,"x")
stats[var][time_stat][season][estat]={ k:v for (k,v) in zip(thekeys,values)}
#stats=stats['ssp585']
filen="%s/%s_%s%s.json"%(outdir,"Hydro_vars_change",data_versions_tag,output_version)
! mkdir -p {outdir}
with open(filen,"w") as f :
json.dump(stats,f,separators=(',', ': '),indent=1,ensure_ascii=True)
print("File written",filen)
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)