SOCAT_20082010_AtlOcean_ADS_05062024.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bc8d30e3",
   "metadata": {},
   "source": [
    "## Observation system simulation experiments in the Atlantic Ocean for enhanced surface ocean $pCO_2$ reconstructions.\n",
    "\n",
    "_Anna Denvil-Sommer (NCAS, University of Reading), Pier Luigi Vidale (NCAS, University of Reading), Laura Cimoli (ICCS, University of Cambridge)_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a339f38a",
   "metadata": {},
   "source": [
    "1. [Introduction](#1.-Introduction)\n",
    "2. [Data](#2.-Data)\n",
    "2. [Features description and preprocessing](#3.-Features-description-and-preprocessing)\n",
    "4. [Exercise ](#4.-Exercise)\n",
    " \n",
    " a. [Read and visualise data](#a.-Read-and-visualise-data)\n",
    " \n",
    " b. [Define the neural network model, train the model and visualise loss function](#b.-Define-the-neural-network-model,-train-the-model-and-visualise-loss-function)\n",
    " \n",
    " c. [Validate the model](#c.-Validate-the-model)\n",
    " \n",
    " d. [Estimate and visualise model's statistic](#d.-Estimate-and-visualise-model's-statistic)\n",
    " \n",
    "5. [Conclusion](#5.-Conclusion)\n",
    "6. [References](#6.-References)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9f67cdb",
   "metadata": {},
   "source": [
    "<img src=\"SummerSchool2024_pCO2_What_Is_It_About.png\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "124a19de",
   "metadata": {},
   "source": [
    "## 1. Introduction\n",
    "\n",
    "To monitor anthropogenic $CO_2$ fluxes to the atmosphere with a good accuracy it is important to reduce the uncertainty in estimates of air-sea $CO_2$ exchanges.\n",
    "\n",
    "Air-sea $CO_2$ fluxes depend on the difference in $CO_2$ partial pressure between the atmosphere and the ocean. The atmosphere being relatively well mixed, the main source of uncertainty in estimating air-sea $CO_2$ fluxes is related to the estimation of ocean surface $CO_2$ partial pressure: $pCO_2$. In practice, $pCO_2$ is estimated from sea surface fugacity ($fCO_2$), which expresses the tendency of $CO_2$ to escape from the water; the partial pressure differs from fugacity by a factor of about 0.996, due to a slight non-ideal gas behavior of $CO_2$.\n",
    "\n",
    "The ocean is a major sink of anthropogenic $CO_2$. For the period 2010–2019 the ocean uptake was 2.5 ± 0.6 GtC yr−1 with a strong intensification (from 1.9 to 3.1 GtC yr−1). Overall, the ocean is contributing to the sequestration of about 30% of the cumulative atmospheric anthropogenic $CO_2$. Most observations contributing to the Surface Ocean CO2 Atlas (SOCAT) (Bakker et al., 2016) are still obtained by surface water sampling systems on board volunteering observing ships. The data density is not homogenous, with southern latitudes being less well sampled in space and time, see figure below. \n",
    "\n",
    "<img src=\"SOCAT.png\" width=\"70%\">\n",
    "<div style=\"text-align: center\">Locations of moorings and tracks of ships and drifters for all data in SOCAT version 2023, observational dataset for 1957-2022.</div>\n",
    "\n",
    "In contrast, measurements of oceanic *state* parameters (like temperature, salinity) are usually much denser in space and time (ARGO datasets, satellite remote sensing data). This has motivated the development of methods that use the statistical relationship between *physical* and *biogeochemical* parameters for extrapolating existing $pCO_2$ data in space and time. In this context, Machine Learning (ML) approaches have gained a lot of attention over recent years. However, sparse data coverage and the lack of observations covering the full seasonal cycle challenge mapping methods and result in noisy reconstructions of surface ocean $pCO_2$ and disagreements between different models. \n",
    "\n",
    "How can we optimise $pCO_2$ sampling strategies and collect observations where they are most needed? In this exercise we will show you how using ML approaches and output of ocean physical-biogeochemical models we can explore design options for a future Atlantic-scale observing system that would optimally combine data streams from various platforms and contribute to reduce the bias in reconstructed surface ocean $pCO_2$ fields and sea–air $CO_2$ fluxes. \n",
    "The use of NEMO/PISCES physical-biogeochemical model's outputs allows us create Argo floats pseudo-observations (i.e. synthetic observations) and identify the most critical regions that require more observations to provide a robust data-driven product based on ML approach. \n",
    "\n",
    "This exercise is based on the work Denvil-Sommer, A., Gehlen, M., and Vrac, M.: Observation system simulation experiments in the Atlantic Ocean for enhanced surface ocean pCO2 reconstructions, Ocean Sci., 17, 1011–1030, https://doi.org/10.5194/os-17-1011-2021, 2021."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a801b604",
   "metadata": {},
   "source": [
    "## 2. Data\n",
    "\n",
    "Three observing platforms were selected for the exercise: (1) volunteering observing ships providing in situ measurements of surface ocean $CO_2$ fugacity ($fCO_2$) (SOCAT database), (2) moorings (SOCAT database), and (3) profilers (Argo). These observations form the dataset of geographical and temporal positions for our experiments. \n",
    "Biogeochemical Argo floats are increasingly equipped with pH sensors, allowing computing $pCO_2$ from pH and SST-based alkalinity. For the design experiments, we considered distributions of physical Argo floats (2008–2011) from Gasparin et al. (2019) and supposed that they were equipped with $pCO_2$ sensors.\n",
    "\n",
    "We have three datasets to estimate the role of Argo floats in accuracy improvement of $pCO_2$ reconstruction using Feed-Forward Neural Network (FFNN): \n",
    "- The first experiment is based on SOCAT data only (i.e. moorings, ship track and drifters data only).\n",
    "- The second experiment is based on SOCAT data and adding synthetic data assuming that 25% of Argo floats over the whole Atlantic basin were equipped with $pCO_2$ sensors <distribution over the whole Atlantic basin>; the 25% choice comes from what was proposed in the original article Denvil-Sommer et al., 2021.\n",
    "- The third experiment is based on SOCAT data and adding 25% of Argo float distribution over the South part of Atlantic basin as synthetic data.\n",
    "  \n",
    "For all experiments, we reconstruct only February month for the period 2008-2010.\n",
    "\n",
    "In this exercise we use the data from NEMO/PISCES physical-biogeochemical model at 5-day resolution. This configuration of the NEMO framework was implemented on a global tripolar grid. \n",
    "\n",
    "The geographical and time positions identified from the observational data are used to create pseudo-observations by sub-sampling NEMO/PISCES model output at sites of real-world observations. Thus, the positions of SOCAT, Argo floats and mooring stations were chosen over 5 d centred on the NEMO/PISCES date and sub-sampled on the model grid.\n",
    "\n",
    "In Observational System Simulation Experiments the model represents a \"truth\" and used to estimate the accuracy of ML approach. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0eb29c3c",
   "metadata": {},
   "source": [
    "## 3. Features description and preprocessing\n",
    "\n",
    "Two sets of data are needed to test the machine learning method: a set of targets and a set of drivers. The drivers represent the input variables to the ML method (here the biological, chemical and environmental variables). The targets represent the variables we are trying to reconstruct.\n",
    "\n",
    "### Target \n",
    "\n",
    "The target of FFNN is NEMO/PISCES $pCO_2$. \n",
    "\n",
    "### Predictors \n",
    "\n",
    "The standard set of variables known to represent the *physical* and *biogeochemical* drivers of surface ocean $pCO_2$ are sea surface salinity (SSS), sea surface temperature (SST), mixed layer depth (MLD), chlorophyll a concentration (CHL), sea surface height (SSH) and the atmospheric $CO_2$ mole fraction ($x CO_{2, atm}$ ) (Takahashi et al., 2009; Landschützer et al., 2013). These variables and their anomalies are proposed as input variables (or predictors) for training ML algorithms. Additionally, we use geographical coordinates, latitude and longitude, as predictors.\n",
    "\n",
    "All data are taken from the NEMO/PISCES model at the position of SOCAT or ARGO float in February 2008-2010. \n",
    "\n",
    "### Normalization\n",
    "\n",
    "We log-transformed MLD and CHL data because of their skewed distribution. It is also worth noting that all datasets (including target) need to be normalized (i.e., centered to zero-mean and reduced to unit standard deviation). Normalization ensures that all predictors fall within a comparable range and avoids giving more weight to predictors with large variability ranges (Kallache et al., 2011), example:\n",
    "\n",
    "$$SSS_n =\\frac{SSS-\\overline{SSS}}{std(SSS)}.$$\n",
    "\n",
    "Latitude and longitude are also normalized using following way:\n",
    "\n",
    "$$lat_n = sin(lat * \\pi/180)$$\n",
    "\n",
    "$$long_{n,1} = sin(long * \\pi/180)$$\n",
    "\n",
    "$$long_{n,2} = cos(long * \\pi/180)$$\n",
    "\n",
    "\n",
    "All data sets are normalized and ready to be used in this exercise. We also provide mean and std values for $pCO_2$ to visualise the $pCO_2$ fields at the end of exercise. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2c19f31",
   "metadata": {},
   "source": [
    "## 4. Exercise \n",
    "\n",
    "The following Python libraries are required to run the exercise successfully. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91219aa5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.optim as optim\n",
    "from torch.utils.data import DataLoader, TensorDataset\n",
    "from sklearn.model_selection import train_test_split\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "168ec9b3",
   "metadata": {},
   "source": [
    "## a. Read and visualise data\n",
    "\n",
    "Before we define a Feed-Forward Neural Network we will have a look on the distribution of available datasets. The data are in .csv files and can be easily read in panda dataframes. \n",
    "We have two sets of data for each case: training and evaluation data. Training datasets are used to train the model, evaluation datasets are also used during training to evaluate the model at each iterative model step. These data were chosen regularly in time and space: every 4th grid point was kept for evaluation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "152259df",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare data\n",
    "# Data_train, data_eval, pCO2_list_train, and pCO2_list_eval are pandas dataframes\n",
    "#Data for SOCAT only experiments\n",
    "data_train_SOCAT      = pd.read_csv('TrainingSetVarn_SOCAT_20082010_Feb.csv')  # Load training data here\n",
    "data_eval_SOCAT       = pd.read_csv('EvaluationSetVarn_SOCAT_20082010_Feb.csv')  # Load evaluation data here\n",
    "pCO2_list_train_SOCAT = pd.read_csv('TrainingSetpCO2_SOCAT_20082010_Feb.csv')  # Load training targets here\n",
    "pCO2_list_eval_SOCAT  = pd.read_csv('EvaluationSetpCO2_SOCAT_20082010_Feb.csv')  # Load evaluation targets here\n",
    "data_train_SOCAT = data_train_SOCAT.drop(columns=['Unnamed: 0'])\n",
    "data_eval_SOCAT = data_eval_SOCAT.drop(columns=['Unnamed: 0'])\n",
    "\n",
    "#Data for SOCAT + ARGO 25% experiments\n",
    "data_train_SOCAT_Argo25      = pd.read_csv('TrainingSetVarn_SOCAT_Argo25_20082010_Feb.csv')  # Load training data here\n",
    "data_eval_SOCAT_Argo25       = pd.read_csv('EvaluationSetVarn_SOCAT_Argo25_20082010_Feb.csv')  # Load evaluation data here\n",
    "pCO2_list_train_SOCAT_Argo25 = pd.read_csv('TrainingSetpCO2_SOCAT_Argo25_20082010_Feb.csv')  # Load training targets here\n",
    "pCO2_list_eval_SOCAT_Argo25  = pd.read_csv('EvaluationSetpCO2_SOCAT_Argo25_20082010_Feb.csv')  # Load evaluation targets here\n",
    "data_train_SOCAT_Argo25 = data_train_SOCAT_Argo25.drop(columns=['Unnamed: 0'])\n",
    "data_eval_SOCAT_Argo25 = data_eval_SOCAT_Argo25.drop(columns=['Unnamed: 0'])\n",
    "\n",
    "#Data for SOCAT + ARGO 25% in the South part of Atlantic basin experiments\n",
    "data_train_SOCAT_Argo25S      = pd.read_csv('TrainingSetVarn_SOCAT_Argo25_South_20082010_Feb.csv')  # Load training data here\n",
    "data_eval_SOCAT_Argo25S       = pd.read_csv('EvaluationSetVarn_SOCAT_Argo25_South_20082010_Feb.csv')  # Load evaluation data here\n",
    "pCO2_list_train_SOCAT_Argo25S = pd.read_csv('TrainingSetpCO2_SOCAT_Argo25_South_20082010_Feb.csv')  # Load training targets herepCO2_list_eval_SOCAT_Argo25S  = pd.read_csv('EvaluationSetpCO2_SOCAT_Argo25_South_20082010_Feb.csv')  # Load evaluation targets here\n",
    "pCO2_list_eval_SOCAT_Argo25S  = pd.read_csv('EvaluationSetpCO2_SOCAT_Argo25_South_20082010_Feb.csv')  # Load evaluation targets here\n",
    "data_train_SOCAT_Argo25S = data_train_SOCAT_Argo25S.drop(columns=['Unnamed: 0'])\n",
    "data_eval_SOCAT_Argo25S = data_eval_SOCAT_Argo25S.drop(columns=['Unnamed: 0'])\n",
    "\n",
    "#Put data in tensor format to further use in ML model\n",
    "data_train_tensor_SOCAT      = torch.tensor(data_train_SOCAT.values, dtype=torch.float32)\n",
    "data_eval_tensor_SOCAT       = torch.tensor(data_eval_SOCAT.values, dtype=torch.float32)\n",
    "pCO2_list_train_tensor_SOCAT = torch.tensor(pCO2_list_train_SOCAT['pCO2_n'].values, dtype=torch.float32)\n",
    "pCO2_list_eval_tensor_SOCAT  = torch.tensor(pCO2_list_eval_SOCAT['pCO2_n'].values, dtype=torch.float32)\n",
    "\n",
    "train_dataset_SOCAT = TensorDataset(data_train_tensor_SOCAT, pCO2_list_train_tensor_SOCAT)\n",
    "eval_dataset_SOCAT  = TensorDataset(data_eval_tensor_SOCAT, pCO2_list_eval_tensor_SOCAT)\n",
    "\n",
    "train_loader_SOCAT = DataLoader(train_dataset_SOCAT, batch_size=20, shuffle=True)\n",
    "eval_loader_SOCAT  = DataLoader(eval_dataset_SOCAT, batch_size=20, shuffle=False)\n",
    "\n",
    "\n",
    "data_train_tensor_SOCAT_Argo25      = torch.tensor(data_train_SOCAT_Argo25.values, dtype=torch.float32)\n",
    "data_eval_tensor_SOCAT_Argo25       = torch.tensor(data_eval_SOCAT_Argo25.values, dtype=torch.float32)\n",
    "pCO2_list_train_tensor_SOCAT_Argo25 = torch.tensor(pCO2_list_train_SOCAT_Argo25['pCO2_n'].values, dtype=torch.float32)\n",
    "pCO2_list_eval_tensor_SOCAT_Argo25  = torch.tensor(pCO2_list_eval_SOCAT_Argo25['pCO2_n'].values, dtype=torch.float32)\n",
    "\n",
    "train_dataset_SOCAT_Argo25 = TensorDataset(data_train_tensor_SOCAT_Argo25, pCO2_list_train_tensor_SOCAT_Argo25)\n",
    "eval_dataset_SOCAT_Argo25  = TensorDataset(data_eval_tensor_SOCAT_Argo25, pCO2_list_eval_tensor_SOCAT_Argo25)\n",
    "\n",
    "train_loader_SOCAT_Argo25 = DataLoader(train_dataset_SOCAT_Argo25, batch_size=20, shuffle=True)\n",
    "eval_loader_SOCAT_Argo25  = DataLoader(eval_dataset_SOCAT_Argo25, batch_size=20, shuffle=False)\n",
    "\n",
    "\n",
    "data_train_tensor_SOCAT_Argo25S      = torch.tensor(data_train_SOCAT_Argo25S.values, dtype=torch.float32)\n",
    "data_eval_tensor_SOCAT_Argo25S       = torch.tensor(data_eval_SOCAT_Argo25S.values, dtype=torch.float32)\n",
    "pCO2_list_train_tensor_SOCAT_Argo25S = torch.tensor(pCO2_list_train_SOCAT_Argo25S['pCO2_n'].values, dtype=torch.float32)\n",
    "pCO2_list_eval_tensor_SOCAT_Argo25S  = torch.tensor(pCO2_list_eval_SOCAT_Argo25S['pCO2_n'].values, dtype=torch.float32)\n",
    "\n",
    "train_dataset_SOCAT_Argo25S = TensorDataset(data_train_tensor_SOCAT_Argo25S, pCO2_list_train_tensor_SOCAT_Argo25S)\n",
    "eval_dataset_SOCAT_Argo25S  = TensorDataset(data_eval_tensor_SOCAT_Argo25S, pCO2_list_eval_tensor_SOCAT_Argo25S)\n",
    "\n",
    "train_loader_SOCAT_Argo25S = DataLoader(train_dataset_SOCAT_Argo25S, batch_size=20, shuffle=True)\n",
    "eval_loader_SOCAT_Argo25S  = DataLoader(eval_dataset_SOCAT_Argo25S, batch_size=20, shuffle=False)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85c38764",
   "metadata": {},
   "source": [
    "Now we will visualise the distribution of three datasets for training and then three datasets for evaluation.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "178a7b0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a figure with three subplots\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 5), subplot_kw={'projection': ccrs.PlateCarree()})\n",
    "\n",
    "# Plot 1: Training Data\n",
    "ax = axes[0]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_list_train_SOCAT['lon'], pCO2_list_train_SOCAT['lat'], c=pCO2_list_train_SOCAT['pCO2'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree())\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('pCO2 Train Data Distribution, SOCAT only')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Plot 2: Evaluation Data\n",
    "ax = axes[1]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_list_train_SOCAT_Argo25['lon'], pCO2_list_train_SOCAT_Argo25['lat'], c=pCO2_list_train_SOCAT_Argo25['pCO2'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree())\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('pCO2 Train Data Distribution, SOCAT + 25% Argo')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Plot 3: Validation Data\n",
    "ax = axes[2]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_list_train_SOCAT_Argo25S['lon'], pCO2_list_train_SOCAT_Argo25S['lat'], c=pCO2_list_train_SOCAT_Argo25S['pCO2'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree())\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('pCO2 Train Data Distribution, SOCAT + 25% Argo South')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Adjust layout to prevent overlap\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e051f460",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a figure with three subplots\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 5), subplot_kw={'projection': ccrs.PlateCarree()})\n",
    "\n",
    "# Plot 1: Training Data\n",
    "ax = axes[0]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_list_eval_SOCAT['lon'], pCO2_list_eval_SOCAT['lat'], c=pCO2_list_eval_SOCAT['pCO2'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree())\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('pCO2 Evaluation Data Distribution, SOCAT only')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Plot 2: Evaluation Data\n",
    "ax = axes[1]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_list_eval_SOCAT_Argo25['lon'], pCO2_list_eval_SOCAT_Argo25['lat'], c=pCO2_list_eval_SOCAT_Argo25['pCO2'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree())\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('pCO2 Evaluation Data Distribution, SOCAT + 25% Argo')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Plot 3: Validation Data\n",
    "ax = axes[2]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_list_eval_SOCAT_Argo25S['lon'], pCO2_list_eval_SOCAT_Argo25S['lat'], c=pCO2_list_eval_SOCAT_Argo25S['pCO2'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree())\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('pCO2 Evaluation Data Distribution, SOCAT + 25% Argo South')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Adjust layout to prevent overlap\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eed79332",
   "metadata": {},
   "source": [
    "## Question:\n",
    "\n",
    "What would you conclude about the future accuracy of ML methods based on the distribution of available data? "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6db8577d",
   "metadata": {},
   "source": [
    "## b. Define the neural network model, train the model and visualise loss function\n",
    "\n",
    "## Feed-Forward Neural Network\n",
    "\n",
    "To adjust the number of FFNN parameters/weights we followed the empirical rule that suggests limiting the number of parameters to the number of training data points divided by 10 to avoid overfitting (Amari et al., 1997). \n",
    "The FFNN has four layers (two hidden layers). The input layer has 15 input nodes (15 predictors) and 20 output nodes that represent the input for the first hidden layer. The first hidden layer has 25 output nodes, and the second hidden layer has 10 output nodes. \n",
    "\n",
    "We use Early Stopping to avoid overfitting. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91702d6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the neural network model\n",
    "class FeedForwardNN(nn.Module):\n",
    "    def __init__(self):\n",
    "        super(FeedForwardNN, self).__init__()\n",
    "        self.fc1 = nn.Linear(15, 20)\n",
    "        self.tanh1 = nn.Tanh()\n",
    "        self.fc2 = nn.Linear(20, 25)\n",
    "        self.tanh2 = nn.Tanh()\n",
    "        self.fc3 = nn.Linear(25, 10)\n",
    "        self.tanh3 = nn.Tanh()\n",
    "        self.fc4 = nn.Linear(10, 1)\n",
    "        self.linear = nn.Identity()\n",
    "\n",
    "    def forward(self, x):\n",
    "        x = self.tanh1(self.fc1(x))\n",
    "        x = self.tanh2(self.fc2(x))\n",
    "        x = self.tanh3(self.fc3(x))\n",
    "        x = self.linear(self.fc4(x))\n",
    "        return x\n",
    "\n",
    "# Early stopping class\n",
    "class EarlyStopping:\n",
    "    def __init__(self, patience=30, verbose=False, delta=0):\n",
    "        self.patience = patience\n",
    "        self.verbose = verbose\n",
    "        self.counter = 0\n",
    "        self.best_score = None\n",
    "        self.early_stop = False\n",
    "        self.val_loss_min = np.Inf\n",
    "        self.delta = delta\n",
    "\n",
    "    def __call__(self, val_loss, model):\n",
    "        score = -val_loss\n",
    "\n",
    "        if self.best_score is None:\n",
    "            self.best_score = score\n",
    "            self.save_checkpoint(val_loss, model)\n",
    "        elif score < self.best_score + self.delta:\n",
    "            self.counter += 1\n",
    "            if self.verbose:\n",
    "                print(f'EarlyStopping counter: {self.counter} out of {self.patience}')\n",
    "            if self.counter >= self.patience:\n",
    "                self.early_stop = True\n",
    "        else:\n",
    "            self.best_score = score\n",
    "            self.save_checkpoint(val_loss, model)\n",
    "            self.counter = 0\n",
    "\n",
    "    def save_checkpoint(self, val_loss, model):\n",
    "        '''Saves model when validation loss decreases.'''\n",
    "        if self.verbose:\n",
    "            print(f'Validation loss decreased ({self.val_loss_min:.6f} --> {val_loss:.6f}).  Saving model ...')\n",
    "        torch.save(model.state_dict(), 'checkpoint.pt')\n",
    "        self.val_loss_min = val_loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a41c40e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the model, loss function, and optimizer\n",
    "model_SOCAT = FeedForwardNN()\n",
    "criterion = nn.MSELoss()\n",
    "optimizer = optim.RMSprop(model_SOCAT.parameters(), lr=0.001, alpha=0.9, eps=1e-08, weight_decay=0.0)\n",
    "\n",
    "# Calculate the number of parameters\n",
    "total_params = sum(p.numel() for p in model_SOCAT.parameters())\n",
    "print(f'Total number of parameters: {total_params}')\n",
    "\n",
    "# Print a detailed summary of the parameters\n",
    "for name, param in model_SOCAT.named_parameters():\n",
    "    if param.requires_grad:\n",
    "        print(f'{name}: {param.numel()} parameters')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5fe33bf7",
   "metadata": {},
   "source": [
    "## Questions:\n",
    "\n",
    "How many parameters in total does the model have? Does the number correspond to the empirical rule for an optimal number of parameters?\n",
    "\n",
    "After how many iterations will the model stop if there is not significant improvement during the training? "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6917b5aa",
   "metadata": {},
   "source": [
    "Now we'll train the model on the dataset that includes only SOCAT data.\n",
    "To follow the evaluation of loss function during the training two list, 'train_losses_SOCAT' and 'val_losses_SOCAT', are created to store the values from train and evaluation datasets, respectively. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ac29097",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Training the model\n",
    "num_epochs = 1600\n",
    "early_stopping = EarlyStopping(patience=30, verbose=True)\n",
    "\n",
    "# Lists to store loss values\n",
    "train_losses_SOCAT = []\n",
    "val_losses_SOCAT = []\n",
    "\n",
    "for epoch in range(num_epochs):\n",
    "    model_SOCAT.train()\n",
    "    train_loss_SOCAT = 0.0\n",
    "    for batch_data, batch_target in train_loader_SOCAT:\n",
    "        optimizer.zero_grad()\n",
    "        outputs = model_SOCAT(batch_data)\n",
    "        loss = criterion(outputs, batch_target.view(-1, 1))\n",
    "        loss.backward()\n",
    "        optimizer.step()        \n",
    "        train_loss_SOCAT += loss.item()\n",
    "            \n",
    "    train_loss_SOCAT /= len(train_loader_SOCAT)\n",
    "    train_losses_SOCAT.append(train_loss_SOCAT)\n",
    "\n",
    "    model_SOCAT.eval()\n",
    "    val_loss_SOCAT = 0.0\n",
    "    with torch.no_grad():\n",
    "        for batch_data, batch_target in eval_loader_SOCAT:\n",
    "            outputs = model_SOCAT(batch_data)\n",
    "            loss = criterion(outputs, batch_target.view(-1, 1))\n",
    "            val_loss_SOCAT += loss.item()\n",
    "\n",
    "    val_loss_SOCAT /= len(eval_loader_SOCAT)\n",
    "    val_losses_SOCAT.append(val_loss_SOCAT)\n",
    "    print(f'Epoch {epoch+1}, Validation Loss: {val_loss_SOCAT:.6f}')\n",
    "\n",
    "    early_stopping(val_loss_SOCAT, model_SOCAT)\n",
    "    if early_stopping.early_stop:\n",
    "        print(\"Early stopping\")\n",
    "        break\n",
    "\n",
    "# Load the best model\n",
    "model_SOCAT.load_state_dict(torch.load('checkpoint.pt'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48db7eb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting the loss\n",
    "plt.figure(figsize=(10,5))\n",
    "plt.plot(train_losses_SOCAT, label='Training Loss')\n",
    "plt.plot(val_losses_SOCAT, label='Evaluation Loss')\n",
    "plt.xlabel('Epochs')\n",
    "plt.ylabel('Loss')\n",
    "plt.legend()\n",
    "plt.title('Training and Evaluation Loss Over Epochs')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebce8ac0",
   "metadata": {},
   "source": [
    "## Question:\n",
    "\n",
    "How would you interpret this figure of loss function? What would you say about overfitting? "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbb579a6",
   "metadata": {},
   "source": [
    "Now we'll repeat same training procedure for two other datasets that include data from Argo floats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bfc0c5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the model, loss function, and optimizer\n",
    "model_SOCAT_Argo25 = FeedForwardNN()\n",
    "criterion = nn.MSELoss()\n",
    "optimizer = optim.RMSprop(model_SOCAT_Argo25.parameters(), lr=0.001, alpha=0.9, eps=1e-08, weight_decay=0.0)\n",
    "\n",
    "# Training the model\n",
    "num_epochs = 1600\n",
    "early_stopping = EarlyStopping(patience=30, verbose=True)\n",
    "\n",
    "# Lists to store loss values\n",
    "train_losses_SOCAT_Argo25 = []\n",
    "val_losses_SOCAT_Argo25 = []\n",
    "\n",
    "for epoch in range(num_epochs):\n",
    "    model_SOCAT_Argo25.train()\n",
    "    train_loss_SOCAT_Argo25 = 0.0\n",
    "    for batch_data, batch_target in train_loader_SOCAT_Argo25:\n",
    "        optimizer.zero_grad()\n",
    "        outputs = model_SOCAT_Argo25(batch_data)\n",
    "        loss = criterion(outputs, batch_target.view(-1, 1))\n",
    "        loss.backward()\n",
    "        optimizer.step()        \n",
    "        train_loss_SOCAT_Argo25 += loss.item()\n",
    "            \n",
    "    train_loss_SOCAT_Argo25 /= len(train_loader_SOCAT_Argo25)\n",
    "    train_losses_SOCAT_Argo25.append(train_loss_SOCAT_Argo25)\n",
    "\n",
    "    model_SOCAT_Argo25.eval()\n",
    "    val_loss_SOCAT_Argo25 = 0.0\n",
    "    with torch.no_grad():\n",
    "        for batch_data, batch_target in eval_loader_SOCAT_Argo25:\n",
    "            outputs = model_SOCAT_Argo25(batch_data)\n",
    "            loss = criterion(outputs, batch_target.view(-1, 1))\n",
    "            val_loss_SOCAT_Argo25 += loss.item()\n",
    "\n",
    "    val_loss_SOCAT_Argo25 /= len(eval_loader_SOCAT_Argo25)\n",
    "    val_losses_SOCAT_Argo25.append(val_loss_SOCAT_Argo25)\n",
    "    print(f'Epoch {epoch+1}, Validation Loss: {val_loss_SOCAT_Argo25:.6f}')\n",
    "\n",
    "    early_stopping(val_loss_SOCAT_Argo25, model_SOCAT_Argo25)\n",
    "    if early_stopping.early_stop:\n",
    "        print(\"Early stopping\")\n",
    "        break\n",
    "\n",
    "# Load the best model\n",
    "model_SOCAT_Argo25.load_state_dict(torch.load('checkpoint.pt'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f29fdb5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the model, loss function, and optimizer\n",
    "model_SOCAT_Argo25S = FeedForwardNN()\n",
    "criterion = nn.MSELoss()\n",
    "optimizer = optim.RMSprop(model_SOCAT_Argo25S.parameters(), lr=0.001, alpha=0.9, eps=1e-08, weight_decay=0.0)\n",
    "\n",
    "# Training the model\n",
    "num_epochs = 1600\n",
    "early_stopping = EarlyStopping(patience=30, verbose=True)\n",
    "\n",
    "# Lists to store loss values\n",
    "train_losses_SOCAT_Argo25S = []\n",
    "val_losses_SOCAT_Argo25S = []\n",
    "\n",
    "for epoch in range(num_epochs):\n",
    "    model_SOCAT_Argo25S.train()\n",
    "    train_loss_SOCAT_Argo25S = 0.0\n",
    "    for batch_data, batch_target in train_loader_SOCAT_Argo25S:\n",
    "        optimizer.zero_grad()\n",
    "        outputs = model_SOCAT_Argo25S(batch_data)\n",
    "        loss = criterion(outputs, batch_target.view(-1, 1))\n",
    "        loss.backward()\n",
    "        optimizer.step()        \n",
    "        train_loss_SOCAT_Argo25S += loss.item()\n",
    "            \n",
    "    train_loss_SOCAT_Argo25S /= len(train_loader_SOCAT_Argo25S)\n",
    "    train_losses_SOCAT_Argo25S.append(train_loss_SOCAT_Argo25S)\n",
    "\n",
    "    model_SOCAT_Argo25S.eval()\n",
    "    val_loss_SOCAT_Argo25S = 0.0\n",
    "    with torch.no_grad():\n",
    "        for batch_data, batch_target in eval_loader_SOCAT_Argo25S:\n",
    "            outputs = model_SOCAT_Argo25S(batch_data)\n",
    "            loss = criterion(outputs, batch_target.view(-1, 1))\n",
    "            val_loss_SOCAT_Argo25S += loss.item()\n",
    "\n",
    "    val_loss_SOCAT_Argo25S /= len(eval_loader_SOCAT_Argo25S)\n",
    "    val_losses_SOCAT_Argo25S.append(val_loss_SOCAT_Argo25S)\n",
    "    print(f'Epoch {epoch+1}, Validation Loss: {val_loss_SOCAT_Argo25S:.6f}')\n",
    "\n",
    "    early_stopping(val_loss_SOCAT_Argo25S, model_SOCAT_Argo25S)\n",
    "    if early_stopping.early_stop:\n",
    "        print(\"Early stopping\")\n",
    "        break\n",
    "\n",
    "# Load the best model\n",
    "model_SOCAT_Argo25S.load_state_dict(torch.load('checkpoint.pt'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20d22a4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a figure with three subplots, three figures of loss functions\n",
    "fig, axes = plt.subplots(1, 3, figsize=(18, 5))\n",
    "\n",
    "# Plot 1: SOCAT Data\n",
    "ax = axes[0]\n",
    "ax.plot(train_losses_SOCAT, label='Training Loss')\n",
    "ax.plot(val_losses_SOCAT, label='Evaluation Loss')\n",
    "ax.set_xlabel('Epochs')\n",
    "ax.set_ylabel('Loss')\n",
    "ax.legend()\n",
    "ax.set_title('Training and Evaluation Loss Over Epochs, SOCAT data')\n",
    "\n",
    "# Plot 2: SOCAT + 25% Argo Data\n",
    "ax = axes[1]\n",
    "ax.plot(train_losses_SOCAT_Argo25, label='Training Loss')\n",
    "ax.plot(val_losses_SOCAT_Argo25, label='Evaluation Loss')\n",
    "ax.set_xlabel('Epochs')\n",
    "ax.set_ylabel('Loss')\n",
    "ax.legend()\n",
    "ax.set_title('Training and Evaluation Loss Over Epochs, SOCAT + 25% Argo data')\n",
    "\n",
    "# Plot 3: SOCAT + 25% Argo South Data\n",
    "ax = axes[2]\n",
    "ax.plot(train_losses_SOCAT_Argo25S, label='Training Loss')\n",
    "ax.plot(val_losses_SOCAT_Argo25S, label='Evaluation Loss')\n",
    "ax.set_xlabel('Epochs')\n",
    "ax.set_ylabel('Loss')\n",
    "ax.legend()\n",
    "ax.set_title('Training and Evaluation Loss Over Epochs, SOCAT + 25% Argo South data')\n",
    "\n",
    "# Adjust layout to prevent overlap\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9eb401e0",
   "metadata": {},
   "source": [
    "## c. Validate the model\n",
    "\n",
    "Trained models will be evaluated now on the validation data. Validation data present the field of $pCO_2$ from NEMO/PISCES for the period of February 2008-1010. \n",
    "The reconstructed $pCO_2$ fields are stored in 'recon_pCO2_SOCAT', 'recon_pCO2_SOCAT_Argo25','recon_pCO2_SOCAT_Argo25S'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f329778a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Evaluate the model\n",
    "data_val_SOCAT      = pd.read_csv('Prediction_Variables_normSOCAT_Feb.csv')  # Load validation data here\n",
    "pCO2_list_val_SOCAT = pd.read_csv('Prediction_pCO2_normSOCAT_Feb.csv')  # Load validation targets here\n",
    "data_val_SOCAT      = data_val_SOCAT.drop(columns=['Unnamed: 0'])\n",
    "\n",
    "data_val_tensor_SOCAT      = torch.tensor(data_val_SOCAT.values, dtype=torch.float32)\n",
    "pCO2_list_val_tensor_SOCAT = torch.tensor(pCO2_list_val_SOCAT['pCO2_n'].values, dtype=torch.float32)\n",
    "\n",
    "val_dataset_SOCAT = TensorDataset(data_val_tensor_SOCAT, pCO2_list_val_tensor_SOCAT)\n",
    "val_loader_SOCAT  = DataLoader(val_dataset_SOCAT, batch_size=20, shuffle=False)\n",
    "\n",
    "model_SOCAT.eval()\n",
    "val_loss_SOCAT = 0.0\n",
    "recon_pCO2_SOCAT = []\n",
    "target_pCO2_SOCAT = []\n",
    "with torch.no_grad():\n",
    "    for batch_data, batch_target in val_loader_SOCAT:\n",
    "        outputs = model_SOCAT(batch_data)\n",
    "        recon_pCO2_SOCAT.append(outputs)\n",
    "        target_pCO2_SOCAT.append(batch_target)\n",
    "        loss = criterion(outputs, batch_target.view(-1, 1))\n",
    "        val_loss_SOCAT += loss.item()\n",
    "\n",
    "val_loss_SOCAT /= len(val_loader_SOCAT)\n",
    "print(f'Validation Loss SOCAT: {val_loss_SOCAT:.6f}')\n",
    "\n",
    "\n",
    "# Evaluate the model\n",
    "\n",
    "model_SOCAT_Argo25.eval()\n",
    "val_loss_SOCAT_Argo25 = 0.0\n",
    "recon_pCO2_SOCAT_Argo25 = []\n",
    "target_pCO2_SOCAT_Argo25 = []\n",
    "with torch.no_grad():\n",
    "    for batch_data, batch_target in val_loader_SOCAT:\n",
    "        outputs = model_SOCAT_Argo25(batch_data)\n",
    "        recon_pCO2_SOCAT_Argo25.append(outputs)\n",
    "        target_pCO2_SOCAT_Argo25.append(batch_target)\n",
    "        loss = criterion(outputs, batch_target.view(-1, 1))\n",
    "        val_loss_SOCAT_Argo25 += loss.item()\n",
    "\n",
    "val_loss_SOCAT_Argo25 /= len(val_loader_SOCAT)\n",
    "print(f'Validation Loss SOCAT + 25% Argo: {val_loss_SOCAT_Argo25:.6f}')\n",
    "\n",
    "\n",
    "# Evaluate the model\n",
    "\n",
    "model_SOCAT_Argo25S.eval()\n",
    "val_loss_SOCAT_Argo25S = 0.0\n",
    "recon_pCO2_SOCAT_Argo25S = []\n",
    "target_pCO2_SOCAT_Argo25S = []\n",
    "with torch.no_grad():\n",
    "    for batch_data, batch_target in val_loader_SOCAT:\n",
    "        outputs = model_SOCAT_Argo25S(batch_data)\n",
    "        recon_pCO2_SOCAT_Argo25S.append(outputs)\n",
    "        target_pCO2_SOCAT_Argo25S.append(batch_target)\n",
    "        loss = criterion(outputs, batch_target.view(-1, 1))\n",
    "        val_loss_SOCAT_Argo25S += loss.item()\n",
    "\n",
    "val_loss_SOCAT_Argo25S /= len(val_loader_SOCAT)\n",
    "print(f'Validation Loss SOCAT + 25% Argo South: {val_loss_SOCAT_Argo25S:.6f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a144ae62",
   "metadata": {},
   "source": [
    "## Question:\n",
    "\n",
    "What would you conclude about the accuracy of these three models? "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0ea0ffe",
   "metadata": {},
   "source": [
    "## d. Estimate and visualise model's statistic\n",
    "\n",
    "It is important to plot the reconstructed fields and evaluate the accuracy regionally as the averaged statistics can be good as a esult of error compensation. In this section we will plot the mean field of three reconstructed $pCO_2$ values, their standard deviation (STD), mean differences and correlation coefficients between reconstructed fields and NEMO/PISCES.\n",
    "\n",
    "### Hints: \n",
    "the output of FFNN is normalized $pCO_2$ values, to work with $pCO_2$ in $\\mu$atm we need to convert it back using mean and std values that can be found in the original files. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07bdc807",
   "metadata": {},
   "outputs": [],
   "source": [
    "# concatenate the FFNN output and target values of pCO2\n",
    "result_SOCAT = torch.cat(recon_pCO2_SOCAT, dim=0) \n",
    "result_target_SOCAT = torch.cat(target_pCO2_SOCAT) \n",
    "\n",
    "result_SOCAT_Argo25 = torch.cat(recon_pCO2_SOCAT_Argo25, dim=0) \n",
    "\n",
    "result_SOCAT_Argo25S = torch.cat(recon_pCO2_SOCAT_Argo25S, dim=0) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c551d2ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create a new dataframe with pCO2 values from NEMO/PISCES (pCO2) and three outputs from three ML experiments in uatm\n",
    "pCO2_results = pCO2_list_val_SOCAT[['year','month','day','lat','lon','pCO2']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c22baab",
   "metadata": {},
   "outputs": [],
   "source": [
    "pCO2_results['pCO2_SOCAT'] = (result_SOCAT.numpy() * pCO2_list_train_SOCAT['pCO2_std'][0] + pCO2_list_train_SOCAT['pCO2_mean'][0]).flatten()\n",
    "pCO2_results['pCO2_SOCAT_Argo25'] = (result_SOCAT_Argo25.numpy() * pCO2_list_train_SOCAT_Argo25['pCO2_std'][0] + pCO2_list_train_SOCAT_Argo25['pCO2_mean'][0]).flatten()\n",
    "pCO2_results['pCO2_SOCAT_Argo25S'] = (result_SOCAT_Argo25S.numpy() * pCO2_list_train_SOCAT_Argo25S['pCO2_std'][0] + pCO2_list_train_SOCAT_Argo25S['pCO2_mean'][0]).flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aaa6bc5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#estimate the difference between NEMO/PISCES and FFNN outputs\n",
    "pCO2_results['Diff_SOCAT'] = pCO2_results['pCO2'] - pCO2_results['pCO2_SOCAT']\n",
    "pCO2_results['Diff_SOCAT_Argo25'] = pCO2_results['pCO2'] - pCO2_results['pCO2_SOCAT_Argo25']\n",
    "pCO2_results['Diff_SOCAT_Argo25S'] = pCO2_results['pCO2'] - pCO2_results['pCO2_SOCAT_Argo25S']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9182f0bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "#estimate mean values of dataframe pCO2_results over each grid point\n",
    "pCO2_results_mean = pCO2_results.groupby(['lat','lon']).mean().reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5dc1174",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot mean values of NEMO/PISCES and three FFNN outputs for February 2008-2010\n",
    "fig, axes = plt.subplots(2, 2, figsize=(15, 10), subplot_kw={'projection': ccrs.PlateCarree()})\n",
    "\n",
    "vmin = 600.\n",
    "vmax = 200.\n",
    "\n",
    "# Subplot 1\n",
    "ax = axes[0, 0]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_results_mean['lon'], pCO2_results_mean['lat'], c=pCO2_results_mean['pCO2'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree(), vmin=vmin, vmax=vmax)\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('NEMO/PISCES pCO2, February 2008-2010')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Subplot 2\n",
    "ax = axes[0, 1]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_results_mean['lon'], pCO2_results_mean['lat'], c=pCO2_results_mean['pCO2_SOCAT'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree(), vmin=vmin, vmax=vmax)\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('FFNN pCO2, SOCAT, February 2008-2010')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Subplot 3\n",
    "ax = axes[1, 0]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_results_mean['lon'], pCO2_results_mean['lat'], c=pCO2_results_mean['pCO2_SOCAT_Argo25'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree(), vmin=vmin, vmax=vmax)\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('FFNN pCO2, SOCAT + 25% Argo, February 2008-2010')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Subplot 4\n",
    "ax = axes[1, 1]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_results_mean['lon'], pCO2_results_mean['lat'], c=pCO2_results_mean['pCO2_SOCAT_Argo25S'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree(), vmin=vmin, vmax=vmax)\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('FFNN pCO2, SOCAT + 25% Argo South, February 2008-2010')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Adjust layout to prevent overlap\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f89f8708",
   "metadata": {},
   "source": [
    "## Question:\n",
    "\n",
    "Could you identify regions with the largest differences?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6b04b8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot mean values of differences betwwen NEMO/PISCES and three FFNN outputs for February 2008-2010\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 5), subplot_kw={'projection': ccrs.PlateCarree()})\n",
    "\n",
    "vmin = -30.\n",
    "vmax = 30.\n",
    "\n",
    "# Subplot 1\n",
    "ax = axes[0]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_results_mean['lon'], pCO2_results_mean['lat'], c=pCO2_results_mean['Diff_SOCAT'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree(), vmin=vmin, vmax=vmax)\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('NEMO/PISCES and SOCAT')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Subplot 2\n",
    "ax = axes[1]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_results_mean['lon'], pCO2_results_mean['lat'], c=pCO2_results_mean['Diff_SOCAT_Argo25'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree(), vmin=vmin, vmax=vmax)\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('NEMO/PISCES and SOCAT + 25% Argo')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Subplot 3\n",
    "ax = axes[2]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_results_mean['lon'], pCO2_results_mean['lat'], c=pCO2_results_mean['Diff_SOCAT_Argo25S'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree(), vmin=vmin, vmax=vmax)\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('NEMO/PISCES and SOCAT + 25% Argo South')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Adjust layout to prevent overlap\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00a353d3",
   "metadata": {},
   "source": [
    "## Question:\n",
    "\n",
    "What would you say about the accuracy of each experiments? Do the data from Argo floats improve the reconstruction over the Atlantic Ocean?  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5931cbfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "#estimate std values of dataframe pCO2_results over each grid point\n",
    "pCO2_results_std = pCO2_results.groupby(['lat','lon']).std().reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "370d7bc0",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#plot std values of NEMO/PISCES and three FFNN outputs for February 2008-2010\n",
    "fig, axes = plt.subplots(2, 2, figsize=(15, 10), subplot_kw={'projection': ccrs.PlateCarree()})\n",
    "\n",
    "vmin = 0.\n",
    "vmax = 40.\n",
    "\n",
    "# Subplot 1\n",
    "ax = axes[0, 0]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_results_std['lon'], pCO2_results_std['lat'], c=pCO2_results_std['pCO2'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree(), vmin=vmin, vmax=vmax)\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('NEMO/PISCES pCO2, February 2008-2010')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Subplot 2\n",
    "ax = axes[0, 1]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_results_std['lon'], pCO2_results_std['lat'], c=pCO2_results_std['pCO2_SOCAT'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree(), vmin=vmin, vmax=vmax)\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('Reconstructed pCO2, SOCAT, February 2008-2010')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Subplot 3\n",
    "ax = axes[1, 0]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_results_std['lon'], pCO2_results_std['lat'], c=pCO2_results_std['pCO2_SOCAT_Argo25'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree(), vmin=vmin, vmax=vmax)\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('Reconstructed pCO2, SOCAT + 25% Argo, February 2008-2010')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Subplot 4\n",
    "ax = axes[1, 1]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_results_std['lon'], pCO2_results_std['lat'], c=pCO2_results_std['pCO2_SOCAT_Argo25S'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree(), vmin=vmin, vmax=vmax)\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('Reconstructed pCO2, SOCAT + 25% Argo South, February 2008-2010')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Adjust layout to prevent overlap\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73d47e9f",
   "metadata": {},
   "source": [
    "## Question:\n",
    "\n",
    "Which experiment shows the spatial varibaility comparable to the NEMO/PISCES (first plot)? Using these plots and plots of mean differences what areas would you identify as needing further improvement? Why do these areas show poor statistics?  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e40fd2a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#estimate correlation coefficients betwenn NEMO/PISCES and three FFNN outputs\n",
    "pCO2_results_mean['Corr_SOCAT'] = pCO2_results.groupby(['lat','lon'])[['pCO2','pCO2_SOCAT']].corr().iloc[0::2,-1].reset_index()['pCO2_SOCAT']\n",
    "pCO2_results_mean['Corr_SOCAT_Argo25'] = pCO2_results.groupby(['lat','lon'])[['pCO2','pCO2_SOCAT_Argo25']].corr().iloc[0::2,-1].reset_index()['pCO2_SOCAT_Argo25']\n",
    "pCO2_results_mean['Corr_SOCAT_Argo25S'] = pCO2_results.groupby(['lat','lon'])[['pCO2','pCO2_SOCAT_Argo25S']].corr().iloc[0::2,-1].reset_index()['pCO2_SOCAT_Argo25S']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71da2869",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot correlation coefficients between NEMO/PISCES and three FFNN outputs for February 2008-2010\n",
    "fig, axes = plt.subplots(1, 3, figsize=(15, 5), subplot_kw={'projection': ccrs.PlateCarree()})\n",
    "\n",
    "vmin = -0.2\n",
    "vmax = 1.\n",
    "\n",
    "# Subplot 1\n",
    "ax = axes[0]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_results_mean['lon'], pCO2_results_mean['lat'], c=pCO2_results_mean['Corr_SOCAT'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree(), vmin=vmin, vmax=vmax)\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('NEMO/PISCES and SOCAT')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Subplot 2\n",
    "ax = axes[1]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_results_mean['lon'], pCO2_results_mean['lat'], c=pCO2_results_mean['Corr_SOCAT_Argo25'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree(), vmin=vmin, vmax=vmax)\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('NEMO/PISCES and (SOCAT + 25% Argo)')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Subplot 3\n",
    "ax = axes[2]\n",
    "ax.coastlines()\n",
    "sc = ax.scatter(pCO2_results_mean['lon'], pCO2_results_mean['lat'], c=pCO2_results_mean['Corr_SOCAT_Argo25S'], cmap='viridis', s=5, edgecolor='none', transform=ccrs.PlateCarree(), vmin=vmin, vmax=vmax)\n",
    "plt.colorbar(sc, ax=ax, label='pCO2, uatm')\n",
    "ax.set_title('NEMO/PISCES and (SOCAT + 25% Argo South)')\n",
    "ax.set_xlabel('Longitude')\n",
    "ax.set_ylabel('Latitude')\n",
    "\n",
    "# Adjust layout to prevent overlap\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd6a219c",
   "metadata": {},
   "source": [
    "## Question:\n",
    "\n",
    "How would you interpret the correlation maps regarding our previous results? \n",
    "\n",
    "Which experiment would you advise as an optimal observational network of $pCO_2$ over the Atlantic Ocean? In your answer you should take into account the accuracy of the experiment's results as well as the realism of data distribution and possible cost of the realisation of this distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "daf3fed3",
   "metadata": {},
   "source": [
    "## 5. Conclusion\n",
    "\n",
    "In this practical session you learnt:\n",
    "\n",
    "1. how to open, work and visualise the geographical data in panda dataframe format;\n",
    "\n",
    "2. how to define a simple Feed-Forward Neural Network to reconstruct ocean surface $pCO_2$ using pyTorch library;\n",
    "\n",
    "3. how to use early stopping function in you Neural Network model to prevent overfitting;\n",
    "\n",
    "4. how to visualise model loss function to evaluate the training;\n",
    "\n",
    "5. how to analyse the environmental data using statistics and plotting."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2dadbf7",
   "metadata": {},
   "source": [
    "## 6. References\n",
    "\n",
    "Amari, S., Murata, N., Müller, K.-R., Finke, M., and Yang, H. H.: Asymptotic Statistical Theory of Overtraining and Cross-Validation, IEEE T. Neural Networ., 8, 985–996, 1997. \n",
    "\n",
    "Bakker, D. C. E., Pfeil, B., Landa, C. S., Metzl, N., O'Brien, K. M., Olsen, A., Smith, K., Cosca, C., Harasawa, S., Jones, S. D., Nakaoka, S., Nojiri, Y., Schuster, U., Steinhoff, T., Sweeney, C., Takahashi, T., Tilbrook, B., Wada, C., Wanninkhof, R., Alin, S. R., Balestrini, C. F., Barbero, L., Bates, N. R., Bianchi, A. A., Bonou, F., Boutin, J., Bozec, Y., Burger, E. F., Cai, W.-J., Castle, R. D., Chen, L., Chierici, M., Currie, K., Evans, W., Featherstone, C., Feely, R. A., Fransson, A., Goyet, C., Greenwood, N., Gregor, L., Hankin, S., Hardman-Mountford, N. J., Harlay, J., Hauck, J., Hoppema, M., Humphreys, M. P., Hunt, C. W., Huss, B., Ibánhez, J. S. P., Johannessen, T., Keeling, R., Kitidis, V., Körtzinger, A., Kozyr, A., Krasakopoulou, E., Kuwata, A., Landschützer, P., Lauvset, S. K., Lefèvre, N., Lo Monaco, C., Manke, A., Mathis, J. T., Merlivat, L., Millero, F. J., Monteiro, P. M. S., Munro, D. R., Murata, A., Newberger, T., Omar, A. M., Ono, T., Paterson, K., Pearce, D., Pierrot, D., Robbins, L. L., Saito, S., Salisbury, J., Schlitzer, R., Schneider, B., Schweitzer, R., Sieger, R., Skjelvan, I., Sullivan, K. F., Sutherland, S. C., Sutton, A. J., Tadokoro, K., Telszewski, M., Tuma, M., van Heuven, S. M. A. C., Vandemark, D., Ward, B., Watson, A. J., and Xu, S.: A multi-decade record of high-quality fCO2 data in version 3 of the Surface Ocean CO2 Atlas (SOCAT), Earth Syst. Sci. Data, 8, 383–413, https://doi.org/10.5194/essd-8-383-2016, 2016. \n",
    "\n",
    "Gasparin, F., Guinehut, S., Mao, C., Mirouze, I., Rémy, E., King, R. R., Hamon, M., Reid, R., Storto, A., Le Traon, P. Y., and Martin, M. J.: Requirements for an integrated in situ Atlantic Ocean observing system from coordinated observing system simulation experiments, Front. Mar. Sci., 6, p. 83, https://doi.org/10.3389/fmars.2019.00083, 2019. \n",
    "\n",
    "Kallache, M., Vrac, M., Naveau, P., and Michelangeli, P.-A.: Non-stationary probabilistic downscaling of extreme precipitation, J. Geophys. Res., 116, D05113, https://doi.org/10.1029/2010JD014892, 2011. \n",
    "\n",
    "Landschützer, P., Gruber, N., Bakker, D. C. E., Schuster, U., Nakaoka, S., Payne, M. R., Sasse, T. P., and Zeng, J.: A neural network-based estimate of the seasonal to inter-annual variability of the Atlantic Ocean carbon sink, Biogeosciences, 10, 7793–7815, https://doi.org/10.5194/bg-10-7793-2013, 2013. \n",
    "\n",
    "Takahashi, T., Sutherland, S. C., Wanninkhof, R., Sweeney, C., Feely, R. A., Chipman, D. W., Hales, B., Friederich, G., Chavez, F., Sabine, C., Watson, A., Bakker, D. C. E., Schuster, U., Metzl, N., Yoshikawa-Inoue, H., Ishii, M., Midorikawa, T., Nojiri, Y., Körtzinger, A., Steinhoff, T., Hoppema, M., Olafsson, J., Arnarson, T. S., Tilbrook, B., Johannessen, T., Olsen, A., Bellerby, R., Wong, C. S., Delille, B., Bates, N. R., and de Baar, H. J. W.: Climatological mean and decadal change in surface ocean pCO2, and net sea-air CO2 flux over the global oceans, Deep.-Sea Res. Pt. II, 56, 554–577, https://doi.org/10.1016/j.dsr2.2008.12.009, 2009. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3d82def",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}