04-create_features.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from tqdm.auto import tqdm\n",
    "from statsmodels.stats.outliers_influence import variance_inflation_factor\n",
    "from itertools import product\n",
    "\n",
    "from scripts.feature_engineering import add_base_features, create_df\n",
    "from scripts.utils import save_fig\n",
    "%reload_ext autoreload\n",
    "%autoreload 2\n",
    "plt.rcParams['figure.figsize'] = (12, 4)\n",
    "font = {'size'   : 18}\n",
    "plt.rc('font', **font)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LIST = ['ERA5', 'METEO-FRANCE']\n",
    "WEATHER_SOURCE = LIST[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "OPTIONS = ['per_vineyard', 'per_region']\n",
    "OPTION = OPTIONS[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LAST_YEAR = 2022"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load prices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prices = pd.read_excel('data/prices/prices_per_vineyard.xlsx', index_col=0)\n",
    "prices.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load weather"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "locations = {\n",
    "    \"Médoc\": [45.12, -0.74],\n",
    "    \"Saint-Emilion\": [44.8, -0.192],\n",
    "    \"Pomerol\": [44.925, -0.198],\n",
    "    \"Pessac-Léognan\": [44.81, -0.62],\n",
    "}\n",
    "locations = pd.read_excel(\n",
    "    \"data/gathered_from_internet/vineyard_locations.xlsx\"\n",
    ").set_index(\"name\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_closest_in_list(distances):\n",
    "    min_list = distances[0]\n",
    "    min_index = 0\n",
    "    for i in range(1, len(distances)):\n",
    "        if distances[i] < min_list:\n",
    "            min_list = distances[i]\n",
    "            min_index = i\n",
    "    return min_index\n",
    "\n",
    "\n",
    "def find_closest_points(coords, dataset):\n",
    "    lon, lat = coords\n",
    "    lons = dataset.variables[\"longitude\"][:]\n",
    "    lats = dataset.variables[\"latitude\"][:]\n",
    "\n",
    "    lon_distances = [np.abs(lon - lons[i]) for i in range(len(lons))]\n",
    "    lat_distances = [np.abs(lat - lats[i]) for i in range(len(lats))]\n",
    "    closest_lon = lons[find_closest_in_list(lon_distances)]\n",
    "    closest_lat = lats[find_closest_in_list(lat_distances)]\n",
    "    return closest_lon, closest_lat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "top_wines = list(locations.keys())\n",
    "\n",
    "dict_weather_by_vineyard = {}\n",
    "pheno = pd.read_excel(\"data/pheno/generated_pheno.xlsx\", index_col=0)\n",
    "\n",
    "if WEATHER_SOURCE == \"ERA5\":\n",
    "    from netCDF4 import Dataset\n",
    "    import xarray\n",
    "    import glob\n",
    "\n",
    "    p = xarray.merge([xarray.open_dataset(f) for f in glob.glob(\"data/ERA5/p/*.nc\")])\n",
    "    skt = xarray.merge(\n",
    "        [xarray.open_dataset(f) for f in glob.glob(\"data/ERA5/skt/*.nc\")]\n",
    "    )\n",
    "    evap = xarray.merge(\n",
    "        [xarray.open_dataset(f) for f in glob.glob(\"data/ERA5/evap/*.nc\")]\n",
    "    )\n",
    "    print(p)\n",
    "    print(skt)\n",
    "    print(evap)\n",
    "\n",
    "    dict_chosen_lon_lat = {}\n",
    "    CONVERT_K_TO_C = 273.15\n",
    "    for vineyard in tqdm(top_wines):\n",
    "        print(vineyard)\n",
    "        target_lat, target_lon = (\n",
    "            locations.loc[vineyard, \"latitude\"],\n",
    "            locations.loc[vineyard, \"longitude\"],\n",
    "        )\n",
    "\n",
    "        method = \"linear\"\n",
    "        df_p = (\n",
    "            p[\"tp\"]\n",
    "            .interp(longitude=target_lon, latitude=target_lat, method=method)\n",
    "            .to_dataframe()\n",
    "        )\n",
    "        chosen_lon, chosen_lat = df_p.iloc[0, :][[\"longitude\", \"latitude\"]]\n",
    "        print(\"Chosen coordinates:\", chosen_lon, chosen_lat)\n",
    "        df_p[\"date\"] = pd.to_datetime(df_p.index.date)\n",
    "\n",
    "        df_weather = (\n",
    "            df_p.groupby(\"date\").sum().copy().rename(columns={\"tp\": \"P\"}) * 1000\n",
    "        )\n",
    "\n",
    "        df_skt = (\n",
    "            skt[\"skt\"]\n",
    "            .interp(longitude=target_lon, latitude=target_lat, method=method)\n",
    "            .to_dataframe()\n",
    "        )\n",
    "        df_skt[\"date\"] = pd.to_datetime(df_skt.index.date)\n",
    "\n",
    "        df_weather[\"Tn\"] = df_skt.groupby(\"date\").min()[\"skt\"].copy() - CONVERT_K_TO_C\n",
    "        df_weather[\"Tx\"] = df_skt.groupby(\"date\").max()[\"skt\"].copy() - CONVERT_K_TO_C\n",
    "        df_weather[\"Tm\"] = df_skt.groupby(\"date\").mean()[\"skt\"].copy() - CONVERT_K_TO_C\n",
    "\n",
    "        df_evap = (\n",
    "            evap[\"pev\"]\n",
    "            .interp(longitude=target_lon, latitude=target_lat, method=method)\n",
    "            .to_dataframe()\n",
    "        )\n",
    "        df_evap[\"date\"] = pd.to_datetime(df_evap.index)\n",
    "\n",
    "        df_weather[\"ETP\"] = -df_evap.groupby(\"date\").sum()[\"pev\"].copy() * 1000\n",
    "        df_weather.index = pd.to_datetime(df_weather.index)\n",
    "        dict_weather_by_vineyard[vineyard] = df_weather.copy()\n",
    "        dict_chosen_lon_lat[vineyard] = chosen_lon, chosen_lat\n",
    "\n",
    "elif WEATHER_SOURCE == \"METEO-FRANCE\":\n",
    "    weather_st_em = pd.read_csv(\n",
    "        \"data/weather/generated_weather_st_em.csv\", parse_dates=[\"Date\"]\n",
    "    ).set_index(\"Date\")\n",
    "    weather_medoc = pd.read_csv(\n",
    "        \"data/weather/generated_weather_pauillac.csv\", parse_dates=[\"Date\"]\n",
    "    ).set_index(\"Date\")\n",
    "    weather_graves = pd.read_csv(\n",
    "        \"data/weather/generated_weather_leognan.csv\", parse_dates=[\"Date\"]\n",
    "    ).set_index(\"Date\")\n",
    "    dict_weather_by_region = {\n",
    "        \"Médoc\": weather_medoc,\n",
    "        \"Saint-Emilion\": weather_st_em,\n",
    "        \"Pomerol\": weather_st_em,\n",
    "        \"Pessac-Léognan\": weather_graves,\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "top_wines = list(prices[\"Vineyard\"].unique())\n",
    "region_per_vineyard = {\n",
    "    vineyard: prices.loc[prices[\"Vineyard\"] == vineyard, \"Appellation\"].unique()[0]\n",
    "    for vineyard in top_wines\n",
    "}\n",
    "all_dfs = []\n",
    "for vineyard in tqdm(top_wines):\n",
    "    df_weather_vineyard = create_df(\n",
    "        add_base_features(dict_weather_by_region[region_per_vineyard[vineyard]]), pheno\n",
    "    )\n",
    "    df_elements = df_weather_vineyard.loc[:LAST_YEAR].copy()\n",
    "    interpolated_values = (\n",
    "        prices.loc[prices[\"Vineyard\"] == vineyard, [\"Vintage\", \"Price\"]]\n",
    "        .set_index(\"Vintage\")\n",
    "        .values\n",
    "    )\n",
    "    df_elements[\"0 - Price\"] = np.concatenate([interpolated_values.ravel(), [np.nan, np.nan]]) # empty values for last 2 years\n",
    "    all_dfs.append(df_elements)\n",
    "\n",
    "df_agg = pd.concat(all_dfs, axis=1)\n",
    "features = list(df_weather_vineyard.columns)\n",
    "prod = list(product(top_wines, features + [\"0 - Price\"]))\n",
    "df_agg.columns = pd.MultiIndex.from_tuples(prod)\n",
    "display(df_agg.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_agg['Château Angélus']['0 - Price'].tail(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_agg.to_csv(f\"data/features/features-{WEATHER_SOURCE}_per_vineyard.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variance inflation factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "example='Petrus'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = df_agg[example][\n",
    "    [\"WD: flowering - harvest\", \"P: flowering\", \"DTR: véraison - harvest\"]\n",
    "]\n",
    "\n",
    "variables = list(range(a.shape[1]))\n",
    "vif = [variance_inflation_factor(a.iloc[:, variables].values, ix) for ix in variables]\n",
    "print(\"VIF for the three variables is: \", vif)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Correlations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "features = [feature for feature in features if \n",
    "    'Tn' not in feature\n",
    "    and 'Tx' not in feature\n",
    "    and 'ETP' not in feature\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "correlations = (\n",
    "    df_agg.loc[1990:2013, ([example], features + [\"0 - Price\"])].T.droplevel(0).corr()\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "correlations.loc[:, '0 - Price'].sort_values(ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 14))\n",
    "cmap = sns.diverging_palette(250, 10, as_cmap=True)\n",
    "sns.heatmap(\n",
    "    df_agg.loc[1990:, ([example], features + [\"0 - Price\"])].T.droplevel(0).T.corr(),\n",
    "    cmap=cmap,\n",
    ")\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"views/heatmap.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot temperature evolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tm_st_em = dict_weather_by_region['Saint-Emilion']['Tm']\n",
    "tm_st_em = tm_st_em[tm_st_em.index.month.isin(range(4, 10))].reset_index()\n",
    "tm_yearly_st_em = tm_st_em.groupby(tm_st_em['Date'].dt.year).mean()\n",
    "tm_yearly_st_em.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rc(\"text\", usetex=False)\n",
    "plt.figure(figsize=(10, 5))\n",
    "font = {\"weight\": \"regular\", \"family\": \"serif\", \"size\": 16}\n",
    "plt.grid(True, which=\"both\", axis=\"both\", alpha=0.5)\n",
    "\n",
    "graycolors = sns.mpl_palette(\"Greys_r\", 4)\n",
    "plt.rc(\"font\", **font)\n",
    "plt.rc(\"xtick\", labelsize=\"14\")\n",
    "plt.rc(\"ytick\", labelsize=\"14\")\n",
    "\n",
    "tm_yearly_st_em[\"Year\"] = tm_yearly_st_em.index\n",
    "sns.lineplot(\n",
    "    data=tm_yearly_st_em,\n",
    "    x=\"Year\",\n",
    "    y=\"Tm\",\n",
    "    marker=\"o\",\n",
    "    markersize=6,\n",
    "    linewidth=1.5,\n",
    "    color=\"black\",\n",
    "    markeredgecolor=\"k\",\n",
    "    linestyle=\"dashed\",\n",
    ")\n",
    "plt.gca().lines[0].set_linestyle(\"--\")\n",
    "plt.ylabel(\"Temperature (°C)\")\n",
    "plt.ylim((14.8, 20))\n",
    "save_fig(\"views/temperature\", \"single\")\n",
    "plt.show()"
   ]
  },
  {
   "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
}