{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# NCL_polar_1.py\nThis script illustrates the following concepts:\n    - Drawing black-and-white contours over a polar stereographic map\n    - Drawing the northern hemisphere of a polar stereographic map\n\nSee following URLs to see the reproduced NCL plot & script:\n    - Original NCL script: https://www.ncl.ucar.edu/Applications/Scripts/polar_1_lg.ncl\n    - Original NCL plot: https://www.ncl.ucar.edu/Applications/Images/polar_1_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 cartopy.feature as cfeature\nimport cartopy.crs as ccrs\nimport matplotlib.pyplot as plt\nimport matplotlib.path as mpath\nimport matplotlib.ticker as mticker\n\nimport geocat.datafiles as gdf\nfrom geocat.viz import util as gvutil"
      ]
    },
    {
      "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\nds = xr.open_dataset(gdf.get(\"netcdf_files/uv300.nc\"))\nU = ds.U[1, :, :]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Fix the artifact of not-shown-data around 0 and 360-degree longitudes\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "wrap_U = gvutil.xr_add_cyclic_longitudes(U, \"lon\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Generate axes, using Cartopy, drawing coastlines, and adding features\nfig = plt.figure(figsize=(10, 10))\nprojection = ccrs.NorthPolarStereo()\nax = plt.axes(projection=projection)\nax.add_feature(cfeature.LAND, facecolor='lightgray')\n\n# Set map boundary to include latitudes between 0 and 40 and longitudes\n# between -180 and 180 only\ngvutil.set_map_boundary(ax, [-180, 180], [0, 40], south_pad=1)\n\n# Set draw_labels to False so that you can manually manipulate it later\ngl = ax.gridlines(ccrs.PlateCarree(),\n                  draw_labels=False,\n                  linestyle=\"--\",\n                  color='black')\n\n# Manipulate latitude and longitude gridline numbers and spacing\ngl.ylocator = mticker.FixedLocator(np.arange(0, 90, 15))\ngl.xlocator = mticker.FixedLocator(np.arange(-180, 180, 30))\n\n# Manipulate longitude labels (0, 30 E, 60 E, ..., 30 W, etc.)\nticks = np.arange(0, 210, 30)\netick = ['0'] + [\n    r'%dE' % tick for tick in ticks if (tick != 0) & (tick != 180)\n] + ['180']\nwtick = [r'%dW' % tick for tick in ticks if (tick != 0) & (tick != 180)]\nlabels = etick + wtick\nxticks = np.arange(0, 360, 30)\nyticks = np.full_like(xticks, -5)  # Latitude where the labels will be drawn\n\nfor xtick, ytick, label in zip(xticks, yticks, labels):\n    if label == '180':\n        ax.text(xtick,\n                ytick,\n                label,\n                fontsize=14,\n                horizontalalignment='center',\n                verticalalignment='top',\n                transform=ccrs.Geodetic())\n    elif label == '0':\n        ax.text(xtick,\n                ytick,\n                label,\n                fontsize=14,\n                horizontalalignment='center',\n                verticalalignment='bottom',\n                transform=ccrs.Geodetic())\n    else:\n        ax.text(xtick,\n                ytick,\n                label,\n                fontsize=14,\n                horizontalalignment='center',\n                verticalalignment='center',\n                transform=ccrs.Geodetic())\n\n# Contour-plot U-data\np = wrap_U.plot.contour(ax=ax,\n                        vmin=-8,\n                        vmax=16,\n                        transform=ccrs.PlateCarree(),\n                        levels=np.arange(-12, 44, 4),\n                        linewidths=0.5,\n                        cmap='black',\n                        add_labels=False)\n\nax.clabel(p, np.arange(-8, 17, 8), fmt='%d', inline=1, fontsize=14)\n\n# Use geocat.viz.util convenience function to add left and right titles\ngvutil.set_titles_and_labels(ax, lefttitle=\"Zonal Wind\", righttitle=\"m/s\")\n\n# Add lower text box\nax.text(1.0,\n        -.10,\n        \"CONTOUR FROM -12 TO 40 BY 4\",\n        horizontalalignment='right',\n        transform=ax.transAxes,\n        bbox=dict(boxstyle='square, pad=0.25',\n                  facecolor='white',\n                  edgecolor='black'))\n\n# Show the plot\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
}