ModVege with EURO-CORDEX data#

  • Jouven, M., Carrère, P., and Baumont, R. (2006a). ‘Model predicting dynamics of biomass, structure and digestibility of herbage in managed permanent pastures. 1. Model description’, Grass and Forage Science, vol. 61, no. 2, pp. 112-124. DOI: 10.1111/j.1365-2494.2006.00515.x.

  • Jouven, M., Carrère, P., and Baumont, R. (2006b). ‘Model predicting dynamics of biomass, structure and digestibility of herbage in managed permanent pastures. 2. Model evaluation’, Grass and Forage Science, vol. 61, no. 2, pp. 125-133. DOI: 10.1111/j.1365-2494.2006.00517.x.

  • Agri4cast (2023). ‘ModVege’. Available at: https://code.europa.eu/agri4cast/modvege (Accessed: 28 August 2023).

import glob
import os
from datetime import datetime, timezone
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr
import climag.plot_configs as cplt
from climag import climag_plot

DATA_DIR = os.path.join("data", "ModVege")

# Ireland boundary
GPKG_BOUNDARY = os.path.join("data", "boundaries", "boundaries_all.gpkg")
ie = gpd.read_file(GPKG_BOUNDARY, layer="NUTS_RG_01M_2021_2157_IE")
ie_bbox = gpd.read_file(
    GPKG_BOUNDARY, layer="NUTS_RG_01M_2021_2157_IE_BBOX_DIFF"
)

# met station coords
LON, LAT = -8.26389, 52.16389  # Moorepark, Fermoy

Historical data#

data = xr.open_mfdataset(
    glob.glob(
        os.path.join(DATA_DIR, "EURO-CORDEX", "historical", "EC-EARTH", "*.nc")
    ),
    chunks="auto",
    decode_coords="all",
)

data

<xarray.Dataset>
Dimensions:       (rlat: 33, rlon: 37, time: 11323, bnds: 2)
Coordinates:
    lat           (rlat, rlon) float64 dask.array<chunksize=(33, 37), meta=np.ndarray>
    lon           (rlat, rlon) float64 dask.array<chunksize=(33, 37), meta=np.ndarray>
  * rlat          (rlat) float64 3.685 3.795 3.905 4.015 ... 6.985 7.095 7.205
  * rlon          (rlon) float64 -17.27 -17.16 -17.05 ... -13.53 -13.41 -13.3
  * time          (time) datetime64[ns] 1975-01-01T12:00:00 ... 2005-12-31T12...
    height        float64 2.0
    rotated_pole  |S1 b''
    time_bnds     (time, bnds) datetime64[ns] dask.array<chunksize=(365, 2), meta=np.ndarray>
    spatial_ref   int64 0
Dimensions without coordinates: bnds
Data variables: (12/24)
    bm_gv         (time, rlat, rlon) float32 dask.array<chunksize=(365, 33, 37), meta=np.ndarray>
    bm_gr         (time, rlat, rlon) float32 dask.array<chunksize=(365, 33, 37), meta=np.ndarray>
    bm_dv         (time, rlat, rlon) float32 dask.array<chunksize=(365, 33, 37), meta=np.ndarray>
    bm_dr         (time, rlat, rlon) float32 dask.array<chunksize=(365, 33, 37), meta=np.ndarray>
    age_gv        (time, rlat, rlon) float32 dask.array<chunksize=(365, 33, 37), meta=np.ndarray>
    age_gr        (time, rlat, rlon) float32 dask.array<chunksize=(365, 33, 37), meta=np.ndarray>
    ...            ...
    sen_gv        (time, rlat, rlon) float32 dask.array<chunksize=(365, 33, 37), meta=np.ndarray>
    sen_gr        (time, rlat, rlon) float32 dask.array<chunksize=(365, 33, 37), meta=np.ndarray>
    abs_dv        (time, rlat, rlon) float32 dask.array<chunksize=(365, 33, 37), meta=np.ndarray>
    abs_dr        (time, rlat, rlon) float32 dask.array<chunksize=(365, 33, 37), meta=np.ndarray>
    omd_gv        (time, rlat, rlon) float32 dask.array<chunksize=(365, 33, 37), meta=np.ndarray>
    omd_gr        (time, rlat, rlon) float32 dask.array<chunksize=(365, 33, 37), meta=np.ndarray>
Attributes:
    creation_date:  2023-03-11 00:16:37.316077+00:00
    contact:        nstreethran@ucc.ie
    frequency:      day
    references:     https://github.com/ClimAg
    input_dataset:  IE_EURO-CORDEX_RCA4_EC-EARTH_historical
# remove the spin-up year
data = data.sel(time=slice("1976", "2005"))

Seasonal averages#

for var in ["gro", "pgro", "bm", "c_bm"]:
    climag_plot.plot_averages(
        data=data,
        var=var,
        averages="season",
        boundary_data=ie_bbox,
        cbar_levels=12,
    )
../_images/266d7a9d26a9b8823828f80253427d828f41e626f81682129037dd334c642f73.png ../_images/0f120ba28db2db4766aa82ff0333efa36714a97dfff7a6dae5df3ee7cba0a7a0.png ../_images/71850d4e080545d65b99699181a096d90290042cfdae2cf7d8ee3e3a488f2e2b.png ../_images/e61c31cfb2c8fac00db73ede9f9f8e6baaaa6863477a1e6d300d7a5c4a4e8e13.png

Point subset#

cds = cplt.rotated_pole_point(data=data, lon=LON, lat=LAT)

