01-generate_pheno.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import datetime\n",
    "import statsmodels.api as sm\n",
    "from scripts.utils import save_fig\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (12, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LAST_YEAR = 2022"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pheno_1 = pd.read_excel(\"data/pheno/pheno_TL.xlsx\", index_col=0).rename(\n",
    "    columns={\n",
    "        \"Pousse de la vigne\": \"Budburst - Source 1\",\n",
    "        \"Première fleur\": \"Beginning flowering - Source 1\",\n",
    "        \"Pleine fleur\": \"End flowering - Source 1\",\n",
    "        \"Premiers grains changés\": \"Beginning véraison - Source 1\",\n",
    "        \"Début vendanges (en Médoc)\": \"Beginning of the first harvests - Source 1\",\n",
    "        \"Pleines vendanges\": \"Beginning of all harvests - Source 1\",\n",
    "        \"Fin vendanges rouges\": \"End harvest - Source 1\",\n",
    "    }\n",
    ")\n",
    "pheno_1 = pheno_1.iloc[:, :7]\n",
    "cols = list(pheno_1.columns)\n",
    "years = list(pheno_1.index)\n",
    "pheno_1_plot = pheno_1.copy()\n",
    "\n",
    "for col in cols:\n",
    "    pheno_1[col] = pd.to_datetime(pheno_1[col], errors=\"coerce\")\n",
    "    for year in years:\n",
    "        try:\n",
    "            pheno_1.loc[year, col] = (\n",
    "                pheno_1.loc[year, col].replace(year=year).normalize()\n",
    "            )\n",
    "        except:\n",
    "            continue\n",
    "    pheno_1_plot[col] = pheno_1[col].dt.dayofyear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pheno_2 = pd.read_excel(\"data/pheno/pheno_latour.xlsx\", index_col=0).rename(\n",
    "    columns={\n",
    "        \"debourrement\": \"Budburst - Château Latour\",\n",
    "        \"half-veraison\": \"Half-véraison - Château Latour\",\n",
    "        \"half-vendanges\": \"Half-harvest - Château Latour\",\n",
    "        \"begin vendanges\": \"Beginning harvest - Château Latour\",\n",
    "        \"end vendanges\": \"End harvest - Château Latour\",\n",
    "    }\n",
    ")\n",
    "pheno_2[\"Half-flowering - Château Latour\"] = (\n",
    "    pheno_2[\"half-floraison M\"]\n",
    "    + (pheno_2[\"half-floraison CS\"] - pheno_2[\"half-floraison M\"]) / 2\n",
    ")\n",
    "pheno_2 = pheno_2.drop([\"half-floraison M\", \"half-floraison CS\"], axis=1)\n",
    "# pheno_2['floraison'] = pheno_2['mi-floraison M']\n",
    "pheno_2_plot = pheno_2.copy()\n",
    "for col in pheno_2.columns:\n",
    "    pheno_2[col] = pd.to_datetime(pheno_2[col], errors=\"coerce\")\n",
    "    for year in years[11:]:\n",
    "        try:\n",
    "            pheno_2.loc[year, col] = (\n",
    "                pheno_2.loc[year, col].replace(year=year).normalize()\n",
    "            )\n",
    "        except:\n",
    "            continue\n",
    "    pheno_2_plot[col] = pheno_2[col].dt.dayofyear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pheno_inrae = pd.read_csv(\n",
    "    \"data/pheno/pheno_INRAE/phenoclim_agroclim_inrae.csv\", sep=\";\", encoding=\"latin1\"\n",
    ").drop(\"source_donnees\", axis=1)\n",
    "pheno_inrae.date = pd.to_datetime(pheno_inrae.date)\n",
    "pheno_inrae[\"dayofyear\"] = pheno_inrae[\"date\"].dt.dayofyear\n",
    "pheno_inrae[\"year\"] = pheno_inrae[\"date\"].dt.year\n",
    "\n",
    "pheno_inrae = pheno_inrae[\n",
    "    [\n",
    "        \"date\",\n",
    "        \"year\",\n",
    "        \"dayofyear\",\n",
    "        \"nom_site\",\n",
    "        \"echelle_phenologique\",\n",
    "        \"description_evenement_phenologique_majeur\",\n",
    "        \"variete\",\n",
    "        \"code_stade_phenologique\",\n",
    "        \"description_stade_phenologique\",\n",
    "        \"degats_gel_de_printemps\",\n",
    "    ]\n",
    "]\n",
    "pheno_inrae = pheno_inrae[\n",
    "    (pheno_inrae[\"nom_site\"].isin([\"INRA Bordeaux - Villenave d'Ornon\"]))\n",
    "    & (pheno_inrae[\"variete\"].isin([\"Merlot\", \"Cabernet-Sauvignon\"]))\n",
    "].drop(\"nom_site\", axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pheno_inrae = (\n",
    "    pheno_inrae.groupby([\"year\", \"description_evenement_phenologique_majeur\"])[\n",
    "        [\"dayofyear\"]\n",
    "    ]\n",
    "    .mean()\n",
    "    .reset_index()\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 6))\n",
    "\n",
    "sns.lineplot(data=pheno_1_plot, palette=\"colorblind\", dashes=False)\n",
    "\n",
    "sns.lineplot(\n",
    "    data=pheno_inrae,\n",
    "    x=\"year\",\n",
    "    y=\"dayofyear\",\n",
    "    hue=\"description_evenement_phenologique_majeur\",\n",
    "    palette=[\"green\", \"orange\", \"brown\"],\n",
    ")\n",
    "\n",
    "palette = sns.color_palette(\"viridis\", 6)\n",
    "sns.lineplot(data=pheno_2_plot, palette=palette)\n",
    "plt.legend(bbox_to_anchor=(1, 1), loc=\"upper left\")\n",
    "plt.ylabel(\"Day of year\")\n",
    "plt.xlabel(\"Vintage\")\n",
    "plt.tight_layout()\n",
    "plt.gca().invert_yaxis()\n",
    "handles, labels = plt.gca().get_legend_handles_labels()\n",
    "# order = [0, 7, 1, -1, 2, 3, 8, 4, -3, 5, -4, -2, 6]\n",
    "plt.legend(\n",
    "    # [handles[idx] for idx in order],\n",
    "    # [labels[idx] for idx in order],\n",
    "    loc=\"upper left\",\n",
    "    bbox_to_anchor=(1, 1),\n",
    ")\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"views/pheno_both_sources.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the above chart, several remarks:\n",
    "\n",
    "- Debourrement dates seem incertain. TL debourrement should always be first, but it is not. Therefore the debourrement will be averaged between the two sources.\n",
    "- \"mi-floraison\" in Source 2 seems to correspond to \"pleine fleur\" in Source 1. The latest date will be chosen for the \"end floraison\" parameter\n",
    "- Missing values in begin and half veraison will be obtained by substracting or adding to the present value the mean difference between series.\n",
    "- Latest and earliest records of any vendanges event will be kept for begin and end of vendanges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pheno_1 = pheno_1.loc[1950:LAST_YEAR]\n",
    "pheno_2 = pheno_2.loc[1950:LAST_YEAR]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "year_range = range(1950, LAST_YEAR+1)\n",
    "pheno_df = pd.DataFrame(np.zeros(len(year_range)), index=year_range).iloc[:, 1:]\n",
    "\n",
    "### Budburst\n",
    "# Calculate mean\n",
    "pheno_df[\"Budburst\"] = (\n",
    "    pheno_1[\"Budburst - Source 1\"]\n",
    "    + (pheno_2[\"Budburst - Château Latour\"] - pheno_1[\"Budburst - Source 1\"]) / 2\n",
    ")\n",
    "# Fill missing values with either only source 1 or only source 2\n",
    "pheno_df[\"Budburst\"] = (\n",
    "    pheno_df[\"Budburst\"]\n",
    "    .fillna(pheno_2[\"Budburst - Château Latour\"])\n",
    "    .fillna(pheno_1[\"Budburst - Source 1\"])\n",
    "    .dt.normalize()\n",
    ")\n",
    "\n",
    "\n",
    "### Floraison\n",
    "pheno_df[\"Beginning flowering\"] = pheno_1[\"Beginning flowering - Source 1\"]\n",
    "\n",
    "# For the end: get the latest known date between end floraison and half-floraison\n",
    "pheno_df[\"End flowering\"] = np.nan\n",
    "for year in pheno_df.index[11:]:\n",
    "    try:\n",
    "        date_1 = pheno_1.loc[year, \"End flowering - Source 1\"]\n",
    "        date_2 = pheno_2.loc[year, \"Half-flowering - Château Latour\"]\n",
    "        pheno_df.loc[year, \"End flowering\"] = max(date_1, date_2)\n",
    "    except:\n",
    "        pheno_df.loc[year, \"End flowering\"] = pheno_2.loc[year, \"Half-flowering - Château Latour\"]\n",
    "pheno_df[\"End flowering\"] = (\n",
    "    pheno_df[\"End flowering\"]\n",
    "    .fillna(pheno_1[\"End flowering - Source 1\"])\n",
    "    .fillna(pheno_2[\"Half-flowering - Château Latour\"])\n",
    "    .dt.normalize()\n",
    ")\n",
    "\n",
    "pheno_df[\"Beginning flowering\"] = pheno_df[\"Beginning flowering\"].fillna((\n",
    "    pheno_df[\"End flowering\"] + (pheno_df[\"Beginning flowering\"] - pheno_df[\"End flowering\"]).mean()\n",
    "))\n",
    "\n",
    "### Veraison\n",
    "average_veraison_delay = np.mean(\n",
    "    pheno_2[\"Half-véraison - Château Latour\"] - pheno_1[\"Beginning véraison - Source 1\"]\n",
    ")\n",
    "pheno_df[\"Beginning véraison\"] = (\n",
    "    pheno_1[\"Beginning véraison - Source 1\"]\n",
    "    .fillna(pheno_2[\"Half-véraison - Château Latour\"] - average_veraison_delay)\n",
    "    .dt.normalize()\n",
    ")\n",
    "pheno_df[\"Half véraison\"] = pheno_2[\"Half-véraison - Château Latour\"]\n",
    "pheno_df[\"Half véraison\"] = (\n",
    "    pheno_df[\"Half véraison\"]\n",
    "    .fillna(pheno_1[\"Beginning véraison - Source 1\"] + average_veraison_delay)\n",
    "    .dt.normalize()\n",
    ")\n",
    "\n",
    "\n",
    "### Vendanges\n",
    "# For the begin: get the earliest known date between the two sources\n",
    "pheno_df[\"Beginning harvest\"] = np.nan\n",
    "\n",
    "for year in pheno_df.index[11:]:\n",
    "    try:\n",
    "        date_1 = pheno_1.loc[year, \"Beginning of the first harvests - Source 1\"]\n",
    "        date_2 = pheno_2.loc[year, \"Beginning harvest - Château Latour\"]\n",
    "        pheno_df.loc[year, \"Beginning harvest\"] = min(date_1, date_2)\n",
    "    except:\n",
    "        pheno_df.loc[year, \"Beginning harvest\"] = pheno_2.loc[year, \"Beginning harvest - Château Latour\"]\n",
    "pheno_df[\"Beginning harvest\"] = pheno_df[\"Beginning harvest\"].fillna(\n",
    "    pheno_1[\"Beginning of all harvests - Source 1\"]\n",
    ")\n",
    "\n",
    "# pheno_df['Beginning harvest'] = pheno_df['Beginning harvest'].dt.normalize()\n",
    "# For the end: get the latest known date between the two sources\n",
    "pheno_df[\"End harvest\"] = np.nan\n",
    "for year in pheno_df.index[11:]:\n",
    "    try:\n",
    "        date_1 = pheno_1.loc[year, \"End harvest - Source 1\"]\n",
    "        date_2 = pheno_2.loc[year, \"End harvest - Château Latour\"]\n",
    "        pheno_df.loc[year, \"End harvest\"] = max(date_1, date_2)\n",
    "    except:\n",
    "        pheno_df.loc[year, \"End harvest\"] = pheno_2.loc[year, \"End harvest - Château Latour\"]\n",
    "pheno_df[\"End harvest\"] = (\n",
    "    pheno_df[\"End harvest\"]\n",
    "    .fillna(pheno_1[\"End harvest - Source 1\"])\n",
    "    .fillna(pheno_2[\"End harvest - Château Latour\"])\n",
    "    .dt.normalize()\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pheno_df_plot = pheno_df.copy()\n",
    "for col in pheno_df_plot.columns:\n",
    "    pheno_df_plot[col] = pheno_df_plot[col].dt.dayofyear\n",
    "font = {'weight' : 'regular',\n",
    "        'family': \"serif\",\n",
    "        'size'   : 16}\n",
    "plt.rc('font', **font)\n",
    "\n",
    "plt.figure(figsize=(10,5))\n",
    "plt.grid(True, which='both', axis='both', alpha=0.5)\n",
    "\n",
    "sns.lineplot(\n",
    "    data=pheno_df_plot.drop(\"Beginning véraison\", axis=1).rename(columns={\"Beginning flowering\": \"First flowers\"}),\n",
    "    markers=[\"s\", \"o\", \"s\", \"s\", \"o\", \"s\"],\n",
    "    dashes=False,\n",
    "    palette=\"colorblind\",\n",
    "    markersize=6,\n",
    "    linewidth=1.5,\n",
    "    markeredgecolor=\"k\",\n",
    ")\n",
    "\n",
    "plt.xlim(1947, LAST_YEAR)\n",
    "plt.xlabel(\"Year\")\n",
    "plt.ylabel(\"Day of year\")\n",
    "\n",
    "handles, labels = plt.gca().get_legend_handles_labels()\n",
    "legend = plt.legend(\n",
    "    reversed(handles),\n",
    "    reversed(labels),\n",
    "    bbox_to_anchor=(1,0.99),\n",
    "    loc='upper left',\n",
    "    frameon=False,\n",
    "    labelspacing=2,\n",
    "    handletextpad=0.8,\n",
    "    #ncol=3,\n",
    "    fontsize=15,\n",
    "    handlelength=0,\n",
    ")\n",
    "\n",
    "save_fig('views/pheno_full', \"single\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pheno_df['Budburst - flowering interval'] = (pheno_df['Beginning flowering'] - pheno_df['Budburst']).dt.days\n",
    "pheno_df['Flowering - véraison interval'] = (pheno_df['Half véraison'] - pheno_df['End flowering']).dt.days\n",
    "pheno_df['Véraison - harvest interval'] = (pheno_df['Beginning harvest'] - pheno_df['Half véraison']).dt.days\n",
    "pheno_df['Growing season duration'] = (pheno_df['Beginning harvest'] - pheno_df['Budburst']).dt.days\n",
    "display(pheno_df.tail())\n",
    "\n",
    "pheno_df.to_excel('data/pheno/generated_pheno.xlsx')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Display average dates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dayofyears = pheno_df[\n",
    "    [\n",
    "        \"Budburst\",\n",
    "        \"Beginning flowering\",\n",
    "        \"End flowering\",\n",
    "        \"Beginning véraison\",\n",
    "        \"Half véraison\",\n",
    "        \"Beginning harvest\",\n",
    "        \"End harvest\",\n",
    "    ]\n",
    "].apply(lambda column: column.dt.dayofyear, axis=1)\n",
    "dayofyears.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.to_datetime(dayofyears.loc[1960:LAST_YEAR].mean(), unit='D', origin=datetime.datetime(year=2001, month=1, day=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "for column in dayofyears.columns:\n",
    "    X = sm.add_constant(dayofyears.index.values.reshape(-1,1))\n",
    "    y = dayofyears[column]\n",
    "    print(column.upper(), sm.OLS(y, X).fit().summary(), \"\\n\\n\\n\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# End of notebook"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.15"
  },
  "vscode": {
   "interpreter": {
    "hash": "0ac38357a05bc42f0ed259dcad053212499e7f4648ddbc814451284e59f63394"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}