Ensemble without assimilating
1 In:
# general python
import warnings
warnings.filterwarnings("ignore", category=UserWarning)
import numpy as np
import os
from pathlib import Path
import yaml
import pandas as pd
import matplotlib.pyplot as plt
from datetime import datetime
import scipy
import xarray as xr
from tqdm import tqdm
import glob
from devtools import pprint
from tqdm import tqdm
2 In:
# general eWC
import ewatercycle
import ewatercycle.forcing
import ewatercycle.models
set up paths
3 In:
path = Path.cwd()
forcing_path = path / "Forcing"
observations_path = path / "Observations"
figure_path = path / "Figures"
output_path = path / "Output"
Marrmot
add forcing
4 In:
from ewatercycle.forcing import sources
5 In:
marrmot_forcing = ewatercycle.forcing.sources["MarrmotForcing"](
directory=forcing_path,
start_time="1989-01-01T00:00:00Z",
end_time="1992-12-31T00:00:00Z",
forcing_file="BMI_testcase_m01_BuffaloRiver_TN_USA.mat",
)
import model
6 In:
ewatercycle.models.sources.MarrmotM14
6 Out:
ewatercycle_marrmot.model.MarrmotM14
7 In:
from ewatercycle.models import MarrmotM14
import DA function
pip install ewatercycle_DA
11 In:
from ewatercycle_DA import DA
12 In:
n_particles = 3
13 In:
ensemble = DA.Ensemble(N=n_particles)
ensemble.setup()
14 In:
maximum_soil_moisture_storage_lst = [10.0, 12.0, 14.0]
15 In:
# values wihch you
setup_kwargs_lst = []
for index in range(n_particles):
setup_kwargs_lst.append({'maximum_soil_moisture_storage': maximum_soil_moisture_storage_lst[index],
'end_time':'1989-02-01T00:00:00Z',
})
16 In:
ensemble.ensemble_list[0].model_name
18 In:
# this initializes the models for all ensemble members. - ignore warnings for now
ensemble.initialize(model_name=["MarrmotM14"]*n_particles,
forcing=[marrmot_forcing]*n_particles,
setup_kwargs=setup_kwargs_lst)
Run
19 In:
ref_model = ensemble.ensemble_list[0].model
20 In:
n_timesteps = int((ref_model.end_time - ref_model.start_time) / ref_model.time_step)
time = []
lst_Q = []
for i in tqdm(range(n_timesteps)):
time.append(pd.Timestamp(ref_model.time_as_datetime.date()))
ensemble.update()
lst_Q.append(ensemble.get_value("flux_out_Q").flatten())
ensemble.finalize()
100%|███████████████████████████████████████████████████████████████████████████████████| 31/31 [00:02<00:00, 14.34it/s]
21 In:
Q_m_arr = np.array(lst_Q).T
df_ensemble = pd.DataFrame(data=Q_m_arr[:,:len(time)].T,index=time,columns=[f'particle {n}' for n in range(n_particles)])
22 In:
ax = df_ensemble[["particle 0"]].plot(lw=2.5)
df_ensemble[["particle 1"]].plot(ax=ax,ls="--",lw=2.5)
df_ensemble[["particle 2"]].plot(ax=ax,ls="-.",lw=1.5)
22 Out:
<Axes: >