afsis-soil-chem-tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": "true"
   },
   "source": [
    "# Table of Contents\n",
    " <p><div class=\"lev1 toc-item\"><a href=\"#AfSIS-Soil-Chemistry-Dataset-Usage-Tutorial\" data-toc-modified-id=\"AfSIS-Soil-Chemistry-Dataset-Usage-Tutorial-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>AfSIS Soil Chemistry Dataset Usage Tutorial</a></div><div class=\"lev2 toc-item\"><a href=\"#Prerequisities\" data-toc-modified-id=\"Prerequisities-11\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>Prerequisities</a></div><div class=\"lev2 toc-item\"><a href=\"#Downloading-the-Soil-Chemistry-Dataset-from-AWS\" data-toc-modified-id=\"Downloading-the-Soil-Chemistry-Dataset-from-AWS-12\"><span class=\"toc-item-num\">1.2&nbsp;&nbsp;</span>Downloading the Soil Chemistry Dataset from AWS</a></div><div class=\"lev2 toc-item\"><a href=\"#Loading-OPUS-spectra\" data-toc-modified-id=\"Loading-OPUS-spectra-13\"><span class=\"toc-item-num\">1.3&nbsp;&nbsp;</span>Loading OPUS spectra</a></div><div class=\"lev2 toc-item\"><a href=\"#3.-Adding-geographical-references\" data-toc-modified-id=\"3.-Adding-geographical-references-14\"><span class=\"toc-item-num\">1.4&nbsp;&nbsp;</span>3. Adding geographical references</a></div><div class=\"lev2 toc-item\"><a href=\"#4.-Wet-chemistry\" data-toc-modified-id=\"4.-Wet-chemistry-15\"><span class=\"toc-item-num\">1.5&nbsp;&nbsp;</span>4. Wet chemistry</a></div><div class=\"lev2 toc-item\"><a href=\"#Modelling\" data-toc-modified-id=\"Modelling-16\"><span class=\"toc-item-num\">1.6&nbsp;&nbsp;</span>Modelling</a></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# AfSIS Soil Chemistry Dataset - Usage Tutorial¶\n",
    "\n",
    "QED | https://qed.ai | info@qed.ai  \n",
    "\n",
    "This tutorial illustrates the kind of data included in the AfSIS Soil Chemistry Dataset hosted at the Registry of Open Data on AWS, and how it can be used in simple descriptive analyses and machine learning models.\n",
    "\n",
    "## Prerequisities\n",
    "\n",
    "In this tutorial, we require python3.6 with the following packages installed: \n",
    "\n",
    "* To download data:\n",
    "    * aws-cli\n",
    "* To parse and manage datasets:\n",
    "    * brukeropusreader\n",
    "    * pandas \n",
    "    * tqdm\n",
    "* To make plots:\n",
    "    * matplotlib\n",
    "    * folium\n",
    "    * seaborn\n",
    "* To train machine learning models:\n",
    "    * xgboost\n",
    "    * scikit-learn\n",
    "\n",
    "You can install all of this packages with one command:\n",
    "\n",
    "    ! pip install awscli brukeropusreader tqdm pandas matplotlib folium seaborn scikit-learn xgboost           "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Downloading the Soil Chemistry Dataset from AWS\n",
    "\n",
    "Download s3 bucket content with the `aws-cli` command line tool."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! aws s3 sync s3://afsis afsis --no-sign-request"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading OPUS spectra\n",
    "\n",
    "OPUS spectra can be opened with the [brukeropusreader](https://github.com/qedsoftware/brukeropusreader) package. \n",
    "\n",
    "The function `brukeropusreader.read_file` parses the binaries and returns a data structure containing information about the wave numbers, absorbance spectra, and file metadata. \n",
    "You can change path in variable `SPECTRA_PATH` to choose which files you want to work with.\n",
    "In this tutorial we will use `Bruker_MPA/*` spectra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "SPECTRA_DIR = 'Bruker_MPA/*'\n",
    "\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from pathlib import Path\n",
    "from brukeropusreader import read_file\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "\n",
    "AFSIS_PATH = Path('afsis/2009-2013/Dry_Chemistry/ICRAF/')\n",
    "names = []\n",
    "spectra = []\n",
    "\n",
    "for path in tqdm(AFSIS_PATH.glob(SPECTRA_DIR)):\n",
    "    if path.is_file():\n",
    "        spect_data = read_file(path)\n",
    "        spectra.append(spect_data[\"AB\"])\n",
    "        names.append(path.stem)\n",
    "wave_nums = spect_data.get_range()\n",
    "\n",
    "    \n",
    "column_names = ['{:.3f}'.format(x) for x in wave_nums]\n",
    "spectra_df = pd.DataFrame(spectra, index=names, columns=column_names)\n",
    "spectra_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we plot a few of the spectra. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = spectra_df.sample(3).T.plot(figsize=(12.5, 10), title='AfSIS Spectra')\n",
    "fig.set_xlabel(\"Wavelengths\")\n",
    "fig.set_ylabel(\"Absorbance\")\n",
    "fig.legend(bbox_to_anchor=(1.2, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Adding geographical references\n",
    "\n",
    "The AfSIS Soil Chemistry Dataset contains georeferences for each spectra. \n",
    "You can use them to locate every measurement on a map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GEOREFS_FILE = 'afsis/2009-2013/Georeferences/georeferences.csv'\n",
    "georefs_df = pd.read_csv(GEOREFS_FILE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import folium\n",
    "m = folium.Map(location=[-3.5, 35.6], tiles=\"stamentoner\", zoom_start=3)\n",
    " \n",
    "for _, row in georefs_df.iterrows():\n",
    "    if row[['Latitude', 'Longitude']].notnull().all():\n",
    "        folium.Marker([row['Latitude'], \n",
    "                       row['Longitude']], \n",
    "                      popup=row['SSN']\n",
    "                     ).add_to(m)\n",
    "\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![afsis_map](img/folium.png \"AfSIS Map\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also count the number of samples in each country:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.value_counts(georefs_df['Country']).plot.bar(title='Measurements per country')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The georeferences table also contains many other pieces of metadata, as shown below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "georefs_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Wet chemistry\n",
    "\n",
    "The AfSIS Soil Chemistry dataset contains wet chemistry data taken at each sampling location.\n",
    "\n",
    "You can select which wet chemistry dataset you want to work with in `WET_CHEM_PATH`variable and which elements you want to plot in `elements` array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "WET_CHEM_PATH = 'afsis/2009-2013/Wet_Chemistry/CROPNUTS/Wet_Chemistry_CROPNUTS.csv'\n",
    "elements = ['M3 Ca', 'M3 K', 'M3 Al']\n",
    "columns_to_load = elements + ['SSN']\n",
    "\n",
    "\n",
    "wet_chem_df = pd.read_csv(WET_CHEM_PATH, index_col='SSN', usecols=columns_to_load)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can plot the elemental distributions of all loaded samples using the `seaborn` library. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sn\n",
    "\n",
    "for e in elements:\n",
    "  sn.distplot(wet_chem_df[e]).set_title(e)\n",
    "  plt.show()    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also check for correlations between any of the wet chemistry pairs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sn.pairplot(wet_chem_df[elements], diag_kind='kde')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modelling\n",
    "\n",
    "Together, the wet chemistry and dry chemistry measurements can serve as a training dataset for performing spectral inference.\n",
    "Namely, we can prepare algorithmic data models to estimate the elemental concentrations that wet chemistry would have produced, based on spectral measurements. Below is a very simple example of how this can work. \n",
    "\n",
    "You can set different wet chemistry path in variable `WET_CHEM_PATH` and different element in variable `Y_COLUMN`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "WET_CHEM_PATH = 'afsis/2009-2013/Wet_Chemistry/CROPNUTS/Wet_Chemistry_CROPNUTS.csv'\n",
    "Y_COLUMN = 'M3 Ca'\n",
    "\n",
    "from sklearn.ensemble import ExtraTreesRegressor\n",
    "from sklearn.model_selection import cross_validate\n",
    "from sklearn.model_selection import KFold\n",
    "\n",
    "import numpy as np\n",
    "from xgboost import XGBRegressor\n",
    "wet_chem_df = pd.read_csv(WET_CHEM_PATH, index_col='SSN')\n",
    "data_df = spectra_df.merge(wet_chem_df, left_index=True, right_index=True)\n",
    "\n",
    "print(f'Training model for {Y_COLUMN}.')\n",
    "print(f'{data_df.shape[0]} samples available.')\n",
    "\n",
    "\n",
    "y = data_df[Y_COLUMN]\n",
    "y_mask = y.notnull()\n",
    "y = y[y_mask]\n",
    "y = np.log1p(y)\n",
    "\n",
    "x = data_df[column_names]\n",
    "x = x[y_mask]\n",
    "x = np.apply_along_axis(np.gradient, 1, x)\n",
    "\n",
    "kfold = KFold(shuffle=True, random_state=0, n_splits=4)\n",
    "model = XGBRegressor(n_estimators=500, min_child_weight=20, n_jobs=-1, objective='reg:squarederror')\n",
    "    \n",
    "cv_result = cross_validate(model, x, y, scoring='r2', cv=kfold)['test_score']\n",
    "print(f'CV r^2 score: {np.mean(cv_result)}')"
   ]
  }
 ],
 "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.6.9"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "136px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": true,
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}