notebooks/baseline_oi.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# II- Demo. Optimal Interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An example of simulated SSH data access is provided in the \"example_data_access_meom.ipynb\" notebook. Here, an example of a mapping technique based on a simple optimal interpolation is proposed. The notebook is structured as follow: \n",
    "\n",
    "    1) set optimal interpolation parameters,\n",
    "    2) reading of pseudo-observations,\n",
    "    3) perform optimal interpolation and,\n",
    "    4) save the results (reconstructed SSH field)\n",
    "\n",
    "\n",
    "Here, we assume a vector of observations, noted $y$ defined as:\n",
    "\n",
    "$$y = H x + \\epsilon $$\n",
    "\n",
    "where $H$ is a linear observation operator between the reconstruction grid space and the observation space\n",
    ", $x$ is the state to estimate and $\\epsilon$ is an independent observation error.\n",
    "\n",
    "The optimal interpolation consists in estimating an analysed state $x_{a}$ in combining the available observations to approximate the real state $x$:\n",
    "\n",
    "$$x_{a} = K y $$\n",
    "where $K$ is the weigth matrix defined as:\n",
    "\n",
    "$$ K = BH^T(HBH^T + R)^{-1} $$\n",
    "\n",
    "$B$ is the covariance matrix of $x$, and $R$ the covariance matrix of the error vector $\\epsilon$ ($^T$ is the transpose operator)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy\n",
    "import warnings\n",
    "import logging\n",
    "import sys\n",
    "import os\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logger = logging.getLogger()\n",
    "logger.setLevel(logging.INFO)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.append('..')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.mod_oi import *\n",
    "from src.mod_inout import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1) set optimal interpolation parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# OI Grid\n",
    "lon_min = 295.                                           # domain min longitude\n",
    "lon_max = 305.                                           # domain max longitude\n",
    "lat_min = 33.                                            # domain min latitude\n",
    "lat_max = 43.                                            # domain max latitude\n",
    "time_min = numpy.datetime64('2017-01-01')                # domain min time\n",
    "time_max = numpy.datetime64('2017-12-31')                # domain max time\n",
    "dx = 0.2                                                 # zonal grid spatial step (in degree)\n",
    "dy = 0.2                                                 # meridional grid spatial step (in degree)\n",
    "dt = numpy.timedelta64(1, 'D')                           # temporal grid step\n",
    "\n",
    "glon = numpy.arange(lon_min, lon_max + dx, dx)           # output OI longitude grid\n",
    "glat = numpy.arange(lat_min, lat_max + dy, dy)           # output OI latitude grid\n",
    "gtime = numpy.arange(time_min, time_max + dt, dt)        # output OI time grid\n",
    "\n",
    "# OI parameters\n",
    "Lx = 1.                                                  # Zonal decorrelation scale (in degree)\n",
    "Ly = 1.                                                  # Meridional decorrelation scale (in degree)\n",
    "Lt = 7.                                                  # Temporal decorrelation scale (in days)\n",
    "noise = 0.05                                             # Noise level (5%)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Open your AVISO+ session: fill the ```<AVISO_LOGIN>``` and ```<AVISO_PWD>``` items below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_aviso_session = rq.Session()\n",
    "my_aviso_session.auth = (\"<AVISO_LOGIN>\", \"<AVISO_PWD>\")\n",
    "url_alongtrack = 'https://tds.aviso.altimetry.fr/thredds/dodsC/2021a-SSH-mapping-OSE-along-track-data'\n",
    "url_map = 'https://tds.aviso.altimetry.fr/thredds/dodsC/2021a-SSH-mapping-OSE-grid-data'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2) reading of pseudo-observations + define output folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inputs = [f'{url_alongtrack}/dt_gulfstream_alg_phy_l3_20161201-20180131_285-315_23-53.nc', \n",
    "          f'{url_alongtrack}/dt_gulfstream_j3_phy_l3_20161201-20180131_285-315_23-53.nc', \n",
    "          f'{url_alongtrack}/dt_gulfstream_s3a_phy_l3_20161201-20180131_285-315_23-53.nc',\n",
    "          f'{url_alongtrack}/dt_gulfstream_h2g_phy_l3_20161201-20180131_285-315_23-53.nc',\n",
    "          f'{url_alongtrack}/dt_gulfstream_j2g_phy_l3_20161201-20180131_285-315_23-53.nc',\n",
    "          f'{url_alongtrack}/dt_gulfstream_j2n_phy_l3_20161201-20180131_285-315_23-53.nc'] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define outputs\n",
    "output_directory = '../results/'\n",
    "if not os.path.exists(output_directory):\n",
    "    os.mkdir(output_directory)  \n",
    "output_oi = f'{output_directory}/OSE_ssh_mapping_BASELINE.nc'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3) perform optimal interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# set OI param & grid\n",
    "ds_oi1_param = oi_param(Lx, Ly, Lt, noise)\n",
    "ds_oi1_grid = oi_grid(glon, glat, gtime)\n",
    "# Read input obs + discard a bit...\n",
    "coarsening = {'time': 5}\n",
    "#ds_oi1_obs = read_obs(inputs, ds_oi1_grid, ds_oi1_param, coarsening)\n",
    "ds_oi1_obs = read_obs_from_aviso(inputs, my_aviso_session, ds_oi1_grid, ds_oi1_param, coarsening)\n",
    "# Run OI (take 1h on my laptop)\n",
    "for it in range(len(gtime)):\n",
    "    oi_core(it, ds_oi1_grid, ds_oi1_param, ds_oi1_obs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4) save the results (reconstructed SSH field)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url_ds_mdt = f'{url_map}/mdt.nc'\n",
    "ds_oi1_grid = reformate_oi_output(ds_oi1_grid, url_ds_mdt, my_aviso_session)\n",
    "ds_oi1_grid.to_netcdf(output_oi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "my_conda_env",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}