{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"http://xarray.pydata.org/en/stable/_static/dataset-diagram-logo.png\" align=\"center\" width=\"30%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 02_Xarray: Plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "1. Work with an Xarray `Dataset`\n",
    "2. Select a variable from the dataset\n",
    "3. Create a contour plot of gridded CFSR reanalysis data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from metpy.units import units\n",
    "import metpy.calc as mpcalc\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Work with an Xarray `Dataset`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An Xarray [Dataset](http://xarray.pydata.org/en/stable/generated/xarray.Dataset.html) consists of one or more `DataArray`s. Let's continue to work with the CFSR NetCDF files from the previous notebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "year = 1993\n",
    "ds_slp = xr.open_dataset(f'/cfsr/data/{year}/pmsl.1993.0p5.anl.nc')\n",
    "ds_g = xr.open_dataset(f'/cfsr/data/{year}/g.1993.0p5.anl.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examine the SLP dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_slp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One attribute of a dataset is its size. We can access that via its `nbytes` attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (f'Size of dataset: {ds_slp.nbytes / 1e6} MB')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select a variable from the dataset and assign it to a `DataArray` object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slp = ds_slp['pmsl'] # Similar to selecting a Series in Pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that `slp` is a 3-dimensional DataArray, with time ('time') as the first dimension."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a contour plot of gridded data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We got a quick plot in the previous notebook. Let's now customize the map, and focus on a particular region. We will use Cartopy and Matplotlib to plot our SLP contours. Matplotlib has two contour methods:\n",
    "1. [Line contours: ax.contour](https://matplotlib.org/api/_as_gen/matplotlib.axes.Axes.contour.html)\n",
    "2. [Filled contours: ax.contourf](https://matplotlib.org/api/_as_gen/matplotlib.axes.Axes.contourf.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will draw contour lines in hPa, with a contour interval of 4."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### As we've done before, let's first define some variables relevant to Cartopy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The CFSR data's x- and y- coordinates are longitude and latitude ... thus, its projection is Plate Carree. However, we wish to display the data on a regional map that uses Lambert Conformal. So we will define two Cartopy projection objects here. We will need to *transform* the data from its native projection to the map's projection when it comes time to draw the contours."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cLon = -80\n",
    "cLat = 35\n",
    "lonW = -100\n",
    "lonE = -60\n",
    "latS = 20\n",
    "latN = 50\n",
    "proj_data = ccrs.PlateCarree() # Our data is lat-lon; thus its native projection is Plate Carree.\n",
    "proj_map = ccrs.LambertConformal(central_longitude=cLon, central_latitude=cLat) # Map projection\n",
    "res = '50m'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Now define the range of our contour values and a contour interval. 4 hPa is standard for a sea-level pressure map on a synoptic-scale region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "minVal = 900\n",
    "maxVal = 1076\n",
    "cint = 4\n",
    "cintervals = np.arange(minVal, maxVal, cint)\n",
    "cintervals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Matplotlib's contour methods require *three arrays* to be passed to them ... x- and y- arrays (longitude and latitude in our case), and a 2-d array (corresponding to x- and y-) of our data variable. So we need to extract the latitude and longitude coordinate variables from our DataArray. We'll also extract the third coordinate value, time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lats = ds_slp.lat\n",
    "lons = ds_slp.lon\n",
    "times = ds_slp.time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use Xarray's `sel` method to specify one time from the `DataArray`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slp_singleTime = slp.sel(time='1993-03-14-18:00:00')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slp_singleTime"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the units are in Pascals. We will exploit MetPy's unit conversion library soon, but for now let's just divide the array by 100."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slp_singleTimeHPA = slp_singleTime / 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### We're set to make our map. After creating our figure, setting the bounds of our map domain, and adding cartographic features, we will plot the contours. We'll assign the output of the contour method to an object, which will then be used to label the contour lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(18,12))\n",
    "ax = plt.subplot(1,1,1,projection=proj_map)\n",
    "ax.set_extent ([lonW,lonE,latS,latN])\n",
    "ax.add_feature(cfeature.COASTLINE.with_scale(res))\n",
    "ax.add_feature(cfeature.STATES.with_scale(res))\n",
    "CL = ax.contour(lons,lats,slp_singleTimeHPA,cintervals,transform=proj_data,linewidths=1.25,colors='green')\n",
    "ax.clabel(CL, inline_spacing=0.2, fontsize=11, fmt='%.0f');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### There we have it! There are definitely some things we can improve about this plot (such as the lack of an informative title, and the fact that our contour labels seem not to be appearing on every contour line), but we'll get to that in the next lesson!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 August 2022 Environment",
   "language": "python",
   "name": "aug22"
  },
  "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}