Notebooks/run_diva2D_MLD.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "This notebook presents an application of the [pydiva2D](../pydiva/pydiva2d.py) module to perform a 2D Diva analysis and then generate figures for the different input and output using the `Basemap` module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import logging\n",
    "from pydiva import pydiva2d\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rcParams\n",
    "from mpl_toolkits.basemap import Basemap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We modify the Matplotlib parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rcParams['figure.dpi'] = 300\n",
    "rcParams['font.size'] = 14.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppress mpl_toolkits warnings due to version mismatch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "import matplotlib.cbook\n",
    "warnings.filterwarnings(\"ignore\", category=matplotlib.cbook.mplDeprecation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If needed, modify the logger level to control the messages displayed on screen (by default: warnings)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "logger = logging.getLogger('pydiva2d')\n",
    "logger.setLevel(logging.WARN)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prepare input files and directories\n",
    "\n",
    "In this example the input files are already created and can be found in the *data* directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "datadir = '../data/'\n",
    "datafile = os.path.join(datadir, 'MLD1.dat')\n",
    "coastfile = os.path.join(datadir, 'coast.cont')\n",
    "paramfile = os.path.join(datadir, 'param.par')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also create a directory *figures* where we will save the plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "figdir = './figures/'\n",
    "if os.path.exists(figdir):\n",
    "    logger.debug(\"Directory {0} already exists\".format(figdir))\n",
    "else:\n",
    "    logger.debug(\"Creating directory {0}\".format(figdir))\n",
    "    os.makedirs(figdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create paths and file names\n",
    "\n",
    "We create the paths for the Diva directories and files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "divadir = \"/home/ctroupin/Software/DIVA/DIVA-diva-4.7.1\"\n",
    "DivaDirs = pydiva2d.DivaDirectories(divadir)\n",
    "DivaFiles = pydiva2d.Diva2Dfiles(DivaDirs.diva2d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read input files\n",
    "\n",
    "Let's read the information from the input files.<br>\n",
    "For each input file, we first creat an object (i.e., **Data**, **Contour** and **Param**) that will allow us to perform operations such as \n",
    "* reading, \n",
    "* writing, \n",
    "* plotting, \n",
    "* getting a description."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Data = pydiva2d.Diva2DData().read_from(datafile)\n",
    "Contour = pydiva2d.Diva2DContours().read_from(coastfile)\n",
    "Param = pydiva2d.Diva2DParameters().read_from(paramfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example with the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Param.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that `ispec=0` (no error field calculated), but we can easily change that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Param.ispec = -101\n",
    "Param.varbak = 1.0\n",
    "Param.write_to(paramfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Make the analysis\n",
    "We simply call the `make` method of the `Diva2DResults` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "analysis2d = pydiva2d.Diva2DResults().make(divadir, datafile=datafile,\n",
    "                                           contourfile=coastfile,\n",
    "                                           paramfile=paramfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mesh is created if necessary, so we can load it afterwards:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh2d = pydiva2d.Diva2DMesh().read_from(DivaFiles.mesh, DivaFiles.meshtopo)\n",
    "mesh2d.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create the plots\n",
    "\n",
    "In this example the plots are created using the Basemap module.<br>\n",
    "All the objects we created have a function *add_to_plot* that takes into account, if it is defined, the selected projection in Basemap.\n",
    "\n",
    "## Initialise the projection\n",
    "\n",
    "We will center the map on the region specified by the parameters and with a Mercator projection.<br>\n",
    "Type `Basemap?` for a list of available projections and options."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "m = Basemap(projection='merc', llcrnrlon=Param.xori, llcrnrlat=Param.yori, \n",
    "            urcrnrlon=Param.xend, urcrnrlat=Param.yend, \n",
    "            resolution='i')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data points\n",
    "\n",
    "We use the method **add_to_plot** with 'm' (the projection) as an optional argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.subplot(111)\n",
    "m.ax = ax\n",
    "dataplot = Data.add_to_plot(m=m, s=10)\n",
    "m.drawmeridians(np.linspace(Param.xori, Param.xend, 4), labels=[0, 0, 0, 1], linewidth=0.2)\n",
    "m.drawparallels(np.linspace(Param.yori, Param.yend, 4), labels=[1, 0, 0, 0], linewidth=0.2)\n",
    "m.drawcoastlines(linewidth=0.2)\n",
    "plt.colorbar(dataplot)\n",
    "plt.title('Data points')\n",
    "# plt.savefig(os.path.join(figdir, 'datapoints.png'))\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.subplot(111)\n",
    "m.ax = ax\n",
    "Contour.add_to_plot(m=m, linewidth=3)\n",
    "m.drawmeridians(np.linspace(Param.xori, Param.xend, 4), labels=[0, 0, 0, 1], linewidth=0.2)\n",
    "m.drawparallels(np.linspace(Param.yori, Param.yend, 4), labels=[1, 0, 0, 0], linewidth=0.2)\n",
    "m.drawcoastlines(linewidth=0.2, zorder=3)\n",
    "plt.title('Contours')\n",
    "# plt.savefig(os.path.join(figdir, 'contours.png'))\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finite element mesh\n",
    "\n",
    "To create the plot, we add the option *m* in the `add_to_plot` method.     \n",
    "Take into account that this plot generation can be relatively slow for a dense mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.subplot(111)\n",
    "m.ax = ax\n",
    "mesh2d.add_to_plot(m, linewidth=0.25, color='k')\n",
    "m.drawmeridians(np.linspace(Param.xori, Param.xend, 4), labels=[0, 0, 0, 1], linewidth=0.2)\n",
    "m.drawparallels(np.linspace(Param.yori, Param.yend, 4), labels=[1, 0, 0, 0], linewidth=0.2)\n",
    "m.drawcoastlines(linewidth=0.2, zorder=3)\n",
    "plt.title(\"Finite-element mesh\")\n",
    "# plt.savefig(os.path.join(figdir, 'mesh.png'))\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysed field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.subplot(111)\n",
    "m.ax = ax\n",
    "FieldPlot = analysis2d.add_to_plot(field='analysis', m=m)\n",
    "m.drawmeridians(np.linspace(Param.xori, Param.xend, 4), labels=[0, 0, 0, 1], linewidth=0.2)\n",
    "m.drawparallels(np.linspace(Param.yori, Param.yend, 4), labels=[1, 0, 0, 0], linewidth=0.2)\n",
    "m.drawcoastlines(linewidth=0.2, zorder=3)\n",
    "plt.colorbar(FieldPlot)\n",
    "plt.title('Analyzed field')\n",
    "# plt.savefig(os.path.join(figdir, 'analysis.png'))\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Error field\n",
    "The `ispec=-101` parameter means we want to obtain the error field in the same grid as the analysis using the so-called *almost exact error calculation*.<br>\n",
    "We can display the whole error field:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.subplot(111)\n",
    "m.ax = ax\n",
    "errorplot = analysis2d.add_to_plot(field='error', m=m, cmap=plt.cm.hot_r)\n",
    "m.drawmeridians(np.linspace(Param.xori, Param.xend, 4), labels=[0, 0, 0, 1], linewidth=0.2)\n",
    "m.drawparallels(np.linspace(Param.yori, Param.yend, 4), labels=[1, 0, 0, 0], linewidth=0.2)\n",
    "m.drawcoastlines(linewidth=0.2, zorder=3)\n",
    "plt.colorbar(errorplot)\n",
    "plt.title('Error field')\n",
    "# plt.savefig(os.path.join(figdir, 'errorfield.png'))\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or use it to mask the analysis for example when the error is above 20%:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "maskedfield = np.ma.masked_where(analysis2d.error>=0.5, analysis2d.analysis)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = plt.subplot(111)\n",
    "m.ax = ax\n",
    "xx, yy = np.meshgrid(analysis2d.x, analysis2d.y)\n",
    "maskedplot = m.pcolormesh(xx, yy, maskedfield, latlon=True)\n",
    "m.drawmeridians(np.linspace(Param.xori, Param.xend, 4), labels=[0, 0, 0, 1], linewidth=0.2)\n",
    "m.drawparallels(np.linspace(Param.yori, Param.yend, 4), labels=[1, 0, 0, 0], linewidth=0.2)\n",
    "m.drawcoastlines(linewidth=0.2, zorder=3)\n",
    "plt.colorbar(maskedplot)\n",
    "plt.title('Masked analysis field')\n",
    "# plt.savefig(os.path.join(figdir, 'maskedanalysis.png'))\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Final remarks\n",
    "\n",
    "* The different elements of the plot can be combined in many ways (see example below).\n",
    "* All the plots can be easily modified by adapting the diffent options: colors, projection, font size etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.subplot(111)\n",
    "m.ax = ax\n",
    "fieldplot = analysis2d.add_to_plot(m=m, vmin=-60, vmax=-5)\n",
    "mesh2d.add_to_plot(m=m, alpha=0.8, linewidth=0.25, color='k', zorder=4)\n",
    "Contour.add_to_plot(m=m, linewidth=3, color='0.5')\n",
    "Data.add_to_plot(m=m, s=3, vmin=-60, vmax=-5, zorder=5)\n",
    "m.drawmeridians(np.linspace(Param.xori, Param.xend, 4), labels=[0, 0, 0, 1], linewidth=0.2)\n",
    "m.drawparallels(np.linspace(Param.yori, Param.yend, 4), labels=[1, 0, 0, 0], linewidth=0.26)\n",
    "m.drawcoastlines(linewidth=0.2, zorder=3)\n",
    "plt.colorbar(fieldplot)\n",
    "plt.title('Analyzed field, mesh,\\n contour and data')\n",
    "#plt.savefig(os.path.join(figdir, 'AnalysisMeshData.png'))\n",
    "plt.show()\n",
    "plt.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Diva-python3.6",
   "language": "python",
   "name": "diva-python3.6"
  },
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}