data_ie = data.sel({"rlon": cds[0], "rlat": cds[1]}, method="nearest")

data_ie

<xarray.Dataset>
Dimensions:       (time: 11323, bnds: 2)
Coordinates:
    lat           float64 dask.array<chunksize=(), meta=np.ndarray>
    lon           float64 dask.array<chunksize=(), meta=np.ndarray>
    rlat          float64 4.235
    rlon          float64 -15.83
  * time          (time) datetime64[ns] 1975-01-01T12:00:00 ... 2005-12-31T12...
    height        float64 2.0
    rotated_pole  |S1 b''
    time_bnds     (time, bnds) datetime64[ns] dask.array<chunksize=(365, 2), meta=np.ndarray>
    spatial_ref   int64 0
Dimensions without coordinates: bnds
Data variables: (12/24)
    bm_gv         (time) float32 dask.array<chunksize=(365,), meta=np.ndarray>
    bm_gr         (time) float32 dask.array<chunksize=(365,), meta=np.ndarray>
    bm_dv         (time) float32 dask.array<chunksize=(365,), meta=np.ndarray>
    bm_dr         (time) float32 dask.array<chunksize=(365,), meta=np.ndarray>
    age_gv        (time) float32 dask.array<chunksize=(365,), meta=np.ndarray>
    age_gr        (time) float32 dask.array<chunksize=(365,), meta=np.ndarray>
    ...            ...
    sen_gv        (time) float32 dask.array<chunksize=(365,), meta=np.ndarray>
    sen_gr        (time) float32 dask.array<chunksize=(365,), meta=np.ndarray>
    abs_dv        (time) float32 dask.array<chunksize=(365,), meta=np.ndarray>
    abs_dr        (time) float32 dask.array<chunksize=(365,), meta=np.ndarray>
    omd_gv        (time) float32 dask.array<chunksize=(365,), meta=np.ndarray>
    omd_gr        (time) float32 dask.array<chunksize=(365,), meta=np.ndarray>
Attributes:
    creation_date:  2023-03-11 00:16:37.316077+00:00
    contact:        nstreethran@ucc.ie
    frequency:      day
    references:     https://github.com/ClimAg
    input_dataset:  IE_EURO-CORDEX_RCA4_EC-EARTH_historical
data_ie = data.sel({"rlon": cds[0], "rlat": cds[1]}, method="nearest")

data_ie_df = pd.DataFrame({"time": data_ie["time"]})
for var in data_ie.data_vars:
    data_ie_df[var] = data_ie[var]

data_ie_df.set_index("time", inplace=True)
data_ie_df = data_ie_df[["pgro", "gro", "h_bm", "i_bm", "bm"]]

# configure plot title
plot_title = []
for var in list(data_ie_df):
    plot_title.append(
        f"{data_ie[var].attrs['long_name']} [{data_ie[var].attrs['units']}]"
    )

data_ie_df.plot(
    subplots=True,
    layout=(5, 1),
    figsize=(12, 14),
    legend=False,
    xlabel="",
    title=plot_title,
    linewidth=1,
)

plt.tight_layout()
plt.show()
../_images/8c78331bb40b368de057d2c03fd002a088e7e5a06b98a67f73169fe14aab9aa9.png 
data_ie = data.sel({"rlon": cds[0], "rlat": cds[1]}, method="nearest").sel(
    time=slice("1990", "1999")
)

data_ie_df = pd.DataFrame({"time": data_ie["time"]})

# configure plot title
plot_title = []

for var in data_ie.data_vars:
    data_ie_df[var] = data_ie[var]
    plot_title.append(
        f"{data_ie[var].attrs['long_name']} [{data_ie[var].attrs['units']}]"
    )

data_ie_df.set_index("time", inplace=True)

data_ie_df.plot(
    subplots=True,
    layout=(9, 3),
    figsize=(18, 18),
    legend=False,
    xlabel="",
    title=plot_title,
    linewidth=1,
)

plt.tight_layout()
plt.show()
../_images/f5b272e5ba17446c6ab09dd24c5a4eeb09295468b091e35048aaccae273aa188.png 
data_ie_y = data_ie.sel(time=slice("1997", "1999"))

data_ie_df = pd.DataFrame({"time": data_ie_y["time"]})

for var in data_ie_y.data_vars:
    data_ie_df[var] = data_ie_y[var]

data_ie_df.set_index("time", inplace=True)
data_ie_df = data_ie_df[["gro"]]

# resample to weekly
data_ie_df = data_ie_df.resample("W-MON").mean()

data_ie_df.reset_index(inplace=True)
mn = data_ie_df.rolling(3, center=True, on="time")["gro"].mean()
data_ie_df["outlier"] = data_ie_df["gro"].sub(mn).abs().gt(10)
data_ie_df["moving_avg"] = mn
data_ie_df.set_index("time", inplace=True)

axs = data_ie_df.plot(y="gro", figsize=(12, 5), xlabel="", label="growth")
data_ie_df.plot(y="moving_avg", ax=axs, color="orange", zorder=1)
if True in list(data_ie_df["outlier"].unique()):
    data_ie_df[data_ie_df["outlier"] == True].plot(
        y="gro",
        ax=axs,
        marker="*",
        linewidth=0,
        color="crimson",
        label="outlier",
    )
plt.xlabel("")
plt.ylabel("Grass growth [kg DM ha⁻¹ day⁻¹]")
plt.tight_layout()
plt.show()
../_images/5db4d94ba910795b04bf82cd255d7d127dcabd7b91ef7e00208b136ee1339804.png