WOFOST模型的卡尔曼滤波算法（EnKF）（MulitiState）

import sys, os
import copy
import datetime as dt
import pcse
from pcse.base import ParameterProvider
from pcse.fileinput import CABOFileReader, YAMLAgroManagementReader, ExcelWeatherDataProvider
from pcse.models import Wofost72_WLP_FD, Wofost71_WLP_FD
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pcse.util import WOFOST71SiteDataProvider

plt.style.use("ggplot")
import matplotlib.pyplot as plt
print("This notebook was built with:")
print("python version: %s " % sys.version)
print("PCSE version: %s" %  pcse.__version__)


states_for_DA = ["LAI", "WST", "WSO"]
observed_states = ["LAI"]
dates_of_observation = [dt.date(2006,4,9), dt.date(2006,5,2), dt.date(2006,6,17),
                        dt.date(2006,8,14), dt.date(2006,9,12)]
observed_lai = np.array([2.2, 3.5, 6.2, 3.3, 2.1])
std_lai = observed_lai * 0.1 # Std. devation is estimated as 10% of observed value
observations_for_DA = []
# Pack them into a convenient format
for d, lai, errlai, in zip(dates_of_observation, observed_lai, std_lai):
    observations_for_DA.append((d, {"LAI":(lai, errlai)}))

ensemble_size = 50
np.random.seed(10000)
data_dir = r'D:\\wofost\\wofost_learn'

cropfile = os.path.join(data_dir, 'sug0601.crop')
cropdata = CABOFileReader(cropfile)


soilfile = os.path.join(data_dir, 'ec3.soil')
soildata = CABOFileReader(soilfile)
sitedata = WOFOST71SiteDataProvider(WAV=100, CO2=360)


parameters = ParameterProvider(cropdata=cropdata, soildata=soildata,sitedata=sitedata)

agromanagement_file = os.path.join(data_dir, 'sugarbeet_calendar.agro')
agromanagement = YAMLAgroManagementReader(agromanagement_file)


wdp = ExcelWeatherDataProvider(os.path.join(data_dir,'nl1.xlsx'))


wofost = Wofost71_WLP_FD(parameters, wdp, agromanagement)
wofost1 = Wofost71_WLP_FD(parameters, wdp, agromanagement)
wofost.run_till_terminate()
df = pd.DataFrame(wofost.get_output()).set_index("day")
df.to_excel("wofost_results_enkf.xlsx")
output = wofost.get_output()

# A container for the parameters that we will override
override_parameters = {}
#Initial conditions
override_parameters["TDWI"] = np.random.normal(150., 50., (ensemble_size))
override_parameters["WAV"] = np.random.normal(4.5, 1.5, (ensemble_size))
# parameters
override_parameters["SPAN"] = np.random.normal(31, 3 ,(ensemble_size))
override_parameters["SMFCF"] = np.random.normal(0.31, 0.03 ,(ensemble_size))

ensemble = []
for i in range(ensemble_size):
    p = copy.deepcopy(parameters)
    for par, distr in override_parameters.items():
        p.set_override(par, distr[i])
    member = Wofost71_WLP_FD(p, wdp, agromanagement)
    ensemble.append(member)
count = 0

while count<5:
    count+=1
    day, obs = observations_for_DA.pop(0)
    for member in ensemble:
        member.run_till(day)
    print("Ensemble now at day %s" % member.day)
    print("%s observations left!" % len(observations_for_DA))

    collected_states = []
    for member in ensemble:
        t = {}
        for state in states_for_DA:
            t[state] = member.get_variable(state)
        collected_states.append(t)
    df_A = pd.DataFrame(collected_states)
    A = np.matrix(df_A).T
    P_e = np.matrix(df_A.cov())

    perturbed_obs = []
    observations = ["LAI"]
    for state in observed_states:
        (value, std) = obs[state]
        d = np.random.normal(value, std, (ensemble_size))
        perturbed_obs.append(d)
    df_perturbed_obs = pd.DataFrame(perturbed_obs).T
    df_perturbed_obs.columns = observations
    D = np.matrix(df_perturbed_obs).T  # Perturbed observations
    R_e = np.matrix(df_perturbed_obs.cov())  # Covariance

    H = np.matrix([1.,0.,0.])
    K1 = P_e * (H.T)
    K2 = (H * P_e) * H.T
    K = K1 * ((K2 + R_e).I)

    # Here we compute the analysed states
    Aa = A + K * (D - (H * A))
    df_Aa = pd.DataFrame(Aa.T, columns=states_for_DA)

    for member, new_states in zip(ensemble, df_Aa.itertuples()):
        member.set_variable("LAI", new_states.LAI)
        member.set_variable("WST", new_states.WST)
        member.set_variable("WSO", new_states.WSO)


for member in ensemble:
    member.run_till_terminate()

results = []
for member in ensemble:
    member_df = pd.DataFrame(member.get_output()).set_index("day")
    results.append(member_df)

df_A.corr()

fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(16,30))
for member_df in results:
    member_df["LAI"].plot(style="k:", ax=axes[0])
    member_df["TWST"].plot(style="k:", ax=axes[1])
    member_df["TWSO"].plot(style="k:", ax=axes[2])
    member_df["TAGP"].plot(style="k:", ax=axes[3])
axes[0].errorbar(dates_of_observation, observed_lai, yerr=std_lai, fmt="o")
axes[0].set_title("Leaf area index")
axes[1].set_title("Weight stems (WST)")
axes[2].set_title("Weight storage organs (WST)")
axes[3].set_title("Total biomass")
fig.autofmt_xdate()
plt.show()