{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# NCL_WRF_interp_3.py\nThis script illustrates the following concepts:\n    - Interpolating a vertical cross-section from a 3D WRF-ARW field.\n    - Recombining two datasets into one usable form\n    - Following best practices when choosing a colormap.\n      More information on colormap best practices can be found `here <https://geocat-examples.readthedocs.io/en/latest/gallery/Colors/CB_Temperature.html#sphx-glr-gallery-colors-cb-temperature-py>`_.\n    \nSee following URLs to see the reproduced NCL plot & script:\n    - Original NCL script: https://www.ncl.ucar.edu/Applications/Scripts/wrf_interp_3.ncl\n    - Original NCL plot: https://www.ncl.ucar.edu/Applications/Images/wrf_interp_3_lg.png\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Import packages\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from netCDF4 import Dataset\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport xarray as xr\nimport os\n\nfrom wrf import (to_np, getvar, CoordPair, vertcross, latlon_coords)\nimport geocat.datafiles as gdf\nimport geocat.viz.util as gvutil"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Read in the data\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Specify the necessary variables needed from the data set in order to use 'z' and 'QVAPOR'\ntoinclude = ['PH', 'P', 'HGT', 'PHB', 'QVAPOR'] \n# Read in necessary datasets\nds = xr.open_mfdataset([gdf.get('netcdf_files/wrfout_d03_2012-04-22_23_00_00_Z.nc'),\n                        gdf.get('netcdf_files/wrfout_d03_2012-04-22_23_00_00_QV.nc')])\n\n# specify a unique output file name to use to read in combined dataset later\nfile3 = 'wrfout_d03_2012-04-22_23.nc'\nmrg = ds[toinclude].to_netcdf(file3)\n\n# Read in the data and extract variables \nwrfin = Dataset(('wrfout_d03_2012-04-22_23.nc'))\n\nz = getvar(wrfin, \"z\")\nqv = getvar(wrfin, \"QVAPOR\")\n# Pull lat/lon coords from QVAPOR data using wrf-python tools\nlats, lons = latlon_coords(qv)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Create vertical cross section using wrf-python tools\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Define start and stop coordinates for cross section\nstart_point = CoordPair(lat=38, lon=-118)\nend_point = CoordPair(lat=40, lon=-115)\n\nqv_cross = vertcross(qv, \n                     z, \n                     wrfin=wrfin,\n                     start_point=start_point, \n                     end_point=end_point,\n                     latlon=True)\n\n# Close 'wrfin' to prevent PermissionError if code is run more than once locally\nwrfin.close()\n# Remove created wrfout file from local directory \nos.remove('wrfout_d03_2012-04-22_23.nc')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the data\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(10,8))\nax = plt.axes()\n\n# Set the x-ticks to use latitude and longitude labels.\ncoord_pairs = to_np(qv_cross.coords[\"xy_loc\"])\nx_ticks = np.arange(coord_pairs.shape[0])\n\n# Plot filled contours\nqv_contours = qv_cross.plot.contourf(ax=ax,\n                                     levels=17,\n                                     cmap='magma',\n                                     vmin=0,\n                                     vmax=0.004,\n                                     zorder=4,\n                                     add_labels=False,\n                                     add_colorbar=False,\n                                     yticks=np.arange(0, 20000, 3000),\n                                     xticks=x_ticks[::20])\n# Add colorbar\nplt.colorbar(qv_contours, \n             ax=ax, \n             ticks=np.arange(0.00025, 0.004, .00025))\n\n# Add minor ticks to the yaxis\ngvutil.add_major_minor_ticks(ax=ax, \n                             x_minor_per_major=1, \n                             y_minor_per_major=3, \n                             labelsize=14)\n\n# Format the xtick labels\nx_labels = [pair.latlon_str(fmt=\"{:.2f}\\N{DEGREE SIGN}N, \\n {:.2f}\\N{DEGREE SIGN}E\")\n            for pair in to_np(coord_pairs)]\nax.set_xticklabels(x_labels[::20], rotation=45, fontsize=12)\n\n# Set the plot titles \nplt.title(\"Cross section from (38,-118) to (40,-115)\", fontsize=18, y=1.07)\nplt.title('Water vapor mixing ratio', loc='left', y=1.02)\nplt.title('kg kg-1', loc='right', y=1.02)\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
}