{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# NCL_conOncon_2.py\nThis script illustrates the following concepts:\n   - Overlaying two sets of contours on a map\n   - Drawing the zero contour line thicker\n   - Changing the center longitude for a cylindrical equidistant projection\n   - Using a blue-white-red color map\n\nSee following URLs to see the reproduced NCL plot & script:\n    - Original NCL script: https://www.ncl.ucar.edu/Applications/Scripts/conOncon_2.ncl\n    - Original NCL plot: https://www.ncl.ucar.edu/Applications/Images/conOncon_2_lg.png\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Import packages:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport xarray as xr\nimport matplotlib.pyplot as plt\nimport cartopy.crs as ccrs\nimport cartopy.feature as cfeature\nfrom cartopy.mpl.gridliner import LongitudeFormatter, LatitudeFormatter\n\nimport geocat.datafiles as gdf\nfrom geocat.viz import util as gvutil\nfrom geocat.viz import cmaps as gvcmaps"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Read in data:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Open a netCDF data file using xarray default engine and load the data into xarrays\nsst = xr.open_dataset(gdf.get(\"netcdf_files/sst8292a.nc\"))\nolr = xr.open_dataset(gdf.get(\"netcdf_files/olr7991a.nc\"))\n\n# Extract data for December 1982\nsst = sst.isel(time=11, drop=True).SSTA\nolr = olr.isel(time=47, drop=True).OLRA\n\n# Fix the artifact of not-shown-data around 0 and 360-degree longitudes\nsst = gvutil.xr_add_cyclic_longitudes(sst, 'lon')\nolr = gvutil.xr_add_cyclic_longitudes(olr, 'lon')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Generate figure and axes\nplt.figure(figsize=(8, 8))\n\n# Set axes projection\nax = plt.axes(projection=ccrs.PlateCarree(central_longitude=-160))\nax.set_extent([100, 300, -60, 60], crs=ccrs.PlateCarree())\n\n# Load in color map and specify contour levels\ncmap = gvcmaps.BlWhRe\nsst_levels = np.arange(-5.5, 6, 0.5)\n# Draw SST contour\ntemp = sst.plot.contourf(ax=ax,\n                         transform=ccrs.PlateCarree(),\n                         cmap=cmap,\n                         levels=sst_levels,\n                         extend='neither',\n                         add_colorbar=False,\n                         add_labels=False,\n                         zorder=0)\nplt.colorbar(temp,\n             ax=ax,\n             orientation='vertical',\n             ticks=np.arange(-5, 6, 1),\n             drawedges=True,\n             shrink=0.5,\n             aspect=10)\n\n# Draw map features on top of filled contour\nax.add_feature(cfeature.LAND, facecolor='lightgray', zorder=1)\nax.add_feature(cfeature.COASTLINE, edgecolor='gray', linewidth=0.5, zorder=1)\n\n# Draw OLR contour\n# Specify contour levels excluding 0\nolr_levels = np.arange(-80, 0, 10)\nolr_levels = np.append(olr_levels, np.arange(10, 50, 10))\n\nrad = olr.plot.contour(ax=ax,\n                       transform=ccrs.PlateCarree(),\n                       levels=olr_levels,\n                       colors='gray',\n                       linewidths=0.5,\n                       add_labels=False)\nax.clabel(rad, [-40, -20, 20], fmt='%d', inline=True, colors='black')\n\n# Plot the zero contour with a black color\nrad = olr.plot.contour(ax=ax,\n                       transform=ccrs.PlateCarree(),\n                       levels=[0],\n                       colors='black',\n                       linewidths=0.5,\n                       add_labels=False)\nax.clabel(rad, [0], fmt='%d', inline=True, colors='black')\n\n# Use geocat.viz.util convenience function to set axes tick values\ngvutil.set_axes_limits_and_ticks(ax,\n                                 ylim=(-60, 60),\n                                 yticks=np.arange(-60, 90, 30),\n                                 xticks=np.arange(-80, 120, 30))\n\n# Use geocat.viz.util convenience function to make plots look like NCL plots by\n# using latitude, longitude tick labels\ngvutil.add_lat_lon_ticklabels(ax)\n\n# Remove the degree symbol from tick labels\nax.yaxis.set_major_formatter(LatitudeFormatter(degree_symbol=''))\nax.xaxis.set_major_formatter(LongitudeFormatter(degree_symbol=''))\n\n# Use geocat.viz.util convenience function to add minor and major tick lines\ngvutil.add_major_minor_ticks(ax,\n                             x_minor_per_major=3,\n                             y_minor_per_major=3,\n                             labelsize=10)\n\ngvutil.set_titles_and_labels(ax,\n                             maintitle=olr.long_name,\n                             maintitlefontsize=14,\n                             lefttitle='degC',\n                             lefttitlefontsize=12,\n                             righttitle='(W m s$^{-2}$)',\n                             righttitlefontsize=12)\n# Add center title\nax.text(0.35, 1.06, 'December 1982', fontsize=12, transform=ax.transAxes)\n\n# Add lower text box\nax.text(1,\n        -0.2,\n        \"CONTOUR FROM -80 TO 40 BY 10\",\n        horizontalalignment='right',\n        transform=ax.transAxes,\n        bbox=dict(boxstyle='square, pad=0.25',\n                  facecolor='white',\n                  edgecolor='black'))\n\nplt.show()"
      ]
    }
  ],
  "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.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}