{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# NCL_panel_13.py\nThis script illustrates the following concepts:\n   - Overlaying a vector field over filled contours\n   - Paneling two plots vertically\n\nSee following URLs to see the reproduced NCL plot & script:\n    - Original NCL script: http://www.ncl.ucar.edu/Applications/Scripts/panel_13.ncl\n    - Original NCL plot: http://www.ncl.ucar.edu/Applications/Images/panel_13_lg.png\n    \nNote:\n    Due to differences in how NCL and Python scale glyphs in vector fields, the\n    smallest vectors in the Python version are much harder to read than in the\n    NCL version. An issue has been opened on the geoCAT examples gallery github\n    so this can be addressed at a later date.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Import packages:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import cartopy.crs as ccrs\nfrom cartopy.mpl.gridliner import LongitudeFormatter, LatitudeFormatter\nimport matplotlib.pyplot as plt\nfrom matplotlib.ticker import FormatStrFormatter\nimport matplotlib.patches as mpatches\nimport numpy as np\nimport xarray as xr\n\nimport geocat.datafiles as gdf\nfrom geocat.viz import cmaps as gvcmaps\nimport geocat.viz.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\"))\n\n# Extract data from second timestep\nds = ds.isel(time=1).drop_vars('time')\n\n# Ensure longitudes range from 0 to 360 degrees\nU = gvutil.xr_add_cyclic_longitudes(ds.U, \"lon\")\nV = gvutil.xr_add_cyclic_longitudes(ds.V, \"lon\")\n\n# Thin data to only include every fourth values\nU = U[::4, ::4]\nV = V[::4, ::4]\n\n# Calculate the magnitude of the winds\nmagnitude = np.sqrt(U.data**2 + V.data**2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Create sublots and specify their projections\nprojection = ccrs.PlateCarree()\nfig, axs = plt.subplots(2,\n                        1,\n                        figsize=(7, 10),\n                        subplot_kw={\"projection\": projection})\nplt.tight_layout(pad=4, h_pad=-5)\n\n# Add coastlines, the zorder keyword specifies the order in which the elements\n# are drawn where elements with lower zorder values are drawn first\naxs[0].coastlines(linewidth=0.5, zorder=1)\naxs[1].coastlines(linewidth=0.5, zorder=1)\n\n# Use geocat.viz.util convenience function to set axes tick values\ngvutil.set_axes_limits_and_ticks(axs[0],\n                                 xlim=[-180, 180],\n                                 ylim=[-90, 90],\n                                 xticks=np.arange(-180, 181, 30),\n                                 yticks=np.arange(-90, 91, 30))\ngvutil.set_axes_limits_and_ticks(axs[1],\n                                 xlim=[-180, 180],\n                                 ylim=[-90, 90],\n                                 xticks=np.arange(-180, 181, 30),\n                                 yticks=np.arange(-90, 91, 30))\n\n# Use geocat.viz.util convenience function to add minor and major tick lines\ngvutil.add_major_minor_ticks(axs[0])\ngvutil.add_major_minor_ticks(axs[1])\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(axs[0])\ngvutil.add_lat_lon_ticklabels(axs[1])\n# Remove the degree symbol from tick labels\naxs[0].yaxis.set_major_formatter(LatitudeFormatter(degree_symbol=''))\naxs[0].xaxis.set_major_formatter(LongitudeFormatter(degree_symbol=''))\naxs[1].yaxis.set_major_formatter(LatitudeFormatter(degree_symbol=''))\naxs[1].xaxis.set_major_formatter(LongitudeFormatter(degree_symbol=''))\n\n# Use geocat.viz.util convenience function to set titles and labels\ngvutil.set_titles_and_labels(axs[0],\n                             lefttitle='Speed',\n                             lefttitlefontsize=10,\n                             righttitle=U.units,\n                             righttitlefontsize=10)\ngvutil.set_titles_and_labels(axs[1],\n                             lefttitle='Wind',\n                             lefttitlefontsize=10,\n                             righttitle=U.units,\n                             righttitlefontsize=10)\n\n# Load in colormap\nnewcmap = gvcmaps.gui_default\n\n# Specify contour levels and contour ticks\nspeed_levels = np.arange(0, 40, 2.5)\nspeed_ticks = np.arange(2.5, 37.5, 2.5)\nwind_levels = np.arange(-16, 44, 4)\nwind_ticks = np.arange(-12, 40, 4)\n\n# Plot filled contours\nspeed = axs[0].contourf(U['lon'],\n                        U['lat'],\n                        magnitude,\n                        levels=speed_levels,\n                        cmap=newcmap,\n                        zorder=0)\nwind = axs[1].contourf(U['lon'],\n                       U['lat'],\n                       U.data,\n                       levels=wind_levels,\n                       cmap=newcmap,\n                       zorder=0)\n\n# Create color bars\nspeed_cbar = plt.colorbar(speed,\n                          ax=axs[0],\n                          orientation='horizontal',\n                          ticks=speed_ticks,\n                          shrink=0.8,\n                          drawedges=True,\n                          pad=0.1)\nplt.colorbar(wind,\n             ax=axs[1],\n             orientation='horizontal',\n             ticks=wind_ticks,\n             shrink=0.8,\n             drawedges=True,\n             pad=0.1)\n# Remove trailing zeros from speed color bar tick labels\nspeed_cbar.ax.xaxis.set_major_formatter(FormatStrFormatter('%g'))\n\n# Plotting vector field\nquiver_speed = axs[0].quiver(U['lon'],\n                             U['lat'],\n                             U.data,\n                             V.data,\n                             scale=400,\n                             width=0.002,\n                             headwidth=6,\n                             headlength=7,\n                             zorder=2)\nquiver_wind = axs[1].quiver(U['lon'],\n                            U['lat'],\n                            U.data,\n                            V.data,\n                            scale=400,\n                            width=0.002,\n                            headwidth=6,\n                            headlength=7,\n                            zorder=2)\n\n# Add white box to go behind reference vector\naxs[0].add_patch(\n    mpatches.Rectangle(xy=[0.775, 0],\n                       width=0.225,\n                       height=0.2,\n                       facecolor='white',\n                       transform=axs[0].transAxes,\n                       zorder=2))\naxs[1].add_patch(\n    mpatches.Rectangle(xy=[0.775, 0],\n                       width=0.225,\n                       height=0.2,\n                       facecolor='white',\n                       transform=axs[1].transAxes,\n                       zorder=2))\n# Add reference vector and label\naxs[0].quiverkey(quiver_speed, 0.8875, 0.1, 20, 20, zorder=2)\naxs[1].quiverkey(quiver_wind, 0.8875, 0.1, 20, 20, zorder=2)\naxs[0].text(0.785,\n            0.025,\n            \"Reference Vector\",\n            transform=axs[0].transAxes,\n            zorder=2)\naxs[1].text(0.785,\n            0.025,\n            \"Reference Vector\",\n            transform=axs[1].transAxes,\n            zorder=2)\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
}