{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import (absolute_import, division, print_function, unicode_literals)\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.cm import get_cmap\n",
    "import cartopy.crs as crs\n",
    "import cartopy.feature as cfeature\n",
    "from netCDF4 import Dataset as nc\n",
    "\n",
    "from wrf import npvalues, getvar, smooth2d\n",
    "\n",
    "# Open the output NetCDF with netcdf-python\n",
    "filename = \"/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00\"\n",
    "ncfile = nc(filename)\n",
    "\n",
    "# Get sea level pressure and cloud top temperature\n",
    "slp = getvar(ncfile, \"slp\")\n",
    "ctt = getvar(ncfile, \"ctt\")\n",
    "\n",
    "# Smooth the SLP\n",
    "smooth_slp = smooth2d(slp, 3)\n",
    "\n",
    "# Extract the latitude and longitude coordinate arrays as regular numpy array instead of xarray.DataArray\n",
    "lons = npvalues(slp.coords[\"XLONG\"])\n",
    "lats = npvalues(slp.coords[\"XLAT\"])\n",
    "\n",
    "# Get the cartopy projection class\n",
    "wrf_proj = slp.attrs[\"projection\"]\n",
    "cart_proj = wrf_proj.cartopy()\n",
    "\n",
    "# Create the figure\n",
    "fig = plt.figure(figsize=(4,4))\n",
    "ax = plt.axes([0.1,0.1,0.8,0.8], projection=cart_proj)\n",
    "\n",
    "# Download and create the states, land, and oceans using cartopy features\n",
    "states = cfeature.NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n",
    "                             name='admin_1_states_provinces_shp')\n",
    "land = cfeature.NaturalEarthFeature(category='physical', name='land', scale='50m', \n",
    "                                    facecolor=cfeature.COLORS['land'])\n",
    "ocean = cfeature.NaturalEarthFeature(category='physical', name='ocean', scale='50m', \n",
    "                                     facecolor=cfeature.COLORS['water'])\n",
    "\n",
    "# Make the pressure contours.\n",
    "contour_levels = [960, 965, 970, 975, 980, 990]\n",
    "c1 = plt.contour(lons, lats, npvalues(smooth_slp), levels=contour_levels, colors=\"white\", \n",
    "            transform=crs.PlateCarree(), zorder=3, linewidths=1.0)\n",
    "\n",
    "# Add pressure contour labels\n",
    "#plt.clabel(c1, contour_levels, inline=True, fmt='%.0f', fontsize=7)\n",
    "\n",
    "# Create the filled cloud top temperature contours\n",
    "contour_levels = [-80, -70, -60, -50, -40, -30, -20, -10, 0, 10]\n",
    "plt.contourf(lons, lats, npvalues(ctt), contour_levels, cmap=get_cmap(\"Greys\"),\n",
    "             transform=crs.PlateCarree(), zorder=2)\n",
    "\n",
    "plt.plot([-80,-77.8], [26.75,26.75], color=\"yellow\", marker=\"o\", transform=crs.PlateCarree(), zorder=3)\n",
    "\n",
    "# Create the label bar for cloud top temperature\n",
    "#cb2 = plt.colorbar(ax=ax, fraction=0.046, pad=0.04)\n",
    "cb2 = plt.colorbar(ax=ax)\n",
    "\n",
    "# Draw the oceans, land, and states\n",
    "ax.add_feature(ocean)\n",
    "ax.add_feature(land)\n",
    "ax.add_feature(states, linewidth=.5, edgecolor=\"black\")\n",
    "\n",
    "# Crop the domain to the region around the hurricane\n",
    "ax.set_extent([-85., -75.0, np.amin(lats), 30.0], crs=crs.PlateCarree())\n",
    "ax.gridlines(crs=crs.PlateCarree(), draw_labels=False)\n",
    "\n",
    "# Add the title and show the image\n",
    "#plt.title(\"Hurricane Matthew Cloud Top Temperature (degC) \")\n",
    "plt.savefig(\"/Users/ladwig/Documents/workspace/wrf_python/doc/source/_static/images/matthew.png\",\n",
    "           transparent=True, bbox_inches=\"tight\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from __future__ import (absolute_import, division, print_function, unicode_literals)\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.cm import get_cmap\n",
    "import cartopy.crs as crs\n",
    "from cartopy.feature import NaturalEarthFeature\n",
    "from Nio import open_file\n",
    "\n",
    "from wrf import npvalues, getvar, smooth2d, ll_to_xy, CoordPair, vertcross, getproj, get_proj_params, to_xy_coords\n",
    "\n",
    "# Open the output NetCDF file with PyNIO\n",
    "filename = \"/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00\"\n",
    "pynio_filename = filename + \".nc\"\n",
    "ncfile = open_file(pynio_filename)\n",
    "\n",
    "# Extract pressure and model height\n",
    "z = getvar(ncfile, \"z\", timeidx=0)\n",
    "dbz = getvar(ncfile, \"dbz\", timeidx=0)\n",
    "\n",
    "wspd =  getvar(ncfile, \"uvmet_wspd_wdir\", units=\"kt\")[0,:]\n",
    "Z = 10**(dbz/10.)\n",
    "\n",
    "start_point = CoordPair(lat=26.75, lon=-80.0)\n",
    "end_point = CoordPair(lat=26.75, lon=-77.8)\n",
    "\n",
    "# Compute the vertical cross-section interpolation.  Also, include the lat/lon points along the cross-section.\n",
    "z_cross = vertcross(Z, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True)\n",
    "wspd_cross = vertcross(wspd, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True)\n",
    "dbz_cross = 10.0 * np.log10(z_cross)\n",
    "\n",
    "# Create the figure\n",
    "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(4,4))\n",
    "#ax = plt.axes([0.1,0.1,0.8,0.8])\n",
    "\n",
    "# Define the contour levels [0, 50, 100, 150, ....]\n",
    "levels = [5 + 5*n for n in range(15)]\n",
    "\n",
    "# Make the contour plot\n",
    "a = axes[0].contourf(npvalues(wspd_cross))\n",
    "# Add the color bar\n",
    "fig.colorbar(a, ax=axes[0])\n",
    "\n",
    "b = axes[1].contourf(npvalues(dbz_cross), levels=levels)\n",
    "fig.colorbar(b, ax=axes[1])\n",
    "\n",
    "# Set the x-ticks to use latitude and longitude labels.\n",
    "coord_pairs = npvalues(dbz_cross.coords[\"xy_loc\"])\n",
    "x_ticks = np.arange(coord_pairs.shape[0])\n",
    "x_labels = [pair.latlon_str() for pair in npvalues(coord_pairs)]\n",
    "axes[0].set_xticks(x_ticks[::20])\n",
    "axes[0].set_xticklabels([], rotation=45)\n",
    "axes[1].set_xticks(x_ticks[::20])\n",
    "axes[1].set_xticklabels(x_labels[::20], rotation=45, fontsize=6) \n",
    "\n",
    "\n",
    "# Set the y-ticks to be height.\n",
    "vert_vals = npvalues(dbz_cross.coords[\"vertical\"])\n",
    "v_ticks = np.arange(vert_vals.shape[0])\n",
    "axes[0].set_yticks(v_ticks[::20])\n",
    "axes[0].set_yticklabels(vert_vals[::20], fontsize=6) \n",
    "axes[1].set_yticks(v_ticks[::20])\n",
    "axes[1].set_yticklabels(vert_vals[::20], fontsize=6) \n",
    "\n",
    "# Set the x-axis and  y-axis labels\n",
    "axes[1].set_xlabel(\"Latitude, Longitude\", fontsize=7)\n",
    "axes[0].set_ylabel(\"Height (m)\", fontsize=7)\n",
    "axes[1].set_ylabel(\"Height (m)\", fontsize=7)\n",
    "\n",
    "# Add a title\n",
    "axes[0].set_title(\"Cross-Section of Wind Speed (kt)\", {\"fontsize\" : 10})\n",
    "axes[1].set_title(\"Cross-Section of Reflectivity (dBZ)\", {\"fontsize\" : 10})\n",
    "\n",
    "plt.savefig(\"/Users/ladwig/Documents/workspace/wrf_python/doc/source/_static/images/matthew_cross.png\",\n",
    "           transparent=True, bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from __future__ import (absolute_import, division, print_function, unicode_literals)\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.cm import get_cmap\n",
    "import cartopy.crs as crs\n",
    "import cartopy.feature as cfeature\n",
    "from Nio import open_file\n",
    "\n",
    "from wrf import getvar, npvalues, vertcross, smooth2d, CoordPair\n",
    "\n",
    "# Open the output NetCDF file with PyNIO\n",
    "filename = \"/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00\"\n",
    "pynio_filename = filename + \".nc\"\n",
    "ncfile = open_file(pynio_filename)\n",
    "\n",
    "# Get the WRF variables\n",
    "slp = getvar(ncfile, \"slp\")\n",
    "smooth_slp = smooth2d(slp, 3)\n",
    "ctt = getvar(ncfile, \"ctt\")\n",
    "z = getvar(ncfile, \"z\", timeidx=0)\n",
    "dbz = getvar(ncfile, \"dbz\", timeidx=0)\n",
    "Z = 10**(dbz/10.)\n",
    "wspd =  getvar(ncfile, \"uvmet_wspd_wdir\", units=\"kt\")[0,:]\n",
    "\n",
    "# Set the start point and end point for the cross section\n",
    "start_point = CoordPair(lat=26.75, lon=-80.0)\n",
    "end_point = CoordPair(lat=26.75, lon=-77.8)\n",
    "\n",
    "# Compute the vertical cross-section interpolation.  Also, include the lat/lon points along the cross-section.\n",
    "z_cross = vertcross(Z, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True)\n",
    "wspd_cross = vertcross(wspd, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True)\n",
    "dbz_cross = 10.0 * np.log10(z_cross)\n",
    "\n",
    "# Extract the latitude and longitude coordinate arrays as regular numpy array instead of xarray.DataArray\n",
    "lons = npvalues(slp.coords[\"XLONG\"])\n",
    "lats = npvalues(slp.coords[\"XLAT\"])\n",
    "\n",
    "# Get the cartopy projection class\n",
    "wrf_proj = slp.attrs[\"projection\"]\n",
    "cart_proj = wrf_proj.cartopy()\n",
    "\n",
    "# Create the figure which will have 3 subplots\n",
    "fig = plt.figure(figsize=(7,5))\n",
    "ax_ctt = fig.add_subplot(1,2,1,projection=cart_proj)\n",
    "ax_wspd = fig.add_subplot(2,2,2)\n",
    "ax_dbz = fig.add_subplot(2,2,4)\n",
    "\n",
    "## Plot the cloud top temperature\n",
    "\n",
    "# Download and create the states, land, and oceans using cartopy features\n",
    "states = cfeature.NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n",
    "                             name='admin_1_states_provinces_shp')\n",
    "land = cfeature.NaturalEarthFeature(category='physical', name='land', scale='50m', \n",
    "                                    facecolor=cfeature.COLORS['land'])\n",
    "ocean = cfeature.NaturalEarthFeature(category='physical', name='ocean', scale='50m', \n",
    "                                     facecolor=cfeature.COLORS['water'])\n",
    "\n",
    "# Make the pressure contours.\n",
    "contour_levels = [960, 965, 970, 975, 980, 990]\n",
    "c1 = ax_ctt.contour(lons, lats, npvalues(smooth_slp), levels=contour_levels, colors=\"white\", \n",
    "            transform=crs.PlateCarree(), zorder=3, linewidths=1.0)\n",
    "\n",
    "# Create the filled cloud top temperature contours\n",
    "contour_levels = [-80.0, -70.0, -60, -50, -40, -30, -20, -10, 0, 10]\n",
    "ctt_contours = ax_ctt.contourf(lons, lats, npvalues(ctt), contour_levels, cmap=get_cmap(\"Greys\"),\n",
    "             transform=crs.PlateCarree(), zorder=2)\n",
    "\n",
    "ax_ctt.plot([start_point.lon, end_point.lon], [start_point.lat, end_point.lat], color=\"yellow\", \n",
    "            marker=\"o\", transform=crs.PlateCarree(), zorder=3)\n",
    "\n",
    "# Create the label bar for cloud top temperature\n",
    "cb_ctt = fig.colorbar(ctt_contours, ax=ax_ctt, shrink=.5)\n",
    "cb_ctt.ax.tick_params(labelsize=5)\n",
    "\n",
    "# Draw the oceans, land, and states\n",
    "ax_ctt.add_feature(ocean)\n",
    "ax_ctt.add_feature(land)\n",
    "ax_ctt.add_feature(states, linewidth=.5, edgecolor=\"black\")\n",
    "\n",
    "# Crop the domain to the region around the hurricane\n",
    "ax_ctt.set_extent([-85., -75.0, np.amin(lats), 30.0], crs=crs.PlateCarree())\n",
    "ax_ctt.gridlines(crs=crs.PlateCarree(), draw_labels=False)\n",
    "\n",
    "## Plot the cross sections\n",
    "\n",
    "# Make the contour plot for wspd\n",
    "wspd_contours = ax_wspd.contourf(npvalues(wspd_cross))\n",
    "# Add the color bar\n",
    "cb_wspd = fig.colorbar(wspd_contours, ax=ax_wspd)\n",
    "cb_wspd.ax.tick_params(labelsize=5)\n",
    "\n",
    "# Make the contour plot for dbz\n",
    "levels = [5 + 5*n for n in range(15)]\n",
    "dbz_contours = ax_dbz.contourf(npvalues(dbz_cross), levels=levels)\n",
    "cb_dbz = fig.colorbar(dbz_contours, ax=ax_dbz)\n",
    "cb_dbz.ax.tick_params(labelsize=5)\n",
    "\n",
    "# Set the x-ticks to use latitude and longitude labels.\n",
    "coord_pairs = npvalues(dbz_cross.coords[\"xy_loc\"])\n",
    "x_ticks = np.arange(coord_pairs.shape[0])\n",
    "x_labels = [pair.latlon_str() for pair in npvalues(coord_pairs)]\n",
    "ax_wspd.set_xticks(x_ticks[::20])\n",
    "ax_wspd.set_xticklabels([], rotation=45)\n",
    "ax_dbz.set_xticks(x_ticks[::20])\n",
    "ax_dbz.set_xticklabels(x_labels[::20], rotation=45, fontsize=4) \n",
    "\n",
    "\n",
    "# Set the y-ticks to be height.\n",
    "vert_vals = npvalues(dbz_cross.coords[\"vertical\"])\n",
    "v_ticks = np.arange(vert_vals.shape[0])\n",
    "ax_wspd.set_yticks(v_ticks[::20])\n",
    "ax_wspd.set_yticklabels(vert_vals[::20], fontsize=4) \n",
    "ax_dbz.set_yticks(v_ticks[::20])\n",
    "ax_dbz.set_yticklabels(vert_vals[::20], fontsize=4) \n",
    "\n",
    "# Set the x-axis and  y-axis labels\n",
    "ax_dbz.set_xlabel(\"Latitude, Longitude\", fontsize=5)\n",
    "ax_wspd.set_ylabel(\"Height (m)\", fontsize=5)\n",
    "ax_dbz.set_ylabel(\"Height (m)\", fontsize=5)\n",
    "\n",
    "# Add a title\n",
    "ax_ctt.set_title(\"Cloud Top Temperature (degC)\", {\"fontsize\" : 7})\n",
    "ax_wspd.set_title(\"Cross-Section of Wind Speed (kt)\", {\"fontsize\" : 7})\n",
    "ax_dbz.set_title(\"Cross-Section of Reflectivity (dBZ)\", {\"fontsize\" : 7})\n",
    "\n",
    "plt.savefig(\"/Users/ladwig/Documents/workspace/wrf_python/doc/source/_static/images/matthew.png\",\n",
    "           transparent=True, bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# SLP\n",
    "from __future__ import (absolute_import, division, print_function, unicode_literals)\n",
    "    \n",
    "from netCDF4 import Dataset   \n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.cm import get_cmap\n",
    "import cartopy.crs as crs\n",
    "from cartopy.feature import NaturalEarthFeature\n",
    "\n",
    "from wrf import npvalues, getvar, smooth2d\n",
    "\n",
    "ncfile = Dataset(\"/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00\")\n",
    "\n",
    "# Get the sea level pressure\n",
    "slp = getvar(ncfile, \"slp\")\n",
    "\n",
    "# Smooth the sea level pressure since it tends to be noisey near the mountains\n",
    "smooth_slp = smooth2d(slp, 3)\n",
    "\n",
    "# Get the numpy array from the XLAT and XLONG coordinates\n",
    "lats = npvalues(slp.coords[\"XLAT\"])\n",
    "lons = npvalues(slp.coords[\"XLONG\"])\n",
    "\n",
    "# Get the wrf.WrfProj object\n",
    "wrf_proj = slp.attrs[\"projection\"]\n",
    "\n",
    "# The cartopy() method returns a cartopy.crs projection object\n",
    "cart_proj = wrf_proj.cartopy()\n",
    "\n",
    "# Create a figure that's 10x10\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "# Get the GeoAxes set to the projection used by WRF\n",
    "ax = plt.axes(projection=cart_proj)\n",
    "\n",
    "# Download and add the states and coastlines\n",
    "states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n",
    "                             name='admin_1_states_provinces_shp')\n",
    "ax.add_feature(states, linewidth=.5)\n",
    "ax.coastlines('50m', linewidth=0.8)\n",
    "\n",
    "# Make the contour outlines and filled contours for the smoothed sea level pressure.\n",
    "# The transform keyword indicates that the lats and lons arrays are lat/lon coordinates and tells \n",
    "# cartopy to transform the points in to grid space.\n",
    "plt.contour(lons, lats, npvalues(smooth_slp), 10, colors=\"black\", transform=crs.PlateCarree())\n",
    "plt.contourf(lons, lats, npvalues(smooth_slp), 10, transform=crs.PlateCarree())\n",
    "\n",
    "# Add a color bar\n",
    "plt.colorbar(ax=ax, shrink=.47)\n",
    "\n",
    "# Set the map limits\n",
    "ax.set_xlim(wrf_proj.cartopy_xlim())\n",
    "ax.set_ylim(wrf_proj.cartopy_ylim())\n",
    "\n",
    "# Add the gridlines\n",
    "ax.gridlines()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}