Climate effects on soybean.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5b396679-1bc0-4006-b365-45d9fc01ebb0",
   "metadata": {},
   "source": [
    "<center><img src='./images/images/GEOSKOP-logo.jpg', alt=\"Geoskop Logo: Climate Expertise and Earth Inisghts for efficient Decision Making\"></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02ec05c8-c751-46ea-b9f0-b36c79f870c7",
   "metadata": {},
   "source": [
    "# <center>How did climate impact on soy agribusiness in North Dakota?</center>\n",
    "## Analysis of the evolution of different temperature related indexes and their effects in the sustainability of a crop investment"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccf84198-14aa-496b-9128-d923af3f64f9",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "## Table of contents\n",
    "1. [Soybean business in North Dakota, US ](#introduction)\n",
    "2. [Data quality](#paragraph1)\n",
    "3. [Growing Degree Days change](#paragraph2)\n",
    "4. [Temperature impact](#paragraph4)\n",
    "5. [3 days frost](#paragraph5)\n",
    "6. [Conclusions](#paragraph6)\n",
    "7. [Disclaimer](#paragraph7)\n",
    "\n",
    "<br>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "eae57d80-620c-4499-b5ac-e047ca2500d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import interact, widgets\n",
    "from IPython.display import IFrame\n",
    "import base64\n",
    "import io\n",
    "# import calendar\n",
    "# import folium\n",
    "# import numpy as np\n",
    "# from branca.colormap import LinearColormap\n",
    "# import matplotlib.pyplot as plt\n",
    "\n",
    "from pathlib import Path\n",
    "\n",
    "# %matplotlib inline\n",
    "# %matplotlib notebook\n",
    "\n",
    "## How to format markodn text:\n",
    "'''\n",
    "<p><font color=\"red\" face=\"Garamond, Verdana, sans-serif\" size=\"3\">Your formatted text goes here</font></p>\n",
    "'''\n",
    "## Configure notebook to allow local htmls..\n",
    "\n",
    "\n",
    "# Intialize some variables\n",
    "# months = [calendar.month_name[i] for i in range(1,13)]\n",
    "years = [i for i in range(1995,2021)]\n",
    "work_dir = Path('./images/')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a82fe1af-f3bd-4ece-bafd-4970b2a6cc75",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "### Soybean business in North Dakota, US<a name=\"introduction\"></a>\n",
    "In recent years, North Dakota soybean business has grown significantly. This growth is shown in all USDA (US Department of Agriculture) figures.<br>\n",
    "To summarize, in 2002 soybean crops occupied 2.6M acres of the total land in North Dakota state (ND), expanding up to 4.7M acres in 2012, becoming 7.1M acres occupied by soybean crop in 2017. Specifically, this growth was more significant in the north and northeast of ND's state, involving the agricultural regions named North East (Towner, Cavalier, Pembina, Ramsey, Walsh, Nelson and Grand Forks), North Central (Botinneau, Rolette, McHenry, Pierce and Benson) and North West (Divide, Burke, Renville, Williams, Mountrail and Ward).<br>\n",
    "Most western regions went from practically not growing soybeans in 2012 to soy becoming the main crop in 2017.<br>\n",
    "<br>\n",
    "\n",
    "County | Agricultural Region | Acres planted in 2012 | Acres planted in 2017 | % difference \n",
    "--|:---------:|:-----------:|:---------:|:-----------:\n",
    "Buttineau|North Central|42k acres|265k acres|531%\n",
    "Rolette|North Central|27k acres|69k acres|155%\n",
    "McHenry|North Central|32k acres|108k acres|237%\n",
    "Pierce|North Central|94k acres|110k acres|17%\n",
    "Benson|North Central|171k acres|202k acres|18%\n",
    "Towner|North East|68k acres|138k acres|103%\n",
    "Cavalier|North East|59k acres|200k acres|239%\n",
    "Pembina|North East|136k acres|218k acres|60%\n",
    "Ramsey|North East|122k acres|185k acres|51%\n",
    "Walsh|North East|97k acres|162k acres|67%\n",
    "Nelson|North East|117k acres|151k acres|29%\n",
    "Grand Forks|North East|162k acres|245k acres|51%\n",
    "Divide|North West|>1k acres|27k acres|2600%\n",
    "Williams|North West|2k acres|15k acres|650%\n",
    "Bruke|North West|>1 acres|31k acres|3000%\n",
    "Ward|North West|43k acres|214k acres|397%\n",
    "Mountrail|North West|2k acres|67k acres|3250%\n",
    "Renville|North West|17k acres|158k acres|829%\n",
    "TOTAL| | 1193k acres|2565k acres|\n",
    "\n",
    "***\n",
    "\n",
    "Table 1. *Evolution of soybean cultivation in the northern agricultural regions. Credits: USDA*\n",
    "\n",
    "***\n",
    "<br>\n",
    "The table above shows that the total growth between the years 2012 and 2017 in the three northern Agricultural Regions was +1.4M acres, while in all other regions of ND the growth was lower, of >1M acres. \n",
    "<br>\n",
    "Looking further, the trade war between the [US and China](https://www.reuters.com/article/us-china-economy-trade-soybeans-idUSKCN1PJ05M) had a negative impact on soybean farmers. The crop plantation in the year 2019 dropped all around the US due to geopolitical causes, gaining momentum in its recovery nowadays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "6cf4d290-348e-40f8-97d5-931a2f88cc03",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "354186e8a16440e597fa37406a810f81",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Dropdown(description='Year', options={'2002': PosixPath('images/images/nd_beans_2002.png…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "files = work_dir.joinpath('images').glob('*beans*')\n",
    "int_dict = {i.stem.split('_')[-1]:i for i in sorted(files)}\n",
    "\n",
    "@widgets.interact(Year=int_dict)\n",
    "def img2widget(Year):\n",
    "    '''\n",
    "    NOAA_station: Fullfilename (path) as string\n",
    "    '''\n",
    "    file = open(Year, \"rb\")\n",
    "    image = file.read()\n",
    "    extension = Path(Year).suffix[1:]\n",
    "    return widgets.Image(value=image, format=extension,\n",
    "                        width=700,\n",
    "                        height=800)\n",
    "\n",
    "# <br>\n",
    "\n",
    "# ![Soybeans harvested in 2002](./images/images/nd_beans_2002.png \"Soybean harvested acres in 2002. Credits: USDA\")\n",
    "# ***\n",
    "\n",
    "# Figure 1. *Total acres of soybeans harvested in 2002. Credits: USDA*\n",
    "\n",
    "# ***\n",
    "\n",
    "# ![Soybeans harvested in 2017](./images/images/nd_beans_2017.png \"Soybean harvested acres in 2017. Credits: USDA\")\n",
    "# ***\n",
    "\n",
    "# Figure 2. *Total acres of soybeans harvested in 2017. Credits: USDA*\n",
    "\n",
    "# ***"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef4accb3-2e5f-4dbd-a315-53632196d9e2",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "Figure 1. *Total acres of soybeans harvested in 2002, 2012 and 2017. Credits: USDA*\n",
    "\n",
    "***\n",
    "<br>\n",
    "On the whole, this growth is primarily due to market interests and infrastructure development in North Dakota state. At Geoskop, however, **we wonder if all this growth is sustainable from a climatological standpoint, or will the climate turn all these million acres of soybean into an inefficient crop.**<br> To do so, we will rigurously evaluate if there has been any trend in air surface temperature that might affect the yield of soybeans. \n",
    "The study will we done following Geoskop's preanalysis methodology, meeting the following criteria.\n",
    "1. It must be present in the last 25 historical years\n",
    "2. Must occur in the growing season\n",
    "3. It must have an impact on some growth rate of the plant (such as Growing Degree Days, Optimal Growth Temperature, 3 days frost).<br>\n",
    "<br>\n",
    "Let's start!\n",
    "<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "id": "679aa280-9e6b-4d05-89d3-0ae9fdc9bc6e",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"700\"\n",
       "            height=\"600\"\n",
       "            src=\"./images/htmls/map1.html\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x7f34c0130f90>"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# stations = {\n",
    "#     'WBAN:94055':(46.358, -104.250),\n",
    "#     'WBAN:00285':(48.480, -99.236),\n",
    "#     'WBAN:24011':(46.7825,-100.757),\n",
    "#     'WBAN:94011':(48.4167, -101.350),\n",
    "#     'WBAN:00321':(46.167, -98.083),\n",
    "#     'WBAN:00368':(48.941, -97.903)}\n",
    "\n",
    "# mapit = folium.Map(location=[47, -100],\n",
    "#                    zoom_start=6,\n",
    "#                    tiles=None,\n",
    "#                    attr='Geoskop SL',\n",
    "#                    width = '100%',\n",
    "#                    height = '100%',\n",
    "#                    min_zoom = 3,\n",
    "#                    max_zoom = 9,\n",
    "#                    scrollWheelZoom=False,\n",
    "#                    zoom_control = True)\n",
    "# folium.TileLayer('OpenStreetMap', name = 'Legend').add_to(mapit)\n",
    "\n",
    "# visual_tip = 'NOAA Station on click'\n",
    "# for st in stations:\n",
    "#     folium.Marker(location=stations[st],\n",
    "#                   popup=f\"<i>{st}</i>\",\n",
    "#                   icon=folium.features.CustomIcon('https://cdn1.iconfinder.com/data/icons/weather-713/32/Weather_forecast_weather_station_radar-512.png', \n",
    "#                                                   icon_size=(40,40)),\n",
    "#                   tooltip=visual_tip,\n",
    "#                   radius=8).add_to(mapit)\n",
    "    \n",
    "\n",
    "    \n",
    "# aoi = folium.vector_layers.Rectangle([(48.9, -102.4), (48.9, -96.9), (47.64, -96.9), (47.64, -102.4)],\n",
    "#                                      popup=None, fill_color='#43d9de',\n",
    "#                                      tooltip='Area of Analysis')\n",
    "\n",
    "# feat_group = folium.FeatureGroup(name='Area of Interest')\n",
    "# feat_group.add_child(aoi)\n",
    "# mapit.add_child(feat_group)\n",
    "\n",
    "# folium.LayerControl(position = 'topleft', collapsed= False).add_to(mapit)\n",
    "# mapit.save(\"./images/htmls/map1.html\")\n",
    "\n",
    "# Just display the html map already saved\n",
    "IFrame(src='./images/htmls/map1.html', width=700, height=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2378a036-dfc8-4d6c-b17e-3f1bfc3e1b3b",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "Map 1. *NOAA meteorological stations used to validate the skill of the historical model. The map also includes the Area of Interest in the north of ND. Credits: NOAA NCEI*\n",
    "\n",
    "***\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "### Data quality<a name=\"paragraph1\"></a>\n",
    "First of all,we always validate the data we will work with. This validation and selection process provides providing insightful feedback to Geoskop and the client on the accuracy / skill of the analysis or prediction.<br> For this report, only historical data will be used, therefore, we will focus our efforts on validating the skill of our historical model.\n",
    "Thanks to previous projects in the United States, we already know which is the most skilled model, let's validate if the accuracy is still good for temperature in North Dakota by comparing the historical model with real observations from NOAA NCEI meteorological stations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "9780eb20-567f-435a-9b1c-2c6124043bbb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "839018774b8e4bf08050149d88c4893d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Dropdown(description='NOAA_station', options={'WBAN:94055': PosixPath('images/case_plots…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "files = work_dir.joinpath('case_plots').glob('*ncei*')\n",
    "int_dict = {i.stem.split('_')[-1]:i for i in files}\n",
    "\n",
    "@widgets.interact(NOAA_station=int_dict)\n",
    "def plot2widget(NOAA_station):\n",
    "    '''\n",
    "    NOAA_station: Fullfilename (path) as string\n",
    "    '''\n",
    "    file = open(NOAA_station, \"rb\")\n",
    "    image = file.read()\n",
    "    extension = Path(NOAA_station).suffix[1:]\n",
    "    return widgets.Image(value=image, format=extension)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a545222-2eb6-4236-9fea-2bb1b5499eab",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "Figure 2. *Validation of the historical model with NOAA's meteorological stations located in ND state. Credits: NOAA, Copernicus C3S and Geoskop*\n",
    "\n",
    "***\n",
    "<br>\n",
    "In figure 3 we can see that some stations do not have data for specific periods, however, the accuracy of the historical model is high in overall terms, with an average percentage RMSE of 5% (as error estimate). The R-squared shows that historical trends are captured by our historical model, with an average r-squared of 0.98. \n",
    "Once we have validated throughout North Dakota that the historical model is representative, we perform our study-case.<br>\n",
    "The final goal is to verify if all the investment in soybean land in the north area of North Dakota's state (+1.4 million acres) is sustainable from an economic point of view, taking into account the impact of the climate.<br>\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "### Growing Degree Days change<a name=\"paragraph2\"></a>\n",
    "\n",
    "First things first, to analyze if there has been any trend in ND we must look at the big picture. The Growing Degree Days (GDD) or Growing Degree Units (GDU) are key indicators to estimate the growth and development of different crops during the growing season, denominated in units of heat.<br>\n",
    "Broadly speaking, the higher GDDs the better development of the crop.<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "d2c0561e-6cc1-4497-8b09-022d71b2de94",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"700\"\n",
       "            height=\"600\"\n",
       "            src=\"./images/htmls/map2.html\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x7f34c0264e50>"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# def colorize_grid(grid_weights, branca_cmap, normalize_grid = True, colorbar_ticks = 128, cmap_name = 'My plt cmap', i_debug = False):   \n",
    "#     from matplotlib.colors import ListedColormap\n",
    "#     colorbar_ticks_col = colorbar_ticks // (len(branca_cmap.colors) - 1) #number of ticks for each color\n",
    "#     cmap_color = []\n",
    "#     for c in range(1, len(branca_cmap.colors)):\n",
    "#         cmap_color.extend(np.linspace(branca_cmap.colors[c-1], branca_cmap.colors[c], colorbar_ticks_col).tolist())\n",
    "#     input_plt_cmp = ListedColormap(cmap_color, name = cmap_name)\n",
    "#     plt_cmp = plt.cm.get_cmap(input_plt_cmp)\n",
    "#     if normalize_grid:\n",
    "#         normed_data = (grid_weights - np.nanmin(grid_weights)) / (np.nanmax(grid_weights) - np.nanmin(grid_weights))\n",
    "#         return plt_cmp(normed_data)\n",
    "#     else: \n",
    "#         return plt_cmp(grid_weights)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# # Define the size of the map using figure\n",
    "# # fig=Figure(width=250,height=100)\n",
    "# mapit = folium.Map(location=[47, -101],\n",
    "#                    zoom_start=5,\n",
    "#                    tiles=None,\n",
    "#                    attr='Geoskop SL',\n",
    "#                    width = '100%',\n",
    "#                    height = '100%',\n",
    "#                    scrollWheelZoom=False,\n",
    "#                    min_zoom = 3,\n",
    "#                    max_zoom = 9,\n",
    "#                    zoom_control = True)\n",
    "# # Whatever tile but with the name I want\n",
    "# folium.TileLayer('OpenStreetMap',\n",
    "#                  name='Accumulated GDD in the growth season').add_to(mapit)\n",
    "\n",
    "# # Import all rasters\n",
    "# rasters = np.load(work_dir.joinpath('yearly_accum_gdd_ND.npy'))\n",
    "# # Branca cmap\n",
    "# br_cmap = LinearColormap(['#0066ff', '#33ccff', '#ccffff', '#ff9900', '#ffff66'],\n",
    "#                          vmin = np.round(rasters.min()),\n",
    "#                          vmax = np.round(rasters.max()))\n",
    "# show_raster = True\n",
    "# for c in range(rasters.shape[0]):\n",
    "#     grid_colors = colorize_grid(rasters[c], br_cmap,\n",
    "#                                 colorbar_ticks = 18)\n",
    "# #                                 , cmap_name = 'Differece Gradient')\n",
    "#     folium.raster_layers.ImageOverlay(grid_colors, \n",
    "#                                       bounds = [(45.54, -96.2), (49.04, -104.4)], \n",
    "#                                       name = f\"{years[c]}\",\n",
    "#                                       overylay = False,\n",
    "#                                       control = True,\n",
    "#                                       mercator_project = True, \n",
    "#                                       show = show_raster, \n",
    "#                                       opacity = .85).add_to(mapit)\n",
    "#     if c == 0:\n",
    "#         show_raster = False\n",
    "\n",
    "\n",
    "# folium.LayerControl(position = 'topleft', collapsed = False).add_to(mapit)\n",
    "\n",
    "# #specify the min and max values of your data\n",
    "# # colormap = br_cmap.to_step(index=[1700, 1850, 2000, 2150,  2300, 2600, 2800])\n",
    "# br_cmap.caption = 'Accumulated GDD/GDU during the growth season in ND'\n",
    "# br_cmap.add_to(mapit)\n",
    "\n",
    "# mapit.save(\"./images/htmls/map2.html\")\n",
    "\n",
    "# Just display the html map already saved\n",
    "IFrame(src='./images/htmls/map2.html', width=700, height=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4bdf715-106d-491a-9c2d-14fd5f5b32e3",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "Map 2. *Accumulated GDD in ND state for every soy growing season (from 1995 to 2020). Warmer colors indicate higher growing units. Credits: Geoskop SL*\n",
    "\n",
    "***\n",
    "<br>\n",
    "When calculating the GDD, data from our historical model is used. The GDD is calculated for the growing season, including the months of May, June, July, August and September. GDD's formula is calculated with daily maximum and minimum temperature within the range of 86ºF and 50ºF. Is thresholds are exceed, no GDD is added for this day.<br>\n",
    "To compare every growing season in terms of GDD, the accumulated GDD for all the growing season is used.<br>\n",
    "The Map 2 allows us to check over the years if there is any obvious regional trend in ND. Unfortunately, things are never that easy. Therefore, more insights of the analysis are required to visualize a trend.<br>\n",
    "<br>\n",
    "\n",
    "##### Fitting GDDs\n",
    "\n",
    "To better analyze the impact of the climate, the study area is  reduced to the northern region of North Dakota (agricultural regions between Northeast and Northwest, see Map 1). Therefore, the results below are **averages for the entire study area in the growing season**.<br>\n",
    "![GDD AOI fit](./images/case_plots/gdd_trend_north_dakota.png \"GDD yearly trend in the growing season in North Dakota. Credits: Copernicus C3S, Geoskop SL\")\n",
    "***\n",
    "\n",
    "Figure 3. *Yearly accumulated GDD in the growing season. Credits: Copernicus C3S, Geoskop SL*\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbd2b69e-7747-4ae5-97a2-1fab81bb4529",
   "metadata": {},
   "source": [
    "Figure 3 display an uptrend, in which soybeans planted in the area of interest enjoy of 9 GDD more each year on average, improving the yield and efficiency of the crop by external causes. <br>\n",
    "Yet, there are some other climate factors that are critical when evaluating the yield of the crop, let's analize more.<br>\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "### Temperature impact<a name=\"paragraph4\"></a>\n",
    "GDD are useful for the purpose of outlining the optimal growth of the soybean, but considering only the surface temperature, there are other indexes to track. For instance, we might we wondering how many times daily surface temperature come out of the optimal growth threshold.<br>\n",
    "Let's count how many days temperature drops the optimal growth threshold for soybean in our area of interest, which for simplicity, will be the values outside of the range between 12ºC and 32ªC.<br>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "b124e9b4-b17e-43e4-ac1a-d27f794b05d8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d36582bbeac14ae49f9b32be7497b576",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Dropdown(description='Plot', options={'Under optimal growth': PosixPath('images/case_plo…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "files = work_dir.joinpath('case_plots').glob('*temp*')\n",
    "int_dict = {i.stem.split('_')[0]:i for i in sorted(files)[::-1]}\n",
    "\n",
    "@widgets.interact(Plot=int_dict)\n",
    "def img2widget(Plot):\n",
    "    '''\n",
    "    NOAA_station: Fullfilename (path) as string\n",
    "    '''\n",
    "    file = open(Plot, \"rb\")\n",
    "    image = file.read()\n",
    "    extension = Path(Plot).suffix[1:]\n",
    "    return widgets.Image(value=image, format=extension)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94702e29-abf3-47e5-ba4f-11ff756a44ac",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "Figure 4. *Number of times that daily surface temperature drops out of optimal growth thresholds (under 12ºC or over 32ªC). Credits: Copernicus C3S, Geoskop SL*\n",
    "\n",
    "***\n",
    "<br>\n",
    "\n",
    "Results displayed in figure 4 are pretty interesting. From one hand, there is a decreasing trend in which temperature falls below 12ºC, meaning that over the years soybean has had days with a temperature higher than 12ºC, and therefore, the positive trend manifested by the GDD is solid. On the other hand, it seems to be an increasing trend in the number of days with temperature over the optimal growth threshold, yet this trend can be confirmed since data is so disperse.<br>\n",
    "<br>\n",
    "By looking at the \"over optimal growth\" plot in figure 4, two different distribution can be seen (orange and blue color), which complicates the goodness-of-fit test, where the Mean Absolute Percentage Deviation is way to high.<br>\n",
    "To sum up, the number of times temperature is under the optimal growth thresholds tends to be reduced, this factor with the positive trend of GDD implies that the yield of the crop increases, hand in hand with temperature.<br>\n",
    "<br>\n",
    "<br>\n",
    "### 3 days frost<a name=\"paragraph5\"></a>\n",
    "After iterating with farmers and different stakeholders in the agriculture business, we have found that one of the biggest headaches are extreme events of temperature. For example, the month of April 2021 in France has been disastrous for many farmers due to the continuous frosts, totally anomalous.<br>\n",
    "This is why we aim to analyze if the climate have had any impact in this kind of extreme events of the temperature in the north of ND.\n",
    "\n",
    "![3 frost days trend in ND](./images/case_plots/frost_days.png \"Evolution of 3-frost days events. Credits: Copernicus C3S, Geoskop SL\")\n",
    "***\n",
    "\n",
    "Figure 5. *Evolution of 3-frost days events during the growing season. Credits: Copernicus C3S, Geoskop SL*\n",
    "\n",
    "***\n",
    "Yet, there is no evidence of trend regarding 3 days frost in the area of interest, perhaps a minimal decreasing slope which is unimportant, given the timescale it represents.\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "### Conclusions<a name=\"paragraph6\"></a>\n",
    "To summarize, this study explains why it is important to take the climate into account when making an investment decision in new farmland. It is important to highlight that the **climate is not static, but dynamic** and can have an impact on the profitability of an investment.\n",
    "Furthermore, assuming that historical climatology can be representative of the future is a common mistake, that is increasingly relevant since **Climate Change is accelerating the climate dynamics and therefore, lowering the predictability of the climate based only on historical behaviour**.<br>\n",
    "Regarding our study area, the northern region of North Dakota, we demonstrate that the trends of change in temperature have been favorable to soybeans growth. This is a proven fact, but it does not always have to be the case. They are considered favorable conditions for the following reasons:\n",
    "1) The heat accumulated by the crop has increased at a rythm of 9 GDD per year\n",
    "2) The presence of temperatures lower than 12ºC has decreased during the growing season\n",
    "3) No tendencies have been found regarding 3 days frost (extreme events)<br>\n",
    "\n",
    "Finally, the goal of this study is to raise awareness to agricultural decision makers and related stakeholders that the climate can negatively impact to financial projections of an investment. The study contains historical data only, yet in Geoskop we generate climate predictions, up to 40 years ahead, so the impact of the climate can be mitigated and the decision making process takes into account new vectors not considered up to date.\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "### Disclaimer<a name=\"paragraph7\"></a>\n",
    "The present report aims to provide awareness on how the climate can impact (positively or negatively) on the profitability of an agribusiness. In this particular, the evolution of air surface temperature and Growing Degree Days and temperature are analyzed demonstrating the existence of trends in the climate regardless of its randomness.<br>\n",
    "The present report does not contain any climate prediction or project, core technology of Geoskop.<br>\n",
    "The present report does not aim to justify climate change nor to denied since we did not carry a comprehensive survey for that purpose.<br>\n",
    "\n",
    "<br>\n",
    "This report (including any enclosures and attachments) has been prepared for informative purpose. The report shows 3 party data from NOAA NCEI and Copernicus C3S mixed with Geoskop algorithms.\n",
    "This report makes no representations or warranties with respect to this site or its contents and hereby disclaims all warranties, express or implied, including without limitation the implied warranties of merchantability and fitness for a particular purpose, or non-infringement relating to this report, its content or any site to which it is linked.\n",
    "The publisher does not assume and hereby disclaims any liability to any party for any loss, damage, or disruption caused by errors or omissions, whether such errors or omissions result from negligence, accident, or any other cause.\n",
    "<br>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "471404fa-ea28-4a70-96d2-45cb601824a1",
   "metadata": {},
   "source": [
    "<center><img src='./images/images/GEOSKOP-logo.jpg', alt=\"Geoskop Logo: Climate Expertise and Earth Inisghts for efficient Decision Making\"></center>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}