diff --git a/PE_CBL.py b/PE_CBL.py
index 72a584b0d782df10a44cd130c67f37e186b84dfa..8f89ef585e53134641a3339523fd67cf4b8a6b5f 100644
--- a/PE_CBL.py
+++ b/PE_CBL.py
@@ -4,446 +4,135 @@ import numpy as np
from matplotlib import pyplot as p
import pickle
import os
-from copy import deepcopy
+import json
# Own modules
from models import CBL
-from ENDA import experiment, transform_log, transform_none, transform_pit
+from ENDA import experiment
from graphics import *
+def load_or_run(settings):
+
+ try:
+ exp = pickle.load(open(settings['path']+'/'+settings['filename'], "rb"))
+ print('loaded experiment %s'%(settings['path']+'/'+settings['filename']))
+ except:
+ print('loading %s failed, running it'%(settings['path']+'/'+settings['filename']))
+ exp = experiment(settings)
+ pickle.dump(exp, open(settings['path']+'/'+settings['filename'], 'wb'))
+
+ return exp
+
if __name__ == '__main__':
# Global settings
- nobs = 16
- nens = 200
+ nens = 20
sigma_b_init = 0.1 # increased by a factor 10 in experiment B
- sigma_o_gn = 0.0 # set equal to sigma_o_as in experiments C and F
+ sigma_o_gn = 0.1 # set equal to sigma_o_as in experiments C and F
sigma_o_as = 0.1 # increased by a factor 10 in experiments C and F
# Define seed for random number generation, to ensure reproducibility
rnseed = 181612
- # Default CBL model settings
- default_cbl_settings ={
- # Physical parameters
- 'g' : 9.806,
- 'f' : 1e-4,
- 'kvonk' : 0.4,
- 'z0' : 0.1,
- 'ug' : 0,
- 'vg' : 10,
- 'theta_0' : 290,
- 'gamma' : 3e-3,
- 'Hmax' : 0.1,
- 'H0_perturbation_ampl_init' : 0.0,
- 'H0_perturbation_ampl_time' : 0.0,
- 'exct' : 0.3,
- 'pfac' : 1.5,
- 'Kmin' : 0.1,
- 'Kmax' : 200,
- # Model formulation
- 'is_bwind' : False,
- 'is_cgrad' : True,
- 'is_cycle' : False,
- # Domain parameters (Kmax determines dt)
- 'dz' : 50,
- 'ztop' : 4000,
- # Ensemble generation parameters (only effective for ensemble runs)
- # Part 1: perturbation of initial state
- # perturbation_type can be "smooth", "random", "uniform"
- # the adjective refers to vertical variability
- 'rnseed' : rnseed,
- 'perturb_ensemble_state' : True,
- 'perturbations_type' : "smooth",
- 'perturbations_theta_amplitude' : sigma_b_init,
- 'perturbations_uv_amplitude' : 0.1,
- 'perturbations_smooth_number' : 11,
- 'perturbations_symmetric' : True,
- 'simulate_error_growth' : False,
- 'error_growth_perturbations_amplitude' : 0.0,
- # Part 2: perturbation of parameters
- 'perturb_ensemble_parameters' : True,
- 'parameter_number' : 1,
- 'parameter_transform' : [ transform_pit ],
- 'parameter_ensemble_min' : np.array([0.5]),
- 'parameter_ensemble_max' : np.array([4.5]),
- 'parameter_true' : np.array([1.5]),
- # In case of parameter estimation
- 'do_parameter_estimation' : True,
- 'parameter_inflation_rtps_alpha' : np.array([0.8]),
- 'return_covariances_increments_and_innovations' : True
- }
- integration_dt = 0.25*default_cbl_settings["dz"]**2/default_cbl_settings["Kmax"]
-
- # Default assimilation settings
- default_da_settings = {'cbl_settings' : default_cbl_settings,
- 'tspinup' : 3600,
- 'trun' : 21600,
- 'assimilation_interval' : 300,
- 'randomize_obs': True,
- 'obs_coordinates' : np.linspace(0,2000,nobs),
- 'obs_kinds' : np.array(['theta']*nobs),
- 'obs_error_sdev_generate' : np.ones(nobs)*0,
- 'obs_error_sdev_assimilate' : np.ones(nobs)*sigma_o_as,
- 'localization_cutoff' : 400,
- 'nens' : nens,
- 'FILTER' : 'EAKF',
- 'inflation_rtps_alpha' : 0.2,
- 'rnseed' : rnseed,
- 'path' : "./runs/seed=%s_enssize=%s_EAKF_6hrs/"%(rnseed,nens)
- }
-
# Safety check on time steps
- assert np.mod(default_da_settings['assimilation_interval'],integration_dt)==0,\
- "Assimilation interval must be an integer multiplier of model dt"
+ #integration_dt = 0.25*default_cbl_settings["dz"]**2/default_cbl_settings["Kmax"]
+ #assert np.mod(default_da_settings['assimilation_interval'],integration_dt)==0,\
+ # "Assimilation interval must be an integer multiplier of model dt"
# Whether or not to run experiments without parameter estimation
- # Applies only to sets B&C; no-PE run is always computed for experiments A&D
+ # Applies only to sets B-C-D; no-PE run is always computed for experiments A
noPE_runs = False
- # Check if output path exists, if not, create it.
- if not os.path.exists(default_da_settings['path']):
- os.makedirs(default_da_settings['path'])
+ # Check if working directory exists, if not, create it.
+ workdir = "./runs/seed=181612_enssize=20_EAKF_6hrs/"
+ if not os.path.exists(workdir):
+ os.makedirs(workdir)
- ############################################################################
- # Experiment A (control)
- # Create a copy of the default settings
- cbl_settings_A = dict(default_cbl_settings)
- da_settings_A = dict(default_da_settings)
+ ##################################################################################
+ # Load the experiments or run them
- # Run it
- try:
- exp_A = pickle.load(open(default_da_settings['path']+"exp_A.pickle", "rb"))
- except:
- exp_A = experiment(da_settings_A)
- setattr(exp_A,'label','A')
- pickle.dump(exp_A, open(default_da_settings['path']+'exp_A.pickle', 'wb'))
-
- # Store the nature run, because it is used also in other experiments/figures
- exp_A.nr.has_parameters_in_state = True
- nature_run = exp_A.nr
- truths = exp_A.truths
- observations = exp_A.observations
-
- # Corresponding experiment without parameter estimation
- cbl_settings_A_noPE = dict(cbl_settings_A)
- da_settings_A_noPE = dict(da_settings_A)
- da_settings_A_noPE['nr'] = deepcopy(nature_run)
- da_settings_A_noPE['truths'] = truths
- da_settings_A_noPE['observations'] = observations
- da_settings_A_noPE['nr'].do_parameter_estimation = False
- cbl_settings_A_noPE['do_parameter_estimation'] = False
- cbl_settings_A_noPE['initial_perturbed_parameters'] = exp_A.da.initial_perturbed_parameters
- da_settings_A_noPE['cbl_settings'] = cbl_settings_A_noPE
-
- # Run it
- try:
- exp_A_noPE = pickle.load(open(default_da_settings['path']+'exp_A_noPE.pickle', "rb"))
- except:
- exp_A_noPE = experiment(da_settings_A_noPE)
- pickle.dump(exp_A_noPE, open(default_da_settings['path']+'exp_A_noPE.pickle', 'wb'))
+ with open('./config_files/da_exp_A1.json', 'r') as fp:
+ da_settings_A1 = json.load(fp)
+ exp_A1 = load_or_run(da_settings_A1)
- ########################################################################
- # Experiment B1
- # Create a copy of the default settings
- # Then re-use the available nature run and derived information
- cbl_settings_B1 = dict(default_cbl_settings)
- da_settings_B1 = dict(default_da_settings)
- da_settings_B1['nr'] = deepcopy(nature_run)
- da_settings_B1['truths'] = truths
- da_settings_B1['observations'] = observations
-
- # Change settings as necessary
- cbl_settings_B1['perturbations_theta_amplitude'] = sigma_b_init*10
- cbl_settings_B1['initial_perturbed_parameters'] = exp_A.da.initial_perturbed_parameters
- da_settings_B1['cbl_settings'] = cbl_settings_B1
- da_settings_B1['obs_error_sdev_assimilate'] = np.ones(nobs)*sigma_o_as
-
- # Run and save to disk
- try:
- exp_B1 = pickle.load(open(default_da_settings['path']+'exp_B1.pickle', "rb"))
- except:
- exp_B1 = experiment(da_settings_B1)
- setattr(exp_B1,'label','B1')
- pickle.dump(exp_B1, open(default_da_settings['path']+'exp_B1.pickle', 'wb'))
+ with open('./config_files/da_exp_A2.json', 'r') as fp:
+ da_settings_A2 = json.load(fp)
+ exp_A2 = load_or_run(da_settings_A2)
- if noPE_runs:
- # Corresponding experiment without parameter estimation
- cbl_settings_B1_noPE = dict(cbl_settings_B1)
- da_settings_B1_noPE = dict(da_settings_B1)
- cbl_settings_B1_noPE['do_parameter_estimation'] = False
- da_settings_B1_noPE['nr'] = deepcopy(nature_run)
- da_settings_B1_noPE['nr'].do_parameter_estimation = False
- da_settings_B1_noPE['cbl_settings'] = cbl_settings_B1_noPE
-
- # Run and save to disk
- try:
- exp_B1_noPE = pickle.load(open(default_da_settings['path']+'exp_B1_noPE.pickle', "rb"))
- except:
- exp_B1_noPE = experiment(da_settings_B1_noPE)
- pickle.dump(exp_B1_noPE, open(default_da_settings['path']+'exp_B1_noPE.pickle', 'wb'))
+ with open('./config_files/da_exp_B1.json', 'r') as fp:
+ da_settings_B1 = json.load(fp)
+ exp_B1 = load_or_run(da_settings_B1)
- ########################################################################
- # Experiment B2
- # Create a copy of the default settings
- # Then re-use the available nature run and derived information
- cbl_settings_B2 = dict(default_cbl_settings)
- da_settings_B2 = dict(default_da_settings)
- da_settings_B2['nr'] = deepcopy(nature_run)
- da_settings_B2['truths'] = truths
- da_settings_B2['observations'] = observations
-
- # Change settings as necessary
- cbl_settings_B2['perturbations_theta_amplitude'] = sigma_b_init
- cbl_settings_B2['initial_perturbed_parameters'] = exp_A.da.initial_perturbed_parameters
- da_settings_B2['cbl_settings'] = cbl_settings_B2
- da_settings_B2['obs_error_sdev_assimilate'] = np.ones(nobs)*sigma_o_as*10
-
- # Run and save to disk
- try:
- exp_B2 = pickle.load(open(default_da_settings['path']+'exp_B2.pickle', "rb"))
- except:
- exp_B2 = experiment(da_settings_B2)
- setattr(exp_B2,'label','B2')
- pickle.dump(exp_B2, open(default_da_settings['path']+'exp_B2.pickle', 'wb'))
+ with open('./config_files/da_exp_B2.json', 'r') as fp:
+ da_settings_B2 = json.load(fp)
+ exp_B2 = load_or_run(da_settings_B2)
- if noPE_runs:
- # Corresponding experiment without parameter estimation
- cbl_settings_B2_noPE = dict(cbl_settings_B2)
- da_settings_B2_noPE = dict(da_settings_B2)
- cbl_settings_B2_noPE['do_parameter_estimation'] = False
- da_settings_B2_noPE['nr'] = deepcopy(nature_run)
- da_settings_B2_noPE['nr'].do_parameter_estimation = False
- da_settings_B2_noPE['cbl_settings'] = cbl_settings_B2_noPE
-
- # Run and save to disk
- try:
- exp_B2_noPE = pickle.load(open(default_da_settings['path']+'exp_B2_noPE.pickle', "rb"))
- except:
- exp_B2_noPE = experiment(da_settings_B2_noPE)
- pickle.dump(exp_B2_noPE, open(default_da_settings['path']+'exp_B2_noPE.pickle', 'wb'))
+ with open('./config_files/da_exp_B3.json', 'r') as fp:
+ da_settings_B3 = json.load(fp)
+ exp_B3 = load_or_run(da_settings_B3)
- ########################################################################
- # Experiment B3
- # Create a copy of the default settings
- # Then re-use the available nature run and derived information
- cbl_settings_B3 = dict(default_cbl_settings)
- da_settings_B3 = dict(default_da_settings)
- da_settings_B3['nr'] = deepcopy(nature_run)
- da_settings_B3['truths'] = truths
- da_settings_B3['observations'] = observations
-
- # Change settings as necessary
- cbl_settings_B3['perturbations_theta_amplitude'] = sigma_b_init*10
- cbl_settings_B3['initial_perturbed_parameters'] = exp_A.da.initial_perturbed_parameters
- da_settings_B3['cbl_settings'] = cbl_settings_B3
- da_settings_B3['obs_error_sdev_assimilate'] = np.ones(nobs)*sigma_o_as*10
-
- # Run and save to disk
- try:
- exp_B3 = pickle.load(open(default_da_settings['path']+'exp_B3.pickle', "rb"))
- except:
- exp_B3 = experiment(da_settings_B3)
- setattr(exp_B3,'label','B3')
- pickle.dump(exp_B3, open(default_da_settings['path']+'exp_B3.pickle', 'wb'))
+ with open('./config_files/da_exp_B4.json', 'r') as fp:
+ da_settings_B4 = json.load(fp)
+ exp_B4 = load_or_run(da_settings_B4)
- if noPE_runs:
- # Corresponding experiment without parameter estimation
- cbl_settings_B3_noPE = dict(cbl_settings_B3)
- da_settings_B3_noPE = dict(da_settings_B3)
- cbl_settings_B3_noPE['do_parameter_estimation'] = False
- da_settings_B3_noPE['nr'] = deepcopy(nature_run)
- da_settings_B3_noPE['nr'].do_parameter_estimation = False
- da_settings_B3_noPE['cbl_settings'] = cbl_settings_B3_noPE
-
- # Run and save to disk
- try:
- exp_B3_noPE = pickle.load(open(default_da_settings['path']+'exp_B3_noPE.pickle', "rb"))
- except:
- exp_B3_noPE = experiment(da_settings_B3_noPE)
- pickle.dump(exp_B3_noPE, open(default_da_settings['path']+'exp_B3_noPE.pickle', 'wb'))
+ with open('./config_files/da_exp_C1.json', 'r') as fp:
+ da_settings_C1 = json.load(fp)
+ exp_C1 = load_or_run(da_settings_C1)
- ########################################################################
- # Experiment B4
- # Create a copy of the default settings
- cbl_settings_B4 = dict(default_cbl_settings)
- da_settings_B4 = dict(default_da_settings)
-
- # Change settings as necessary
- # Changes include generation of observations, so the existing nature run
- # can't be reused.
- cbl_settings_B4['initial_perturbed_parameters'] = exp_A.da.initial_perturbed_parameters
- da_settings_B4['cbl_settings'] = cbl_settings_B4
- da_settings_B4['obs_error_sdev_generate'] = np.ones(nobs)*sigma_o_as*5
- da_settings_B4['obs_error_sdev_assimilate'] = np.ones(nobs)*sigma_o_as*10
-
- # Run and save to disk
- try:
- exp_B4 = pickle.load(open(default_da_settings['path']+'exp_B4.pickle', "rb"))
- except:
- exp_B4 = experiment(da_settings_B4)
- setattr(exp_B4,'label','B4')
- pickle.dump(exp_B4, open(default_da_settings['path']+'exp_B4.pickle', 'wb'))
+ with open('./config_files/da_exp_C2.json', 'r') as fp:
+ da_settings_C2 = json.load(fp)
+ exp_C2 = load_or_run(da_settings_C2)
- if noPE_runs:
- # Corresponding experiment without parameter estimation
- cbl_settings_B4_noPE = dict(cbl_settings_B4)
- da_settings_B4_noPE = dict(da_settings_B4)
- cbl_settings_B4_noPE['do_parameter_estimation'] = False
- da_settings_B4_noPE['cbl_settings'] = cbl_settings_B4_noPE
-
- # Run and save to disk
- try:
- exp_B4_noPE = pickle.load(open(default_da_settings['path']+'exp_B4_noPE.pickle', "rb"))
- except:
- exp_B4_noPE = experiment(da_settings_B4_noPE)
- pickle.dump(exp_B4_noPE, open(default_da_settings['path']+'exp_B4_noPE.pickle', 'wb'))
+ with open('./config_files/da_exp_D1.json', 'r') as fp:
+ da_settings_D1 = json.load(fp)
+ exp_D1 = load_or_run(da_settings_D1)
- ########################################################################
- # Experiment C1
- # Create a copy of the default settings
- cbl_settings_C1 = dict(default_cbl_settings)
- da_settings_C1 = dict(default_da_settings)
-
- # Re-use the available nature run and derived information
- da_settings_C1['nr'] = deepcopy(nature_run)
- da_settings_C1['truths'] = truths
- da_settings_C1['observations'] = observations
-
- # Change settings as necessary
- cbl_settings_C1['initial_perturbed_parameters'] = exp_A.da.initial_perturbed_parameters
- cbl_settings_C1['is_cgrad'] = False
- da_settings_C1['cbl_settings'] = cbl_settings_C1
-
- # Run and save to disk
- try:
- exp_C1 = pickle.load(open(default_da_settings['path']+'exp_C1.pickle', "rb"))
- except:
- exp_C1 = experiment(da_settings_C1)
- setattr(exp_C1,'label','C1')
- pickle.dump(exp_C1, open(default_da_settings['path']+'exp_C1.pickle', 'wb'))
+ with open('./config_files/da_exp_D2.json', 'r') as fp:
+ da_settings_D2 = json.load(fp)
+ exp_D2 = load_or_run(da_settings_D2)
if noPE_runs:
- # Corresponding experiment without parameter estimation
- cbl_settings_C1_noPE = dict(cbl_settings_C1)
- da_settings_C1_noPE = dict(da_settings_C1)
- cbl_settings_C1_noPE['do_parameter_estimation'] = False
- da_settings_C1_noPE['cbl_settings'] = cbl_settings_C1_noPE
-
- # Run and save to disk
- try:
- exp_C1_noPE = pickle.load(open(default_da_settings['path']+'exp_C1_noPE.pickle', "rb"))
- except:
- exp_C1_noPE = experiment(da_settings_C1_noPE)
- pickle.dump(exp_C1_noPE, open(default_da_settings['path']+'exp_C1_noPE.pickle', 'wb'))
+ with open('./config_files/da_exp_B1_noPE.json', 'r') as fp:
+ da_settings_B1_noPE = json.load(fp)
+ exp_B1_noPE = load_or_run(da_settings_B1_noPE)
########################################################################
- # Experiment C2
- # Create a copy of the default settings
- cbl_settings_C2 = dict(default_cbl_settings)
- da_settings_C2 = dict(default_da_settings)
-
- # Re-use the available nature run and derived information
- da_settings_C2['nr'] = deepcopy(nature_run)
- da_settings_C2['truths'] = truths
- da_settings_C2['observations'] = observations
-
- # Change settings as necessary
- cbl_settings_C2['initial_perturbed_parameters'] = exp_A.da.initial_perturbed_parameters
- cbl_settings_C2['Hmax'] = 0.15
- da_settings_C2['cbl_settings'] = cbl_settings_C2
-
- # Run and save to disk
- try:
- exp_C2 = pickle.load(open(default_da_settings['path']+'exp_C2.pickle', "rb"))
- except:
- exp_C2 = experiment(da_settings_C2)
- setattr(exp_C2,'label','C2')
- pickle.dump(exp_C2, open(default_da_settings['path']+'exp_C2.pickle', 'wb'))
-
- if noPE_runs:
- # Corresponding experiment without parameter estimation
- cbl_settings_C2_noPE = dict(cbl_settings_C2)
- da_settings_C2_noPE = dict(da_settings_C2)
- cbl_settings_C2_noPE['do_parameter_estimation'] = False
- da_settings_C2_noPE['cbl_settings'] = cbl_settings_C2_noPE
-
- # Run and save to disk
- try:
- exp_C2_noPE = pickle.load(open(default_da_settings['path']+'exp_C2_noPE.pickle', "rb"))
- except:
- exp_C2_noPE = experiment(da_settings_C2_noPE)
- pickle.dump(exp_C2_noPE, open(default_da_settings['path']+'exp_C2_noPE.pickle', 'wb'))
-
- ########################################################################
- # Experiment D
- # Create a copy of the default settings
- cbl_settings_D = dict(default_cbl_settings)
- da_settings_D = dict(default_da_settings)
-
- # Change settings as necessary
- # Changes include generation of observations, so the existing nature run
- # can't be reused.
- cbl_settings_D['initial_perturbed_parameters'] = exp_A.da.initial_perturbed_parameters
- cbl_settings_D['perturbations_theta_amplitude'] = sigma_b_init*10
- cbl_settings_D['Hmax'] = 0.15
- cbl_settings_D['is_cgrad'] = False
- cbl_settings_D['simulate_error_growth'] = True
- cbl_settings_D['error_growth_perturbations_amplitude'] = sigma_b_init*5
- da_settings_D['cbl_settings'] = cbl_settings_D
- da_settings_D['obs_error_sdev_generate'] = np.ones(nobs)*sigma_o_as*5
- da_settings_D['obs_error_sdev_assimilate'] = np.ones(nobs)*sigma_o_as*10
-
- # Run and save to disk
- try:
- exp_D = pickle.load(open(default_da_settings['path']+'exp_D.pickle', "rb"))
- except:
- exp_D = experiment(da_settings_D)
- setattr(exp_D,'label','D')
- pickle.dump(exp_D, open(default_da_settings['path']+'exp_D.pickle', 'wb'))
-
- # Experiment matching D, but without parameter estimation
- if noPE_runs:
- cbl_settings_D_noPE = dict(cbl_settings_D)
- da_settings_D_noPE = dict(da_settings_D)
- cbl_settings_D_noPE['do_parameter_estimation'] = False
- da_settings_D_noPE['cbl_settings'] = cbl_settings_D_noPE
-
- try:
- exp_D_noPE = pickle.load(open(default_da_settings['path']+'exp_D_noPE.pickle', "rb"))
- except:
- exp_D_noPE = experiment(da_settings_D_noPE)
- pickle.dump(exp_A_noPE, open(default_da_settings['path']+'exp_D_noPE.pickle', 'wb'))
-
- ########################################################################
-
- # Decide what figures to plot
- fig101 = False
- fig102 = False
- fig103 = False
- fig104 = False
- fig105 = False
- fig106 = False
- fig107 = False
- fig108 = False
- fig109 = False
+ # Manuscript figures
+ fig01 = False # CHECKED OK # Sensitivity to p
+ fig02 = False # CHECKED OK # Successful PE experiment
+ fig03 = False # CHECKED OK # DA diagnostics
+ fig04 = False # CHECKED OK # Experiments with varying weights
+ fig05 = False # CHECKED OK # Experiments with systematic model errors
+ fig06 = False # CHECKED OK # Analysis increments in model error experiments
+ fig07 = False # CHECKED OK # Experiment with error growth
+ fig08 = False # CHECKED OK # Comparison successful vs unsuccessful PE
+ fig09 = False # Comparison EAKF vs LETKF
- # Other possible, optional, plots
- opt01 = False # assimilation of a single observation
- opt02 = False # assimilation of profiles at two times
- opt03 = False # assimilation of real observations
- opt04 = True # assimilation of real observations, different spinup times
+ # Other figures
+ opt01 = False # CHECKED OK # assimilation of a single observation
+ opt02 = False # CHECKED OK # assimilation of profiles at two times
+ opt03 = False # CHECKED OK # NOT EXPECTED TO WORK # assimilation of real observations
+ opt04 = True # NOT EXPECTED TO WORK # assimilation of real observations, different spinup times
+ opt05 = False # CHECKED OK # Check impact of localization cutoff
+ opt06 = False # CHECKED OK # Check impact of RTPS variance inflation factor for parameters
+ opt07 = False # CHECKED OK # Check impact of RTPS variance inflation factor for state vector
+ opt08 = False # CHECKED OK # Check sensitivity to observation sorting
+ opt09 = False # CHECKED OK # Plot covariance matrix
+ opt10 = False # CHECKED OK # Consistency checks
+ opt11 = False # CHECKED OK # Consistency checks
+
+ # A test to check if results are exactly reproducible
+ check_reproducibility = False
- if fig101:
+ if fig01:
# Create a copy of the default settings
- cbl_settings = dict(default_cbl_settings)
+ # and update as necessary
+ with open('./config_files/default_cbl.json', 'r') as fp:
+ cbl_settings = json.load(fp)
maxtime = 21600
cbl_settings['maxtime'] = maxtime
- # Disable parameter estimation (not used here)
- cbl_settings['do_parameter_estimation'] = False
-
# Panel a) deterministic run using default settings
cbl_det = CBL(cbl_settings)
cbl_det.initialize(1)
@@ -466,17 +155,19 @@ if __name__ == '__main__':
# Panel c) spread induced by p
# Do a free ensemble run (ensemble size set expliclity)
- cbl_settings_free = dict(default_cbl_settings)
+ # Set size to 200 members
+ with open('./config_files/default_cbl.json', 'r') as fp:
+ cbl_settings_free = json.load(fp)
cbl_settings_free['maxtime'] = maxtime
cbl_settings_free['perturb_ensemble_state'] = False
cbl_free = CBL(cbl_settings_free)
- cbl_free.initialize(nens)
+ cbl_free.initialize(200)
cbl_free.run(output_full_history=True)
# Make plots
ncont = 13
- fig, [[ax4, ax2],[ax1, ax3]] = p.subplots(2,2,constrained_layout=True)
- fig.set_size_inches(6,6)
+ fig, [[ax4, ax2],[ax1, ax3], [ax5, ax6] ] = p.subplots(3,2,constrained_layout=True)
+ fig.set_size_inches(6,9)
c1 = ax1.pcolormesh(cbl_det.history['time']/3600,cbl_det.zt,cbl_det.history['theta'],vmin=290,vmax=296)
ax1.set_ylim([0,zmax])
@@ -499,7 +190,7 @@ if __name__ == '__main__':
ax2.legend(loc=4,frameon=False)
ax2.set_title(r'b) $\overline{\theta}$ sensitivity to $p$')
- ax3,c3 = plot_spread(cbl_free,ax=ax3)
+ ax3,c3 = plot_free_ensemble(cbl_free,show='spread',ax=ax3)
ax3.set_ylabel(r'Height (m)')
ax3.set_title(r'd) $\sigma_\theta$ (K)')
ax3.set_xlabel('Time (h)')
@@ -516,85 +207,101 @@ if __name__ == '__main__':
ax4.set_xlim([0,0.5])
ax4.legend(loc=4,frameon=False)
+ ax5,c5 = plot_free_ensemble(cbl_free,show='cov',ax=ax5,cmap='RdBu_r')
+ ax5.set_ylabel(r'Height (m)')
+ ax5.set_title(r'e) cov$(p,\overline{\theta})$ (K)')
+ ax5.set_xlabel('Time (h)')
+ ax5.set_xticks(np.arange((maxtime+1)/3600))
+ p.colorbar(c5,orientation='horizontal')
+
+ ax6 = plot_linearity(cbl_free,ax=ax6)
+ ax6.set_ylabel(r'$\overline{\theta}$ (K)')
+ ax6.set_title(r'f) $p$ vs $\overline{\theta}$, $t=6$ h')
+ ax6.set_xlabel(r'$p$')
+
#p.setp(ax2.get_yticklabels(), visible=False)
#p.setp(ax3.get_yticklabels(), visible=False)
fig.savefig('fig01.png',format='png',dpi=300)
p.close(fig)
- if fig102:
-
- exp_A = pickle.load(open(default_da_settings['path']+"exp_A.pickle", "rb"))
+ if fig02:
- fig, [[ax0, ax1], [ax2, ax3]] = p.subplots(2,2,constrained_layout=True)
- fig.set_size_inches(6,6)
- #
- [ax0,ax1,ax2],c0,c1,c2 = plot_CBL_identifiability(exp_A,da_settings_A['obs_error_sdev_assimilate'][0],ax=[ax0,ax1,ax2])
- ax0.set_title(r'a) Exp. A$_1$, $\rho(p\prime\prime,y_b}$)')
- ax0.set_xlabel('Time (h)')
- ax0.set_ylabel('Height (m)')
- ax1.set_title(r'b) Exp. A$_1$, $\delta\overline{y}\cdot(\sigma_{p\prime\prime}/\sigma_{y^b})$')
- ax1.set_xlabel('Time (h)')
- ax1.set_ylabel('Height (m)')
- ax2.set_title(r'c) Exp. A$_1$, $\delta\overline{p}\prime\prime$')
- ax2.set_xlabel('Time (h)')
- ax2.set_ylabel('Height (m)')
- ax3 = plot_CBL_PE(exp_A,None,ax=ax3)
- ax3.set_title(r'd) Exp. A$_1$, evolution of $p$')
- ax3.set_xlabel('Time (h)')
- ax3.set_yticks([0,1,2,3,4,5])
- maxtime = exp_A.trun
- ax0.set_xlim([0,maxtime/3600])
- ax1.set_xlim([0,maxtime/3600])
- ax2.set_xlim([0,maxtime/3600])
- ax3.set_xlim([0,maxtime/3600])
- ax0.set_xticks(np.arange((maxtime+1)/3600))
- ax1.set_xticks(np.arange((maxtime+1)/3600))
- ax2.set_xticks(np.arange((maxtime+1)/3600))
- ax3.set_xticks(np.arange((maxtime+1)/3600))
- p.colorbar(c0,orientation='horizontal')
- p.colorbar(c1,orientation='horizontal')
- p.colorbar(c2,orientation='horizontal')
- #
- fig.savefig('fig02.png',format='png',dpi=300)
- p.close(fig)
+ exp_A1 = pickle.load(open(workdir+"exp_A1.pickle", "rb"))
+ plot_overview(exp_A1,"A$_1$",'fig02.png')
- if fig103:
+ if fig03:
- exp_A = pickle.load(open(default_da_settings['path']+"exp_A.pickle", "rb"))
- exp_A_noPE = pickle.load(open(default_da_settings['path']+"exp_A_noPE.pickle", "rb"))
+ exp1 = pickle.load(open(workdir+"exp_A1.pickle", "rb"))
+ exp2 = pickle.load(open(workdir+"exp_A2.pickle", "rb"))
- experiments_pe = [exp_A]
- experiments_nope = [exp_A_noPE]
- labels = ["Exp. A"]
+ fig, [[ax1, ax2], [ax3, ax4]] = p.subplots(2,2,constrained_layout=True)
+ fig.set_size_inches(6,4)
+ plot_diagnostics(exp1,show='errors',axes=[ax1,ax2,ax3,ax4],label='with PE',linecolor='#1f77b4')
+ plot_diagnostics(exp2,show='errors',axes=[ax1,ax2,ax3,ax4],label='without PE',linecolor='#9EC9E9')
+ ax1.set_title('a) Analysis error')
+ ax1.set_xlabel(r'RMSE$^a_\theta$')
+ ax2.set_title('b) First-guess error')
+ ax2.set_xlabel(r'RMSE$^b_\theta$')
+ ax3.set_title('c) Error reduction')
+ ax3.set_xlabel(r'RMSE$^b_\theta-$RMSE$^a_\theta$')
+ ax4.set_title('d) Spread-error consistency')
+ ax4.set_xlabel(r'$\sigma^b_\theta$/RMSE$^b_\theta$')
+ ax1.set_ylabel('height (m)')
+ ax3.set_ylabel('height (m)')
+ #ax1.set_xlim([0,0.2])
+ #ax2.set_xlim([0,0.2])
+ #ax3.set_xlim([0,0.02])
+ #ax4.set_xlim([0,10])
+ fig.savefig('fig03_with_errors.png',format='png',dpi=300)
+ p.close(fig)
- plot_diagnostics(experiments_pe,experiments_nope,labels,'fig03.png')
+ fig, [[ax1, ax2], [ax3, ax4]] = p.subplots(2,2,constrained_layout=True)
+ fig.set_size_inches(6,4)
+ plot_diagnostics(exp1,show='deviations',axes=[ax1,ax2,ax3,ax4],label='with PE',linecolor=u'#1f77b4')
+ plot_diagnostics(exp2,show='deviations',axes=[ax1,ax2,ax3,ax4],label='without PE',linecolor=u'#9EC9E9')
+ ax1.set_title('a) Analysis error')
+ ax1.set_xlabel(r'RMSD$^a_\theta$')
+ ax2.set_title('b) First-guess error')
+ ax2.set_xlabel(r'RMSD$^b_\theta$')
+ ax3.set_title('c) Error reduction')
+ ax3.set_xlabel(r'RMSD$^b_\theta-$RMSD$^a_\theta$')
+ ax4.set_title('d) Spread-error consistency')
+ ax4.set_xlabel(r'$\sigma^b_\theta$/RMSD$^b_\theta$')
+ ax1.set_ylabel('height (m)')
+ ax3.set_ylabel('height (m)')
+ #ax1.set_xlim([0,0.2])
+ #ax2.set_xlim([0,0.2])
+ #ax3.set_xlim([0,0.02])
+ #ax4.set_xlim([0,10])
+ fig.savefig('fig03_with_deviations.png',format='png',dpi=300)
+ p.close(fig)
- if fig104:
+ if fig04:
- exp_B1 = pickle.load(open(default_da_settings['path']+"exp_B1.pickle", "rb"))
- exp_B2 = pickle.load(open(default_da_settings['path']+"exp_B2.pickle", "rb"))
- exp_B3 = pickle.load(open(default_da_settings['path']+"exp_B3.pickle", "rb"))
- exp_B4 = pickle.load(open(default_da_settings['path']+"exp_B4.pickle", "rb"))
+ exp_B1 = pickle.load(open(workdir+"exp_B1.pickle", "rb"))
+ exp_B2 = pickle.load(open(workdir+"exp_B2.pickle", "rb"))
+ exp_B3 = pickle.load(open(workdir+"exp_B3.pickle", "rb"))
+ exp_B4 = pickle.load(open(workdir+"exp_B4.pickle", "rb"))
fig, [[ax1, ax2], [ax3, ax4]] = p.subplots(2,2,constrained_layout=True)
fig.set_size_inches(6,4)
#
- ax1 = plot_CBL_PE(exp_B1,None,ax=ax1)
+ ax1 = plot_CBL_PE(exp_B1,None,ax=ax1,show_true_parameter=True)
ax1.set_title(r'a) Exp. B$_1$')
ax1.set_xlabel('Time (h)')
ax1.set_ylabel(r'$p$')
#
- ax2 = plot_CBL_PE(exp_B2,None,ax=ax2)
+ ax2 = plot_CBL_PE(exp_B2,None,ax=ax2,show_true_parameter=True)
ax2.set_title(r'b) Exp. B$_2$')
ax2.set_xlabel('Time (h)')
ax2.set_ylabel(r'$p$')
#
- ax3 = plot_CBL_PE(exp_B3,None,ax=ax3)
+ ax3 = plot_CBL_PE(exp_B3,None,ax=ax3,show_true_parameter=True)
ax3.set_title(r'c) Exp. B$_3$')
ax3.set_xlabel('Time (h)')
ax3.set_ylabel(r'$p$')
#
- ax4 = plot_CBL_PE(exp_B4,None,ax=ax4)
+ ax4 = plot_CBL_PE(exp_B4,None,ax=ax4,show_true_parameter=True)
ax4.set_title(r'd) Exp. B$_4$')
ax4.set_xlabel('Time (h)')
ax4.set_ylabel(r'$p$')
@@ -607,20 +314,20 @@ if __name__ == '__main__':
fig.savefig('fig04.png',format='png',dpi=300)
p.close(fig)
- if fig105:
+ if fig05:
- exp_C1 = pickle.load(open(default_da_settings['path']+"exp_C1.pickle", "rb"))
- exp_C2 = pickle.load(open(default_da_settings['path']+"exp_C2.pickle", "rb"))
+ exp_C1 = pickle.load(open(workdir+"exp_C1.pickle", "rb"))
+ exp_C2 = pickle.load(open(workdir+"exp_C2.pickle", "rb"))
fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
fig.set_size_inches(6,2)
#
- ax1 = plot_CBL_PE(exp_C1,None,ax=ax1)
+ ax1 = plot_CBL_PE(exp_C1,None,ax=ax1,show_true_parameter=True)
ax1.set_title(r'a) Exp. C$_1$')
ax1.set_xlabel('Time (h)')
ax1.set_ylabel(r'$p$')
#
- ax2 = plot_CBL_PE(exp_C2,None,ax=ax2)
+ ax2 = plot_CBL_PE(exp_C2,None,ax=ax2,show_true_parameter=True)
ax2.set_title(r'b) Exp. C$_2$')
ax2.set_xlabel('Time (h)')
ax2.set_ylabel(r'$p$')
@@ -631,10 +338,10 @@ if __name__ == '__main__':
fig.savefig('fig05.png',format='png',dpi=300)
p.close(fig)
- if fig106:
+ if fig06:
- exp_C1 = pickle.load(open(default_da_settings['path']+"exp_C1.pickle", "rb"))
- exp_C2 = pickle.load(open(default_da_settings['path']+"exp_C2.pickle", "rb"))
+ exp_C1 = pickle.load(open(workdir+"exp_C1.pickle", "rb"))
+ exp_C2 = pickle.load(open(workdir+"exp_C2.pickle", "rb"))
npar2 = exp_C1.nr.parameter_number
npar3 = exp_C2.nr.parameter_number
@@ -667,18 +374,18 @@ if __name__ == '__main__':
fig.savefig('fig06.png',format='png',dpi=300)
p.close(fig)
- if fig107:
+ if fig07:
- exp_D = pickle.load(open(default_da_settings['path']+"exp_D.pickle", "rb"))
+ exp_D = pickle.load(open(workdir+"exp_D2.pickle", "rb"))
fig, [ax0, ax1] = p.subplots(1,2,constrained_layout=True)
fig.set_size_inches(6,3)
#
- ax0,c0 = plot_CBL_identifiability(exp_D,da_settings_D['obs_error_sdev_assimilate'][0],ax=ax0)
+ ax0,c0 = plot_CBL_identifiability(exp_D,da_settings_D2['obs_error_sdev_assimilate'],ax=ax0)
ax0.set_title(r'a) Exp. D, $\rho(p\prime\prime,y_b}$)')
ax0.set_xlabel('Time (h)')
ax0.set_ylabel('Height (m)')
- ax1 = plot_CBL_PE(exp_D,None,ax=ax1)
+ ax1 = plot_CBL_PE(exp_D,None,ax=ax1,show_true_parameter=True)
ax1.set_title(r'b) Exp. D, evolution of $p$')
ax1.set_xlabel('Time (h)')
ax1.set_yticks([0,1,2,3,4,5])
@@ -692,15 +399,15 @@ if __name__ == '__main__':
fig.savefig('fig07.png',format='png',dpi=300)
p.close(fig)
- if fig108:
+ if fig08:
- exp_A = pickle.load(open(default_da_settings['path']+"exp_A.pickle", "rb"))
- exp_D = pickle.load(open(default_da_settings['path']+"exp_D.pickle", "rb"))
+ exp_A = pickle.load(open(workdir+"exp_A1.pickle", "rb"))
+ exp_D = pickle.load(open(workdir+"exp_D2.pickle", "rb"))
fig, [[ax0, ax1],[ax2,ax3]] = p.subplots(2,2,constrained_layout=True)
fig.set_size_inches(6,6)
#
- [ax1,ax0],c0,c1 = plot_CBL_identifiability(exp_A,da_settings_A['obs_error_sdev_assimilate'][0],ax=[ax1,ax0])
+ [ax1,ax0],c0,c1 = plot_CBL_identifiability(exp_A,da_settings_A1['obs_error_sdev_assimilate'],ax=[ax1,ax0])
ax0.set_title(r'a) Exp. A$_1$, $\delta\overline{y}$ (K)')
ax0.set_xlabel('Time (h)')
ax0.set_ylabel('Height (m)')
@@ -712,11 +419,11 @@ if __name__ == '__main__':
ax1.set_xlim([0,maxtime/3600])
ax0.set_xticks(np.arange((maxtime+1)/3600))
ax1.set_xticks(np.arange((maxtime+1)/3600))
- [ax3,ax2],c2,c3 = plot_CBL_identifiability(exp_D,da_settings_D['obs_error_sdev_assimilate'][0],ax=[ax3,ax2])
- ax2.set_title(r'c) Exp. D, $\delta\overline{y}$ (K)')
+ [ax3,ax2],c2,c3 = plot_CBL_identifiability(exp_D,da_settings_D2['obs_error_sdev_assimilate'],ax=[ax3,ax2])
+ ax2.set_title(r'c) Exp. D$_2$, $\delta\overline{y}$ (K)')
ax2.set_xlabel('Time (h)')
ax2.set_ylabel('Height (m)')
- ax3.set_title(r'd) Exp. D, $\sigma_{p\prime\prime}/\sigma_{y^b}$ (K$^{-1}$)')
+ ax3.set_title(r'd) Exp. D$_2$, $\sigma_{p\prime\prime}/\sigma_{y^b}$ (K$^{-1}$)')
ax3.set_xlabel('Time (h)')
ax3.set_ylabel('Height (m)')
maxtime = exp_D.trun
@@ -732,36 +439,40 @@ if __name__ == '__main__':
fig.savefig('fig08.png',format='png',dpi=300)
p.close(fig)
- if fig109:
+ if fig09:
- path_EAKF = "./runs/seed=181612_enssize=200_EAKF_6hrs/"
- path_LETKF = "./runs/seed=181612_enssize=200_LETKF_6hrs/"
+ path_EAKF = ""
+ path_LETKF = "./runs/seed=181612_enssize=20_LETKF_6hrs/"
- exp_Aa = pickle.load(open(path_EAKF+"exp_A.pickle", "rb"))
+ exp_Aa = pickle.load(open(path_EAKF+"exp_A1.pickle", "rb"))
exp_B1a = pickle.load(open(path_EAKF+"exp_B1.pickle", "rb"))
exp_B2a = pickle.load(open(path_EAKF+"exp_B2.pickle", "rb"))
exp_B3a = pickle.load(open(path_EAKF+"exp_B3.pickle", "rb"))
exp_B4a = pickle.load(open(path_EAKF+"exp_B4.pickle", "rb"))
exp_C1a = pickle.load(open(path_EAKF+"exp_C1.pickle", "rb"))
exp_C2a = pickle.load(open(path_EAKF+"exp_C2.pickle", "rb"))
- exp_Da = pickle.load(open(path_EAKF+"exp_D.pickle", "rb"))
+ exp_D1a = pickle.load(open(path_EAKF+"exp_D1.pickle", "rb"))
+ exp_D2a = pickle.load(open(path_EAKF+"exp_D2.pickle", "rb"))
- exp_Ab = pickle.load(open(path_LETKF+"exp_A.pickle", "rb"))
+ exp_Ab = pickle.load(open(path_LETKF+"exp_A1.pickle", "rb"))
exp_B1b = pickle.load(open(path_LETKF+"exp_B1.pickle", "rb"))
exp_B2b= pickle.load(open(path_LETKF+"exp_B2.pickle", "rb"))
exp_B3b = pickle.load(open(path_LETKF+"exp_B3.pickle", "rb"))
exp_B4b = pickle.load(open(path_LETKF+"exp_B4.pickle", "rb"))
exp_C1b = pickle.load(open(path_LETKF+"exp_C1.pickle", "rb"))
exp_C2b = pickle.load(open(path_LETKF+"exp_C2.pickle", "rb"))
- exp_Db = pickle.load(open(path_LETKF+"exp_D.pickle", "rb"))
+ exp_D1b = pickle.load(open(path_LETKF+"exp_D1.pickle", "rb"))
+ exp_D2b = pickle.load(open(path_LETKF+"exp_D2.pickle", "rb"))
fig, [ax1, ax2] = p.subplots(1,2)
fig.set_size_inches(6,3)
#
- cc = 0
- linecolors = ['darkviolet','chocolate','darkorange','orange','gold','darkgreen','seagreen','royalblue']
- labels = [r'Exp. A',r'Exp. B$_1$',r'Exp. B$_2$',r'Exp. B$_3$',r'Exp. B$_4$',r'Exp. C$_1$',r'Exp. C$_2$',r'Exp. D']
- for exp in [exp_Aa,exp_B1a,exp_B2a,exp_B3a,exp_B4a,exp_C1a,exp_C2a,exp_Da]:
+ linecolors = ['darkviolet','chocolate','darkorange','orange','gold','darkgreen','seagreen','royalblue','black']
+ labels = [r'Exp. A',r'Exp. B$_1$',r'Exp. B$_2$',r'Exp. B$_3$',r'Exp. B$_4$',r'Exp. C$_1$',r'Exp. C$_2$',r'Exp. D$_1$',r'Exp. D$_2$']
+ experiments1 = [exp_Aa,exp_B1a,exp_B2a,exp_B3a,exp_B4a,exp_C1a,exp_C2a,exp_D1a,exp_D2a]
+ experiments2 = [exp_Ab,exp_B1b,exp_B2b,exp_B3b,exp_B4b,exp_C1b,exp_C2b,exp_D1b,exp_D2b]
+ #
+ for exp,label,linecolor in zip(experiments1,labels,linecolors):
# Experiment settings
true_value = exp.nr.pfac
parameter_number = exp.nr.parameter_number
@@ -773,14 +484,12 @@ if __name__ == '__main__':
t = exp.da.time/3600.
# Initial parameter distribution
- ax1.step([0,t[0]],np.mean(initpar_phys,axis=1)*np.array([1,1]),color=linecolors[cc],label=labels[cc])
+ ax1.step([0,t[0]],np.mean(initpar_phys,axis=1)*np.array([1,1]),color=linecolor,label=label)
# Later times
- ax1.step(t,np.mean(par_phys,axis=1),color=linecolors[cc])
- cc += 1
+ ax1.step(t,np.mean(par_phys,axis=1),color=linecolor)
- cc = 0
- for exp in [exp_Ab,exp_B1b,exp_B2b,exp_B3b,exp_B4b,exp_C1b,exp_C2b,exp_Db]:
+ for exp,label,linecolor in zip(experiments2,labels,linecolors):
# Experiment settings
true_value = exp.nr.pfac
parameter_number = exp.nr.parameter_number
@@ -792,11 +501,10 @@ if __name__ == '__main__':
t = exp.da.time/3600.
# Initial parameter distribution
- ax2.step([0,t[0]],np.mean(initpar_phys,axis=1)*np.array([1,1]),color=linecolors[cc])
+ ax2.step([0,t[0]],np.mean(initpar_phys,axis=1)*np.array([1,1]),color=linecolor)
# Later times
- ax2.step(t,np.mean(par_phys,axis=1),color=linecolors[cc])
- cc += 1
+ ax2.step(t,np.mean(par_phys,axis=1),color=linecolor)
# Other stuff
ax1.axhline(y=true_value,linestyle='--',color='black')
@@ -822,351 +530,290 @@ if __name__ == '__main__':
if opt01:
- da_settings = {'cbl_settings' : dict(default_cbl_settings),
- 'tspinup' : 10800,
- 'trun' : 3600,
- 'assimilation_interval' : 3600,
- 'obs_coordinates' : np.array([433.]),
- 'obs_kinds' : np.array(['theta']),
- 'obs_error_sdev_generate' : np.array([.2]),
- 'obs_error_sdev_assimilate' : np.array([1.]),
- 'localization_cutoff' : 400,
- 'nens' : 50,
- 'FILTER' : 'EAKF',
- 'inflation_rtps_alpha' : 0.1,
- }
- exp = experiment(da_settings)
+ with open('./default_da.json', 'r') as fp:
+ da_settings = json.load(fp)
+ exp = load_or_run(da_settings)
+
plot_CBL_assimilation(exp,'single_assimilation.png')
if opt02:
- # Make plots
- fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
- fig.set_size_inches(6,3)
+ plot_profiles(exp_A1,da_settings_A1,'profiles.png')
- which_cycle = 0
- tt1 = (which_cycle+1)*da_settings_A['assimilation_interval']
- ax1 = plot_CBL_assimilation(exp_A,None,which_cycle=which_cycle,ax=ax1)
- ax1.set_ylim([0,2000])
- ax1.set_xlim([288,298])
- ax1.set_ylabel(r'height (m)')
- ax1.set_title(r'a) Analysis at $t=$%u s'%tt1)
- ax1.legend(frameon=False)
-
- which_cycle = 31
- tt2 = (which_cycle+1)*da_settings_A['assimilation_interval']
- ax2 = plot_CBL_assimilation(exp_A,None,which_cycle=which_cycle,ax=ax2)
- ax2.set_ylim([0,2000])
- ax2.set_xlim([288,298])
- ax2.set_ylabel(r'height (m)')
- ax2.set_title(r'b) Analysis at $t=$%u s'%tt2)
+ if opt03:
- fig.savefig('profiles.png',format='png',dpi=300)
- p.close(fig)
+ # Control experiment (an OSSE with nature run)
+ with open('./runs/real/da_exp_RC.json', 'r') as fp:
+ da_settings_RC = json.load(fp)
+ exp_RC = load_or_run(da_settings_RC)
- if opt03:
# Test with "real" observations (no nature run)
-
- nens = 20
- ############################################################################
- # Control experiment (perfect-model)
- # Create a copy of the default settings
- # and change some of the settings as needed
- cbl_settings_RC = dict(default_cbl_settings)
- cbl_settings_RC['theta_0'] = 290
- cbl_settings_RC['gamma'] = 3e-3
- cbl_settings_RC['Hmax'] = 0.12
- #cbl_settings_RC['parameter_ensemble_min'] = np.array([0.5])
- #cbl_settings_RC['parameter_ensemble_max'] = np.array([2.5])
- da_settings_RC = dict(default_da_settings)
- da_settings_RC['cbl_settings'] = dict(cbl_settings_RC)
- da_settings_RC['nens'] = nens
- da_settings_RC['type'] = 'OSSE'
- da_settings_RC['path'] = "./runs/real_init_exp_seed=%s_enssize=%s_EAKF_6hrs/"%(rnseed,nens)
-
- # Run and save to disk
- try:
- exp_RC = pickle.load(open(da_settings_RC['path']+"exp_RC.pickle", "rb"))
- except:
- exp_RC = experiment(da_settings_RC)
- setattr(exp_RC,'label','RC')
- pickle.dump(exp_RC, open(da_settings_RC['path']+'exp_RC.pickle', 'wb'))
-
- ############################################################################
- # Imperfect-model experiment
- # Create a copy of the default settings
- # and change some of the settings as needed
- cbl_settings_R = dict(default_cbl_settings)
- cbl_settings_R['theta_0'] = 290
- cbl_settings_R['gamma'] = 3e-3
- cbl_settings_R['Hmax'] = 0.12
- #cbl_settings_R['parameter_ensemble_min'] = np.array([0.5])
- #cbl_settings_R['parameter_ensemble_max'] = np.array([2.5])
- da_settings_R = dict(default_da_settings)
- da_settings_R['cbl_settings'] = dict(cbl_settings_R)
- da_settings_R['nens'] = nens
- da_settings_R['type'] = 'real'
- da_settings_R['observation_files'] = './LES_data/Theta/*csv'
- da_settings_R['obs_coordinates'] = None
- da_settings_R['obs_error_sdev_generate'] = None
- da_settings_R['path'] = "./runs/real_init_exp_seed=%s_enssize=%s_EAKF_6hrs/"%(rnseed,nens)
-
- # Run and save to disk
- try:
- exp_R = pickle.load(open(da_settings_R['path']+'exp_R.pickle', "rb"))
- except:
- exp_R = experiment(da_settings_R)
- setattr(exp_R,'label','R')
- pickle.dump(exp_R, open(da_settings_R['path']+'exp_R.pickle', 'wb'))
-
- ############################################################################
- # Run a free ensemble if desired
- freens = False
- if freens:
- maxtime = 3600+21600 # spinup + analysis run
- # Do a free ensemble run (ensemble size set expliclity)
- cbl_settings_free = dict(cbl_settings_R)
- cbl_settings_free['maxtime'] = maxtime
- cbl_settings_free['perturb_ensemble_state'] = False
- cbl_settings_free['parameter_ensemble_min'] = np.array([0.5])
- cbl_settings_free['parameter_ensemble_max'] = np.array([2.5])
- cbl_free = CBL(cbl_settings_free)
- cbl_free.initialize(nens)
- cbl_free.run(output_full_history=True)
- # Make a plot about the free ensemble
- fig, ax1 = p.subplots(1,1,constrained_layout=True)
- fig.set_size_inches(4,3)
- ax1,c1 = plot_spread(cbl_free,plot='mean',ax=ax1)
- ax1.set_ylabel(r'Height (m)')
- ax1.set_title(r'$\overline{\theta}$ (K)')
- ax1.set_xlabel('Time (h)')
- ax1.set_xticks(np.arange((maxtime+1)/3600))
- c1.set_clim([290,296])
- p.colorbar(c1,orientation='horizontal')
- fig.savefig('expR_free_ens.png',format='png',dpi=300)
- p.close(fig)
+ with open('./runs/real/da_exp_R.json', 'r') as fp:
+ da_settings_R = json.load(fp)
+ exp_R = load_or_run(da_settings_R)
#########################################################################################à
- # Make plots: perfect-model experiment (RC)
- fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
- fig.set_size_inches(6,3)
+ # Plots
+ for exp, label in zip([exp_RC, exp_R],['RC','R']):
+ plot_overview(exp,label,'exp_%s.png'%label)
- which_cycle = 0
- tt1 = (which_cycle+1)*da_settings_RC['assimilation_interval']
- ax1 = plot_CBL_assimilation(exp_RC,None,which_cycle=which_cycle,ax=ax1)
- ax1.set_ylim([0,2000])
- ax1.set_xlim([288,298])
- ax1.set_ylabel(r'height (m)')
- ax1.set_title(r'a) Analysis at $t=$%u s'%tt1)
- ax1.legend(frameon=False)
-
- which_cycle = 31
- tt2 = (which_cycle+1)*da_settings_RC['assimilation_interval']
- ax2 = plot_CBL_assimilation(exp_RC,None,which_cycle=which_cycle,ax=ax2)
- ax2.set_ylim([0,2000])
- ax2.set_xlim([288,298])
- ax2.set_ylabel(r'height (m)')
- ax2.set_title(r'b) Analysis at $t=$%u s'%tt2)
+ plot_profiles(exp_RC,da_settings_RC,'expRC_profiles.png')
+ plot_profiles(exp_R,da_settings_R,'expR_profiles.png')
- fig.savefig('expRC_profiles.png',format='png',dpi=300)
- p.close(fig)
+ if opt04:
- fig, [[ax0, ax1], [ax2, ax3]] = p.subplots(2,2,constrained_layout=True)
- fig.set_size_inches(6,6)
- #
- [ax0,ax1,ax2],c0,c1,c2 = plot_CBL_identifiability(exp_RC,da_settings_RC['obs_error_sdev_assimilate'][0],ax=[ax0,ax1,ax2])
- ax0.set_title(r'a) Exp. RC, $\rho(p\prime\prime,y_b}$)')
- ax0.set_xlabel('Time (h)')
- ax0.set_ylabel('Height (m)')
- ax1.set_title(r'b) Exp. RC, $\delta\overline{y}\cdot(\sigma_{p\prime\prime}/\sigma_{y^b})$')
- ax1.set_xlabel('Time (h)')
- ax1.set_ylabel('Height (m)')
- ax2.set_title(r'c) Exp. RC, $\delta\overline{p}\prime\prime$')
- ax2.set_xlabel('Time (h)')
- ax2.set_ylabel('Height (m)')
- ax3 = plot_CBL_PE(exp_RC,None,ax=ax3)
- ax3.set_title(r'd) Exp. RC, evolution of $p$')
- ax3.set_xlabel('Time (h)')
- ax3.set_yticks([0,1,2,3,4,5])
- maxtime = exp_R.trun
- ax0.set_xlim([0,maxtime/3600])
- ax1.set_xlim([0,maxtime/3600])
- ax2.set_xlim([0,maxtime/3600])
- ax3.set_xlim([0,maxtime/3600])
- ax0.set_xticks(np.arange((maxtime+1)/3600))
- ax1.set_xticks(np.arange((maxtime+1)/3600))
- ax2.set_xticks(np.arange((maxtime+1)/3600))
- ax3.set_xticks(np.arange((maxtime+1)/3600))
+ # Test with "real" observations (no nature run)
+ with open('./runs/real/da_exp_R.json', 'r') as fp:
+ da_settings_S0 = json.load(fp)
+ exp_S0 = load_or_run(da_settings_S0)
+ exp_S0.label = 'S0'
+
+ # Additional tests with longer spinupt time
+ with open('./runs/real/da_exp_S1.json', 'r') as fp:
+ da_settings_S1 = json.load(fp)
+ exp_S1 = load_or_run(da_settings_S1)
+
+ with open('./runs/real/da_exp_S2.json', 'r') as fp:
+ da_settings_S2 = json.load(fp)
+ exp_S2 = load_or_run(da_settings_S2)
+
+ with open('./runs/real/da_exp_S3.json', 'r') as fp:
+ da_settings_S3 = json.load(fp)
+ exp_S3 = load_or_run(da_settings_S3)
+
+ for exp,label in zip([exp_S0,exp_S1,exp_S2,exp_S3],['S0','S1','S2','S3']):
+ plot_overview(exp,label,'exp_%s.png'%label)
+
+ if opt05:
+
+ with open('./config_files/da_exp_L1.json', 'r') as fp:
+ da_settings_localization_1 = json.load(fp)
+ exp_localization_1 = load_or_run(da_settings_localization_1)
+
+ with open('./config_files/da_exp_A1.json', 'r') as fp:
+ da_settings_localization_2 = json.load(fp)
+ exp_localization_2 = load_or_run(da_settings_localization_2)
+ exp_localization_2.label = 'L2'
+
+ with open('./config_files/da_exp_L3.json', 'r') as fp:
+ da_settings_localization_3 = json.load(fp)
+ exp_localization_3 = load_or_run(da_settings_localization_3)
+
+ with open('./config_files/da_exp_L4.json', 'r') as fp:
+ da_settings_localization_4 = json.load(fp)
+ exp_localization_4 = load_or_run(da_settings_localization_4)
+
+ for exp,settings,label in zip([exp_localization_1,exp_localization_2,exp_localization_3,exp_localization_4],\
+ [da_settings_localization_1,da_settings_localization_2,da_settings_localization_3,da_settings_localization_4],
+ [r"L$_1$",r"L$_2$",r"L$_3$",r"L$_4$"]):
+ print(exp.localization_cutoff)
+
+ plot_overview(exp,label,'exp_%s.png'%exp.label)
+
+ if opt06:
+
+ with open('./config_files/da_exp_I1.json', 'r') as fp:
+ da_settings_inflation_1 = json.load(fp)
+ exp_inflation_1 = load_or_run(da_settings_inflation_1)
+
+ with open('./config_files/da_exp_I2.json', 'r') as fp:
+ da_settings_inflation_2 = json.load(fp)
+ exp_inflation_2 = load_or_run(da_settings_inflation_2)
+
+ with open('./config_files/da_exp_A1.json', 'r') as fp:
+ da_settings_inflation_3 = json.load(fp)
+ exp_inflation_3 = load_or_run(da_settings_inflation_3)
+ exp_inflation_3.label = 'I3'
+
+ with open('./config_files/da_exp_I4.json', 'r') as fp:
+ da_settings_inflation_4 = json.load(fp)
+ exp_inflation_4 = load_or_run(da_settings_inflation_4)
+
+ for exp,settings,label in zip([exp_inflation_1,exp_inflation_2,exp_inflation_3,exp_inflation_4],\
+ [da_settings_inflation_1,da_settings_inflation_2,da_settings_inflation_3,da_settings_inflation_4],
+ [r"I$_1$",r"I$_2$",r"I$_3$",r"I$_4$"]):
+ print(exp.cbl_settings["parameter_inflation_rtps_alpha"])
+
+ plot_overview(exp,label,'exp_%s.png'%exp.label,show_spread=True,show_transformed_parameter=True)
+
+ if opt07:
+
+ with open('./config_files/da_exp_I5.json', 'r') as fp:
+ da_settings_inflation_5 = json.load(fp)
+ exp_inflation_5 = load_or_run(da_settings_inflation_5)
+
+ with open('./config_files/da_exp_A1.json', 'r') as fp:
+ da_settings_inflation_6 = json.load(fp)
+ exp_inflation_6 = load_or_run(da_settings_inflation_6)
+ exp_inflation_6.label = 'I6'
+
+ with open('./config_files/da_exp_I7.json', 'r') as fp:
+ da_settings_inflation_7 = json.load(fp)
+ exp_inflation_7 = load_or_run(da_settings_inflation_7)
+
+ with open('./config_files/da_exp_I8.json', 'r') as fp:
+ da_settings_inflation_8 = json.load(fp)
+ exp_inflation_8 = load_or_run(da_settings_inflation_8)
+
+ for exp,settings,label in zip([exp_inflation_5,exp_inflation_6,exp_inflation_7,exp_inflation_8],\
+ [da_settings_inflation_5,da_settings_inflation_6,da_settings_inflation_7,da_settings_inflation_8],
+ [r"I$_5$",r"I$_6$",r"I$_7$",r"I$_8$"]):
+ print(exp.inflation_rtps_alpha)
+
+ plot_overview(exp,label,'exp_%s.png'%exp.label)
+
+ if opt08:
+
+ with open('./config_files/da_exp_A1.json', 'r') as fp:
+ da_settings_O1 = json.load(fp)
+ exp_O1 = load_or_run(da_settings_O1)
+ exp_O1.label = 'O1'
+
+ with open('./config_files/da_exp_O2.json', 'r') as fp:
+ da_settings_O2 = json.load(fp)
+ exp_O2 = load_or_run(da_settings_O2)
+
+ # Plots
+ for exp,settings,label in zip([exp_O1,exp_O2],\
+ [da_settings_O1,da_settings_O2],\
+ [r"O$_1$",r"O$_2$"]):
+ print(exp.randomize_obs)
+
+ plot_overview(exp,label,'exp_%s.png'%exp.label)
+
+ # Difference between two experiments
+ nobs = exp_O1.nobs
+ maxtime = exp_O1.nr.maxtime
+ ntimes = exp_O1.da.ncycles
+ me1 = np.zeros(exp_O1.da.model_equivalents.shape)+np.nan
+ me2 = np.zeros(exp_O2.da.model_equivalents.shape)+np.nan
+
+ # To get the plot right, we need to re-sort the
+ # model equivalents and observation coordinates
+ zobs1 = exp_O1.da.obs_coordinates
+ zobs2 = exp_O2.da.obs_coordinates
+ for k in range(ntimes):
+ if is_sorted(zobs1[k,:]):
+ indices = np.arange(nobs)
+ else:
+ indices = np.argsort(zobs1[k,:])
+ zobs1[k,:] = zobs1[k,indices]
+ me1[k,:,:] = exp_O1.da.model_equivalents[k,indices,:]
+ for k in range(ntimes):
+ if is_sorted(zobs2[k,:]):
+ indices = np.arange(nobs)
+ else:
+ indices = np.argsort(zobs2[k,:])
+ zobs2[k,:] = zobs2[k,indices]
+ me2[k,:,:] = exp_O2.da.model_equivalents[k,indices,:]
+
+ times = exp_O1.da.time
+ times = times[:,None]+np.zeros((ntimes,nobs))
+ delta_theta = me1.mean(axis=2) - me2.mean(axis=2)
+
+ fig, ax = p.subplots(1,1,constrained_layout=True)
+ fig.set_size_inches(4,3)
+ c0 = ax.pcolormesh(times/3600.,zobs1,delta_theta,cmap='RdBu_r')
+ ax.set_xlim([0,maxtime/3600])
+ ax.set_xticks(np.arange((maxtime+1)/3600))
+ ax.set_xlabel('Time (h)')
+ ax.set_ylabel('Height (m)')
+ ax.set_title(r'$\overline{y}_b$(O2)$ - \overline{y}_b$(O1) (K)')
p.colorbar(c0,orientation='horizontal')
- p.colorbar(c1,orientation='horizontal')
- p.colorbar(c2,orientation='horizontal')
- #
- fig.savefig('expRC_PE.png',format='png',dpi=300)
+ fig.savefig('impact_of_random_obs.png',format='png',dpi=300)
p.close(fig)
-
- #########################################################################################à
- # Make plots: imperfect-model experiment (R)
- fig, [ax1, ax2] = p.subplots(1,2,constrained_layout=True)
+
+ if opt09:
+
+ plot_B_matrix(exp_A1,timeid=71)
+
+ if opt10:
+
+ # Plots: Scatter plots of Desroziers statistics
+ #experiments = [exp_A1, exp_B1, exp_B2, exp_B3, exp_B4, exp_C1, exp_C2, exp_D1, exp_D2 ]
+ #labels = [r'A$_1$', r'B$_1$', r'B$_2$', r'B$_3$', r'B$_4$', r'C$_1$', r'C$_2$', r'D$_1$', r'D$_2$' ]
+ experiments1 = [exp_A1, exp_B1, exp_B4, exp_C1, exp_C2 ]
+ labels1 = [r'A$_1$', r'B$_1$', r'B$_4$', r'C$_1$', r'C$_2$' ]
+ experiments2 = [exp_B2, exp_B3, exp_D1, exp_D2 ]
+ labels2 = [r'B$_2$', r'B$_3$', r'D$_1$', r'D$_2$' ]
+
+ fig = p.figure(0)
fig.set_size_inches(6,3)
+ ax1 = fig.add_subplot(1,2,1)
+ ax2 = fig.add_subplot(1,2,2)
- which_cycle = 0
- tt1 = (which_cycle+1)*da_settings_R['assimilation_interval']
- ax1 = plot_CBL_assimilation(exp_R,None,which_cycle=which_cycle,ax=ax1)
- ax1.set_ylim([0,2000])
- ax1.set_xlim([288,298])
- ax1.set_ylabel(r'height (m)')
- ax1.set_title(r'a) Analysis at $t=$%u s'%tt1)
- ax1.legend(frameon=False)
-
- which_cycle = 31
- tt2 = (which_cycle+1)*da_settings_R['assimilation_interval']
- ax2 = plot_CBL_assimilation(exp_R,None,which_cycle=which_cycle,ax=ax2)
- ax2.set_ylim([0,2000])
- ax2.set_xlim([288,298])
- ax2.set_ylabel(r'height (m)')
- ax2.set_title(r'b) Analysis at $t=$%u s'%tt2)
+ ax1 = plot_desroziers([i.dg for i in experiments1],labels=labels1,ax=ax1)
+ ax2 = plot_desroziers([i.dg for i in experiments2],labels=labels2,ax=ax2)
- fig.savefig('expR_profiles.png',format='png',dpi=300)
+ fig.tight_layout()
+ fig.savefig('Diag_desroziers.png',format='png',dpi=150)
p.close(fig)
- fig, [[ax0, ax1], [ax2, ax3]] = p.subplots(2,2,constrained_layout=True)
- fig.set_size_inches(6,6)
- #
- [ax0,ax1,ax2],c0,c1,c2 = plot_CBL_identifiability(exp_R,da_settings_R['obs_error_sdev_assimilate'][0],ax=[ax0,ax1,ax2])
- ax0.set_title(r'a) Exp. R, $\rho(p\prime\prime,y_b}$)')
- ax0.set_xlabel('Time (h)')
- ax0.set_ylabel('Height (m)')
- ax1.set_title(r'b) Exp. R, $\delta\overline{y}\cdot(\sigma_{p\prime\prime}/\sigma_{y^b})$')
- ax1.set_xlabel('Time (h)')
- ax1.set_ylabel('Height (m)')
- ax2.set_title(r'c) Exp. R, $\delta\overline{p}\prime\prime$')
- ax2.set_xlabel('Time (h)')
- ax2.set_ylabel('Height (m)')
- ax3 = plot_CBL_PE(exp_R,None,ax=ax3)
- ax3.set_title(r'd) Exp. R, evolution of $p$')
- ax3.set_xlabel('Time (h)')
- ax3.set_yticks([0,1,2,3,4,5])
- maxtime = exp_R.trun
- ax0.set_xlim([0,maxtime/3600])
- ax1.set_xlim([0,maxtime/3600])
- ax2.set_xlim([0,maxtime/3600])
- ax3.set_xlim([0,maxtime/3600])
- ax0.set_xticks(np.arange((maxtime+1)/3600))
- ax1.set_xticks(np.arange((maxtime+1)/3600))
- ax2.set_xticks(np.arange((maxtime+1)/3600))
- ax3.set_xticks(np.arange((maxtime+1)/3600))
- p.colorbar(c0,orientation='horizontal')
- p.colorbar(c1,orientation='horizontal')
- p.colorbar(c2,orientation='horizontal')
- #
- fig.savefig('expR_PE.png',format='png',dpi=300)
- p.close(fig)
+ if opt11:
- if opt04:
+ # Plots
+ plot_overview(exp_D1,r"D$_1$",'exp_D1.png')
+ plot_overview(exp_D2,r"D$_2$",'exp_D2.png')
+ plot_overview(exp_B3,r"B$_3$",'exp_B3.png')
- # Use only 20 ensemble members for these experiments
- nens = 20
-
- # Base settings
- cbl_settings_base = dict(default_cbl_settings)
- cbl_settings_base['theta_0'] = 290
- cbl_settings_base['gamma'] = 3e-3
- cbl_settings_base['Hmax'] = 0.12
- da_settings_base = dict(default_da_settings)
- da_settings_base['cbl_settings'] = dict(cbl_settings_base)
- da_settings_base['nens'] = nens
- da_settings_base['type'] = 'real'
- da_settings_base['observation_files'] = './LES_data/Theta/*csv'
- da_settings_base['obs_coordinates'] = None
- da_settings_base['obs_error_sdev_generate'] = None
-
- # Base experiment (in opt03 has tspinup = 1 hour). Try other lengths.
- da_settings_spinup_1 = dict(da_settings_base)
- da_settings_spinup_1['tspinup'] = 7200
- da_settings_spinup_1['trun'] = 25200-7200
- da_settings_spinup_1['path'] = "./runs/spinup_seed=%s_enssize=%s_EAKF_6hrs/"%(rnseed,nens)
-
- da_settings_spinup_2 = dict(da_settings_base)
- da_settings_spinup_2['tspinup'] = 10800
- da_settings_spinup_2['trun'] = 25200-10800
- da_settings_spinup_2['path'] = "./runs/spinup_seed=%s_enssize=%s_EAKF_6hrs/"%(rnseed,nens)
-
- da_settings_spinup_3 = dict(da_settings_base)
- da_settings_spinup_3['tspinup'] = 14400
- da_settings_spinup_3['trun'] = 25200-14400
- da_settings_spinup_3['path'] = "./runs/spinup_seed=%s_enssize=%s_EAKF_6hrs/"%(rnseed,nens)
-
- # Run and save to disk
- try:
- exp_spinup_1 = pickle.load(open(da_settings_spinup_1['path']+'exp_spinup_1.pickle', "rb"))
- except:
- exp_spinup_1 = experiment(da_settings_spinup_1)
- setattr(exp_spinup_1,'label','spinup_1')
- pickle.dump(exp_spinup_1, open(da_settings_spinup_1['path']+'exp_spinup_1.pickle', 'wb'))
- try:
- exp_spinup_2 = pickle.load(open(da_settings_spinup_2['path']+'exp_spinup_2.pickle', "rb"))
- except:
- exp_spinup_2 = experiment(da_settings_spinup_2)
- setattr(exp_spinup_2,'label','spinup_2')
- pickle.dump(exp_spinup_2, open(da_settings_spinup_2['path']+'exp_spinup_2.pickle', 'wb'))
- try:
- exp_spinup_3 = pickle.load(open(da_settings_spinup_3['path']+'exp_spinup_3.pickle', "rb"))
- except:
- exp_spinup_3 = experiment(da_settings_spinup_3)
- setattr(exp_spinup_3,'label','spinup_3')
- pickle.dump(exp_spinup_3, open(da_settings_spinup_3['path']+'exp_spinup_3.pickle', 'wb'))
-
- for exp,settings in zip([exp_spinup_1,exp_spinup_2,exp_spinup_3],[da_settings_spinup_1,da_settings_spinup_2,da_settings_spinup_3]):
- label = exp.label
- maxtime = exp.trun
-
- fig, [[ax0, ax1], [ax2, ax3]] = p.subplots(2,2,constrained_layout=True)
- fig.set_size_inches(6,6)
- #
- [ax0,ax1,ax2],c0,c1,c2 = plot_CBL_identifiability(exp,settings['obs_error_sdev_assimilate'][0],ax=[ax0,ax1,ax2])
- ax0.set_title(r'a) Exp. %s, $\rho(p\prime\prime,y_b}$)'%label)
- ax0.set_xlabel('Time (h)')
- ax0.set_ylabel('Height (m)')
- ax1.set_title(r'b) Exp. %s, $\delta\overline{y}\cdot(\sigma_{p\prime\prime}/\sigma_{y^b})$'%label)
- ax1.set_xlabel('Time (h)')
- ax1.set_ylabel('Height (m)')
- ax2.set_title(r'c) Exp. %s, $\delta\overline{p}\prime\prime$'%label)
- ax2.set_xlabel('Time (h)')
- ax2.set_ylabel('Height (m)')
- ax3 = plot_CBL_PE(exp,None,ax=ax3)
- ax3.set_title(r'd) Exp. %s, evolution of $p$'%label)
- ax3.set_xlabel('Time (h)')
- ax3.set_yticks([0,1,2,3,4,5])
- ax0.set_xlim([0,maxtime/3600])
- ax1.set_xlim([0,maxtime/3600])
- ax2.set_xlim([0,maxtime/3600])
- ax3.set_xlim([0,maxtime/3600])
- ax0.set_xticks(np.arange((maxtime+1)/3600))
- ax1.set_xticks(np.arange((maxtime+1)/3600))
- ax2.set_xticks(np.arange((maxtime+1)/3600))
- ax3.set_xticks(np.arange((maxtime+1)/3600))
- p.colorbar(c0,orientation='horizontal')
- p.colorbar(c1,orientation='horizontal')
- p.colorbar(c2,orientation='horizontal')
- #
- fig.savefig('exp_%s_PE.png'%label,format='png',dpi=300)
- p.close(fig)
+ fig, [[ax1, ax2], [ax3, ax4]] = p.subplots(2,2,constrained_layout=True)
+ fig.set_size_inches(6,4)
+ plot_diagnostics(exp_D1,show='errors',axes=[ax1,ax2,ax3,ax4],label='with PE',linecolor='#1f77b4')
+ plot_diagnostics(exp_D2,show='errors',axes=[ax1,ax2,ax3,ax4],label='without PE',linecolor='#9EC9E9')
+ plot_diagnostics(exp_B3,show='errors',axes=[ax1,ax2,ax3,ax4],label='without PE',linecolor='orange')
+ ax1.set_title('a) Analysis error')
+ ax1.set_xlabel(r'RMSE$^a_\theta$')
+ ax2.set_title('b) First-guess error')
+ ax2.set_xlabel(r'RMSE$^b_\theta$')
+ ax3.set_title('c) Error reduction')
+ ax3.set_xlabel(r'RMSE$^b_\theta-$RMSE$^a_\theta$')
+ ax4.set_title('d) Spread-error consistency')
+ ax4.set_xlabel(r'$\sigma^b_\theta$/RMSE$^b_\theta$')
+ ax1.set_ylabel('height (m)')
+ ax3.set_ylabel('height (m)')
+ #ax1.set_xlim([0,0.2])
+ #ax2.set_xlim([0,0.2])
+ #ax3.set_xlim([0,0.02])
+ #ax4.set_xlim([0,10])
+ fig.savefig('profiles_errors.png',format='png',dpi=300)
+ p.close(fig)
+
+ fig, [[ax1, ax2], [ax3, ax4]] = p.subplots(2,2,constrained_layout=True)
+ fig.set_size_inches(6,4)
+ plot_diagnostics(exp_D1,show='deviations',axes=[ax1,ax2,ax3,ax4],label='with PE',linecolor=u'#1f77b4')
+ plot_diagnostics(exp_D2,show='deviations',axes=[ax1,ax2,ax3,ax4],label='without PE',linecolor=u'#9EC9E9')
+ plot_diagnostics(exp_B3,show='deviations',axes=[ax1,ax2,ax3,ax4],label='without PE',linecolor='orange')
+ ax1.set_title('a) Analysis error')
+ ax1.set_xlabel(r'RMSD$^a_\theta$')
+ ax2.set_title('b) First-guess error')
+ ax2.set_xlabel(r'RMSD$^b_\theta$')
+ ax3.set_title('c) Error reduction')
+ ax3.set_xlabel(r'RMSD$^b_\theta-$RMSD$^a_\theta$')
+ ax4.set_title('d) Spread-error consistency')
+ ax4.set_xlabel(r'$\sigma^b_\theta$/RMSD$^b_\theta$')
+ ax1.set_ylabel('height (m)')
+ ax3.set_ylabel('height (m)')
+ #ax1.set_xlim([0,0.2])
+ #ax2.set_xlim([0,0.2])
+ #ax3.set_xlim([0,0.02])
+ #ax4.set_xlim([0,10])
+ fig.savefig('profiles_deviations.png',format='png',dpi=300)
+ p.close(fig)
- check_reproducibility = False
if check_reproducibility:
- path = default_da_settings['path']
- for i in ["exp_A.pickle",
- "exp_A_noPE.pickle",
+
+ for i in ["exp_A1.pickle",
+ "exp_A2.pickle",
"exp_B1.pickle",
"exp_B2.pickle",
"exp_B3.pickle",
"exp_B4.pickle",
"exp_C1.pickle",
"exp_C2.pickle",
- "exp_D.pickle"]:
- experiment = pickle.load(open(path+i, "rb"))
+ "exp_D1.pickle",
+ "exp_D2.pickle"]:
+ experiment = pickle.load(open(workdir+i, "rb"))
print(i)
print("OBS :",experiment.da.observations[:,0])
print("PARAMETERS:",experiment.da.initial_perturbed_parameters[0,:5])
- print("STATE :",experiment.da.initial_state[0,0,:5])
-
+ print("STATE :",experiment.da.initial_state[0,0,:5])
\ No newline at end of file