{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# NCL_WRF_zoom_1_2.py\nThis script illustrates the following concepts:\n    - Plotting WRF data on native grid\n    - Subsetting data to 'zoom in' on an area\n    - Plotting data using wrf python functions\n    - Following best practices when choosing a colormap\n    \nSee following URLs to see the reproduced NCL plot & script:\n    - Original NCL script: https://www.ncl.ucar.edu/Applications/Scripts/wrf_zoom_1.ncl\n    - Original NCL plot: https://www.ncl.ucar.edu/Applications/Images/wrf_zoom_1_2_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 matplotlib.ticker as mticker\nimport cartopy.crs as ccrs\nfrom cartopy.feature import NaturalEarthFeature\n\nfrom wrf import (getvar, to_np, latlon_coords, get_cartopy)\nimport geocat.datafiles as gdf"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Read in the data\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "wrfin = Dataset(gdf.get(\"netcdf_files/wrfout_d03_2012-04-22_23_00_00_subset.nc\"))\ntd2 = getvar(wrfin, \"td2\")\n\n# Set attributes for creating plot titles later\nwe = getattr(wrfin, 'WEST-EAST_GRID_DIMENSION')\nsn = getattr(wrfin, 'SOUTH-NORTH_GRID_DIMENSION')\nlvl = getattr(wrfin, 'BOTTOM-TOP_GRID_DIMENSION')\ndis = getattr(wrfin, 'DY')/1000 # Divide by 1000 to go from m to km\nphys = getattr(wrfin,'MP_PHYSICS')\npbl = getattr(wrfin, 'BL_PBL_PHYSICS')\ncu = getattr(wrfin, 'CU_PHYSICS')\ns_date = getattr(wrfin, 'START_DATE') \nstr_format = \"WE={}; SN={}; Levels={}; Dis={}km; Phys Opt={}; PBL Opt={}; Cu Opt={}\"\nsd_frmt = \"Init: {}\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Create a subset of the data for zoomed in projection\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dims = td2.shape\n\ny_start = int(dims[0]/2)\ny_end = int(dims[0]-1)\n\nx_start = int(0)\nx_end = int(dims[1]/2)\n\ntd2_zoom = td2[y_start:y_end, x_start:x_end]\n\n# Define the latitude and longitude coordinates\nlats, lons = latlon_coords(td2_zoom)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the data\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# The `get_cartopy` wrf function will automatically find and use the \n# intended map projection for this dataset \ncart_proj = get_cartopy(td2_zoom)\n\nfig = plt.figure(figsize=(12,12))\nax = plt.axes(projection=cart_proj)\n\n# Add features to the projection\nstates = NaturalEarthFeature(category=\"cultural\",\n                             scale=\"50m\",\n                             facecolor=\"none\",\n                             name=\"admin_1_states_provinces_shp\")\n\nax.add_feature(states, linewidth=0.5, edgecolor=\"black\")\nax.coastlines('50m', linewidth=0.8)\n\n# Add filled dew point temperature contours\nplt.contourf(to_np(lons),\n             to_np(lats),\n             to_np(td2_zoom),\n             levels=13, cmap=\"magma\",\n             transform=ccrs.PlateCarree(),\n             vmin=-8,\n             vmax=18)\n\n# Add a colorbar\ncbar = plt.colorbar(ax=ax,\n                    orientation=\"horizontal\",\n                    ticks=np.arange(-6,18,2),\n                    drawedges=True,\n                    extendrect=True,\n                    pad=0.08,\n                    shrink=0.75,\n                    aspect=30)\n\n# Format location of colorbar text to look like NCL version\ncbar.ax.text(0.5,\n             1.5,\n             '2m Dewpoint Temperature (C)',\n             horizontalalignment='center',\n             verticalalignment='center',\n             transform=cbar.ax.transAxes)\n\n# Format colorbar ticks and labels \ncbar.ax.tick_params(labelsize=10)\ncbar.ax.get_xaxis().labelpad = -48\n\n# Draw gridlines\ngl = ax.gridlines(crs=ccrs.PlateCarree(),\n                  draw_labels=True,\n                  dms=False,\n                  x_inline=False,\n                  y_inline=False,\n                  linewidth=1,\n                  color=\"k\",\n                  alpha=0.25)\n\n# Manipulate latitude and longitude gridline numbers and spacing\ngl.top_labels = False\ngl.right_labels = False\ngl.xlocator = mticker.FixedLocator([-120, -121, -122, -123, -124])\ngl.ylocator = mticker.FixedLocator([38, 39, 40, 41, 42])\ngl.xlabel_style = {\"rotation\": 0, \"size\": 15}\ngl.ylabel_style = {\"rotation\": 0, \"size\": 15}\ngl.xlines = False\ngl.ylines = False\n\n# Add titles to the plot\nplt.title(\"Zoomed in plot\", loc='center', x=.13, y=1.1, size=15)\nplt.title(\"2m Dewpoint Temperature (C)\", loc='left', y=1.02, size=10)\nplt.title(sd_frmt.format(s_date), loc='right', y=1.1, size=10)\n\n# Add lower text using attributes from the dataset\nfig.text(0.25, 0.1, getattr(wrfin, 'TITLE'), size=12)\nfig.text(0.252, 0.08, \n         str_format.format(we, sn, lvl, dis, phys, pbl, cu),\n         size=12)\n\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
}