FV3_FCST_skew-T.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examples of Plotting FV3 Forecast Files\n",
    "\n",
    "FV3 forecasts are written out in two separate files at each output time. Files prefixed `dynf` contain 3D dynamics variables, while files prefixed `phyf` contain physics forecast information, mainly at the surface."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This first cell is to run the common blocks before calling plotting, with input for the forecast hours - fhr. This way the whole notebook looks cleaner.\n",
    "\n",
    "Main functions are to import packages, read in and define data, define plotting functions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import os\n",
    "import sys\n",
    "from argparse import Namespace\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.basemap import Basemap, shiftgrid\n",
    "import numpy as np\n",
    "from netCDF4 import Dataset\n",
    "from matplotlib.colors import LinearSegmentedColormap\n",
    "from datetime import date, datetime, timedelta\n",
    "import glob\n",
    "import metpy.calc as mpcalc\n",
    "import metpy\n",
    "from metpy.plots import Hodograph, SkewT\n",
    "from metpy.units import pandas_dataframe_to_unit_arrays, units\n",
    "# plt.rcParams.update({'figure.max_open_warning': 0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Add units to the data arrays\n",
    "# p = p * units.mbar\n",
    "# T = T * units.degC\n",
    "# Td = Td * units.degC\n",
    "# spd = spd * units.knot\n",
    "# direc = direc * units.deg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "station='SLC'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initdate='2019043018'\n",
    "# fhr=6\n",
    "# stn={'LZK':(34.84, -92.26), 'SGF':(37.23, -93.38), 'OUN':(35.18, -97.44), 'FWD':(32.83, -97.3)}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initdate='2017121306'\n",
    "fhr=6\n",
    "stn={'SLC':(40.77, -111.95), 'LKN':(40.86, -115.73)}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get station lat/lon\n",
    "lat_obs, long_obs=stn[station]\n",
    "if long_obs < 0:\n",
    "    long_obs = 360 + long_obs \n",
    "print(lat_obs, long_obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inithr=initdate[8:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_validdate(initdate, fhr):\n",
    "    validdate=(datetime.strptime(initdate,\"%Y%m%d%H\")+timedelta(hours=fhr)).strftime(\"%Y%m%d%H\")\n",
    "    print(validdate)\n",
    "    validday=validdate[0:8]\n",
    "    print(validday)\n",
    "    validhr=validdate[8:10]\n",
    "    print(validhr)\n",
    "    return validdate, validday, validhr\n",
    "    \n",
    "validdate, validday, validhr=get_validdate(initdate, fhr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Provide a file path to a forecast directory. \n",
    "# The example below creates a dictionary containing 2 experiments, expt, through the first 4 forecast hours (including 0)\n",
    "\n",
    "exptlist=['k2n3.no_cu_gf', 'k6n2.no_cu_gf.L82_2mb']\n",
    "expt=exptlist[1]\n", 
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_path = f'/path/to/FV3SAR/{expt}/{initdate}'\n",
    "dynfilenc = os.path.join(file_path, f'dynf{fhr:03d}.nc')\n",
    "phyfilenc = os.path.join(file_path, f'phyf{fhr:03d}.nc')\n",
    "# print(f'dynf{fhr:03d}.nc')\n",
    "akbkfilenc = os.path.join(file_path, 'atmos_static.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dynfile = Dataset(dynfilenc,'r')\n",
    "phyfile = Dataset(phyfilenc,'r')\n",
    "akbkfile = Dataset(akbkfilenc,'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "akbkfile = Dataset(akbkfilenc, 'r')\n",
    "for v, info in akbkfile.variables.items():\n",
    "    print(v, ':', info.long_name, info.units, info.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "dynfile = Dataset(dynfilenc,'r')\n",
    "for v, info in dynfile.variables.items():\n",
    "    print(v, ':', info.long_name, info.units, info.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t3d=np.squeeze(dynfile['tmp'][::])\n",
    "u3d=np.squeeze(dynfile['ugrd'][::])\n",
    "v3d=np.squeeze(dynfile['vgrd'][::])\n",
    "psfc=np.squeeze(dynfile['pressfc'][::])\n",
    "spfh3d=np.squeeze(dynfile['spfh'][::])\n",
    "\n",
    "print(np.shape(t3d))\n",
    "print('T:', np.amin(t3d), np.amax(t3d))\n",
    "print(np.shape(u3d))\n",
    "print(np.shape(v3d))\n",
    "print(np.shape(spfh3d))\n",
    "print(np.shape(psfc))\n",
    "\n",
    "ny=np.shape(psfc)[0]\n",
    "nx=np.shape(psfc)[1]\n",
    "print(ny)\n",
    "print(nx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ak=np.squeeze(akbkfile['pk'][::])\n",
    "bk=np.squeeze(akbkfile['bk'][::])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This cell calculates pressure (phalf) based on ak + bk * Psfc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nlev=np.shape(ak)[0]\n",
    "pres = np.array([ak[i] + bk[i] * psfc for i in range(nlev)]) \n",
    "print(np.shape(pres))\n",
    "# print(pres[:,5,5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This cell calculates pfull from phalf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.shape(pres))\n",
    "presf = (pres[1:,:,:] + pres[0:-1, :,: ])/2\n",
    "print(np.shape(presf))\n",
    "print(presf[:,5,5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lat = np.squeeze(dynfile['grid_yt'][::]) * 180 / math.pi\n",
    "lon = np.squeeze(dynfile['grid_xt'][::]) * 180 / math.pi\n",
    "print(np.shape(lat))\n",
    "print(np.amax(lat), np.amin(lat))\n",
    "print(np.shape(lon))\n",
    "print(np.amax(lon), np.amin(lon))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This cell serves the important role of finding index of lat/lon in the dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y, x = np.unravel_index((np.abs(lat - lat_obs) + np.abs(lon - long_obs)).argmin(), lat.shape)\n",
    "print(y,x)\n",
    "print(lat[y,x], lon[y,x])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This cell get the 1-D profile of T, p, u, v, spfh at obs location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "T=t3d[:, y, x] - 273.15\n",
    "# print(T)\n",
    "T=T*units.degC\n",
    "p=presf[:, y, x]/100 * units.hPa\n",
    "# print(p)\n",
    "u=u3d[:, y, x]*1.94384 * units.knots\n",
    "v=v3d[:, y, x]*1.94384 * units.knots\n",
    "spfh=spfh3d[:, y, x] * units.kg/units.kg\n",
    "# print(u)\n",
    "# print(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculates dewpoint temperature.\n",
    "metpy.calc.dewpoint_from_specific_humidity(pressure, temperature, specific_humidity)\n",
    "https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.dewpoint_from_specific_humidity.html\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Td=mpcalc.dewpoint_from_specific_humidity(pressure=p,temperature=T,specific_humidity=spfh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_skewT(p, T, Td, u, v):\n",
    "    # Change default to be better for skew-T\n",
    "    fig = plt.figure(figsize=(9, 11))\n",
    "\n",
    "    # Initiate the skew-T plot type from MetPy class loaded earlier\n",
    "    skew = SkewT(fig, rotation=45)\n",
    "\n",
    "    # Plot the data using normal plotting functions, in this case using\n",
    "    # log scaling in Y, as dictated by the typical meteorological plot\n",
    "    skew.plot(p, T, 'r')\n",
    "    skew.plot(p, Td, 'g')\n",
    "#     skew.plot_barbs(p[::3], u[::3], v[::3], y_clip_radius=0.03)\n",
    "    skew.plot_barbs(p, u, v, y_clip_radius=0.03)\n",
    "\n",
    "    # Set some appropriate axes limits for x and y\n",
    "    skew.ax.set_xlim(-30, 40)\n",
    "    skew.ax.set_ylim(1020, 100)\n",
    "\n",
    "    # Add the relevant special lines to plot throughout the figure\n",
    "    skew.plot_dry_adiabats(t0=np.arange(233, 533, 10) * units.K,\n",
    "                           alpha=0.25, color='orangered')\n",
    "    skew.plot_moist_adiabats(t0=np.arange(233, 400, 5) * units.K,\n",
    "                             alpha=0.25, color='tab:green')\n",
    "    skew.plot_mixing_lines(p=np.arange(1000, 99, -20) * units.hPa,\n",
    "                           linestyle='dotted', color='tab:blue')\n",
    "\n",
    "    # Add some descriptive titles\n",
    "#     plt.title('{} Sounding'.format(station), loc='left')\n",
    "    plt.title(f'Sounding at {station}  (lat={lat_obs}, lon={np.around(long_obs-360,3)}) \\n')\n",
    "    plt.title(f'Valid Time: {validdate}     ', loc='right')\n",
    "    plt.title(expt, loc='left')\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The following sections make the plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_skewT(p, T, Td, u, v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}