bulk_shear_calc_issue.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate Bulk Wind Shear for JSON Sounding Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import metpy.calc as mpcalc\n",
    "from metpy.units import units\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "from metpy.calc.tools import get_layer, _get_bound_pressure_height\n",
    "from metpy.calc.tools import _less_or_close, _greater_or_close, log_interpolate_1d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pres</th>\n",
       "      <th>hght</th>\n",
       "      <th>tmpc</th>\n",
       "      <th>dwpt</th>\n",
       "      <th>wdir</th>\n",
       "      <th>wspd</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>[[995, 984, 972, 939.7, 925, 911, 908.4, 877.6...</td>\n",
       "      <td>[[46, 93, 194, 305, 610, 752, 888, 914, 1219, ...</td>\n",
       "      <td>[[36, 35, 34, 31.1, 29.8, 28.6, 28.4, 26.1, 24...</td>\n",
       "      <td>[[23, 25, 24.4, 22.6, 21.8, 21.6, 21.4, 19.1, ...</td>\n",
       "      <td>[[210, 210, 210, 210, 215, 223, 225, 270, 280,...</td>\n",
       "      <td>[[7, 9, 12, 12, 10, 8, 8, 9, 12, 15, 18, 21, 2...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                pres  \\\n",
       "0  [[995, 984, 972, 939.7, 925, 911, 908.4, 877.6...   \n",
       "\n",
       "                                                hght  \\\n",
       "0  [[46, 93, 194, 305, 610, 752, 888, 914, 1219, ...   \n",
       "\n",
       "                                                tmpc  \\\n",
       "0  [[36, 35, 34, 31.1, 29.8, 28.6, 28.4, 26.1, 24...   \n",
       "\n",
       "                                                dwpt  \\\n",
       "0  [[23, 25, 24.4, 22.6, 21.8, 21.6, 21.4, 19.1, ...   \n",
       "\n",
       "                                                wdir  \\\n",
       "0  [[210, 210, 210, 210, 215, 223, 225, 270, 280,...   \n",
       "\n",
       "                                                wspd  \n",
       "0  [[7, 9, 12, 12, 10, 8, 8, 9, 12, 15, 18, 21, 2...  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_json(\"data_IAD.json\")\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = df['pres'].values[0][0] * units.hPa\n",
    "t = df['tmpc'].values[0][0] * units.degC\n",
    "td = df['dwpt'].values[0][0] * units.degC\n",
    "dir = df['wdir'].values[0][0] * units.degrees\n",
    "spd = df['wspd'].values[0][0] * units.knots\n",
    "heights = df['hght'].values[0][0] * units.meter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "u, v = mpcalc.wind_components(df[\"wspd\"].values[0][0] * units.knot, df[\"wdir\"].values[0][0] * units.deg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quick function to get bulk shear at desired layer (meters) for JSON files\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_bulk_shear(sound_file,depth):\n",
    "    \"\"\"Get bulk shear for desired layer depth based on JSON sounding file\n",
    "    \n",
    "    Args\n",
    "    ----\n",
    "    sound_file : str\n",
    "        JSON sounding file name\n",
    "    \n",
    "    depth : int\n",
    "        layer depth desired in meters\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    Prints u, v, speed, and direction for bilk shear values\n",
    "    \n",
    "    u_bulk_shear : pint.quantity.build_quantity_class.<locals>.Quantity\n",
    "        u-component of layer bulk shear\n",
    "        \n",
    "    v_bulk_shear : pint.quantity.build_quantity_class.<locals>.Quantity\n",
    "        v-component of layer bulk shear\n",
    "        \n",
    "    bulk_shear_speed : pint.quantity.build_quantity_class.<locals>.Quantity\n",
    "        layer bulk shear wind speed\n",
    "        \n",
    "    bulk_shear_dir : pint.quantity.build_quantity_class.<locals>.Quantity\n",
    "        layer bulk shear wind direction\n",
    "        \n",
    "    \"\"\"\n",
    "    printmd(f\"\\n**Sounding Location: {sound_file}**\")\n",
    "    print(f\"Desired layer: {depth/1000}km\\n\"+\\\n",
    "         \"---------------------------------\")\n",
    "    \n",
    "    df = pd.read_json(sound_file)\n",
    "    p = df['pres'].values[0][0] * units.hPa\n",
    "    Z = df['hght'].values[0][0] * units.meter\n",
    "    \n",
    "    def replace_empty_str(col):\n",
    "        for i in range(len(df[col][0][0][:])):\n",
    "            if df[col][0][0][i] == '':\n",
    "                df[col][0][0][i] = 0\n",
    "        return df\n",
    "    \n",
    "    for i in df.columns:\n",
    "        replace_empty_str(i)\n",
    "        \n",
    "    u, v = mpcalc.wind_components(df[\"wspd\"].values[0][0] * units.knot, df[\"wdir\"].values[0][0] * units.deg)\n",
    "    u_bulk_shear,v_bulk_shear = mpcalc.bulk_shear(p,u,v,heights=Z,depth=depth * units.meter)\n",
    "    print(f\"u-bulk shear: {u_bulk_shear}\\nv-bulk shear: {v_bulk_shear}\")\n",
    "    \n",
    "    bulk_shear_speed = np.sqrt(u_bulk_shear**2 + v_bulk_shear**2)\n",
    "    bulk_shear_dir = mpcalc.wind_direction(u_bulk_shear,v_bulk_shear)\n",
    "    print(f\"bulk shear speed: {bulk_shear_speed}\\nbulk shear direction: {bulk_shear_dir}\")\n",
    "    \n",
    "    return u_bulk_shear, v_bulk_shear, bulk_shear_speed, bulk_shear_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 322,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_layer(pressure, *args, heights=None, bottom=None, depth=100 * units.hPa,\n",
    "              interpolate=True):\n",
    "    r\"\"\"Return an atmospheric layer from upper air data with the requested bottom and depth.\n",
    "\n",
    "    This function will subset an upper air dataset to contain only the specified layer. The\n",
    "    bottom of the layer can be specified with a pressure or height above the surface\n",
    "    pressure. The bottom defaults to the surface pressure. The depth of the layer can be\n",
    "    specified in terms of pressure or height above the bottom of the layer. If the top and\n",
    "    bottom of the layer are not in the data, they are interpolated by default.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    pressure : array-like\n",
    "        Atmospheric pressure profile\n",
    "    args : array-like\n",
    "        Atmospheric variable(s) measured at the given pressures\n",
    "    heights: array-like, optional\n",
    "        Atmospheric heights corresponding to the given pressures. Defaults to using\n",
    "        heights calculated from ``p`` assuming a standard atmosphere [NOAA1976]_.\n",
    "    bottom : `pint.Quantity`, optional\n",
    "        The bottom of the layer as a pressure or height above the surface pressure. Defaults\n",
    "        to the highest pressure or lowest height given.\n",
    "    depth : `pint.Quantity`, optional\n",
    "        The thickness of the layer as a pressure or height above the bottom of the layer.\n",
    "        Defaults to 100 hPa.\n",
    "    interpolate : bool, optional\n",
    "        Interpolate the top and bottom points if they are not in the given data. Defaults\n",
    "        to True.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    `pint.Quantity, pint.Quantity`\n",
    "        The pressure and data variables of the layer\n",
    "\n",
    "    \"\"\"\n",
    "    # If we get the depth kwarg, but it's None, set it to the default as well\n",
    "    if depth is None:\n",
    "        depth = 100 * units.hPa\n",
    "\n",
    "    # Make sure pressure and datavars are the same length\n",
    "    for datavar in args:\n",
    "        if len(pressure) != len(datavar):\n",
    "            raise ValueError('Pressure and data variables must have the same length.')\n",
    "\n",
    "    # If the bottom is not specified, make it the surface pressure\n",
    "    if bottom is None:\n",
    "        bottom = np.nanmax(pressure.m) * pressure.units\n",
    "\n",
    "    bottom_pressure, bottom_height = _get_bound_pressure_height(pressure, bottom,\n",
    "                                                                heights=heights,\n",
    "                                                                interpolate=interpolate)\n",
    "\n",
    "    # Calculate the top if whatever units depth is in\n",
    "    if depth.dimensionality == {'[length]': -1.0, '[mass]': 1.0, '[time]': -2.0}:\n",
    "        top = bottom_pressure - depth\n",
    "    elif depth.dimensionality == {'[length]': 1}:\n",
    "        top = bottom_height + depth\n",
    "    else:\n",
    "        raise ValueError('Depth must be specified in units of length or pressure')\n",
    "\n",
    "    top_pressure, _ = _get_bound_pressure_height(pressure, top, heights=heights,\n",
    "                                                 interpolate=interpolate)\n",
    "\n",
    "    ret = []  # returned data variables in layer\n",
    "\n",
    "    # Ensure pressures are sorted in ascending order\n",
    "    sort_inds = np.argsort(pressure)\n",
    "    pressure = pressure[sort_inds]\n",
    "\n",
    "    # Mask based on top and bottom pressure\n",
    "    inds = (_less_or_close(pressure, bottom_pressure)\n",
    "            & _greater_or_close(pressure, top_pressure))\n",
    "    p_interp = pressure[inds]\n",
    "\n",
    "    # Interpolate pressures at bounds if necessary and sort\n",
    "    if interpolate:\n",
    "        # If we don't have the bottom or top requested, append them\n",
    "        if not np.any(np.isclose(top_pressure, p_interp)):\n",
    "            p_interp = np.sort(np.append(p_interp.m, top_pressure.m)) * pressure.units\n",
    "        if not np.any(np.isclose(bottom_pressure, p_interp)):\n",
    "            p_interp = np.sort(np.append(p_interp.m, bottom_pressure.m)) * pressure.units\n",
    "\n",
    "    ret.append(p_interp[::-1])\n",
    "\n",
    "    for datavar in args:\n",
    "        # Ensure that things are sorted in ascending order\n",
    "        datavar = datavar[sort_inds]\n",
    "\n",
    "        if interpolate:\n",
    "            # Interpolate for the possibly missing bottom/top values\n",
    "            datavar_interp = log_interpolate_1d(p_interp, pressure, datavar)\n",
    "            datavar = datavar_interp\n",
    "        else:\n",
    "            datavar = datavar[inds]\n",
    "\n",
    "        ret.append(datavar[::-1])\n",
    "    return ret"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from metpy import constants as mpconsts\n",
    "\n",
    "t0 = 288. * units.kelvin\n",
    "p0 = 1013.25 * units.hPa\n",
    "\n",
    "def height_to_pressure_std(height):\n",
    "    r\"\"\"Convert height data to pressures using the U.S. standard atmosphere [NOAA1976]_.\n",
    "\n",
    "    The implementation inverts the formula outlined in [Hobbs1977]_ pg.60-61.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    height : `pint.Quantity`\n",
    "        Atmospheric height\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    `pint.Quantity`\n",
    "        The corresponding pressure value(s)\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    .. math:: p = p_0 e^{\\frac{g}{R \\Gamma} \\text{ln}(1-\\frac{Z \\Gamma}{T_0})}\n",
    "\n",
    "    \"\"\"\n",
    "    \n",
    "    gamma = 6.5 * units('K/km')\n",
    "    return p0 * (1 - (gamma / t0) * height) ** (mpconsts.g / (mpconsts.Rd * gamma))\n",
    "\n",
    "def pressure_to_height_std(pressure):\n",
    "    r\"\"\"Convert pressure data to heights using the U.S. standard atmosphere [NOAA1976]_.\n",
    "\n",
    "    The implementation uses the formula outlined in [Hobbs1977]_ pg.60-61.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    pressure : `pint.Quantity`\n",
    "        Atmospheric pressure\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    `pint.Quantity`\n",
    "        The corresponding height value(s)\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    .. math:: Z = \\frac{T_0}{\\Gamma}[1-\\frac{p}{p_0}^\\frac{R\\Gamma}{g}]\n",
    "\n",
    "    \"\"\"\n",
    "    gamma = 6.5 * units('K/km')\n",
    "    return (t0 / gamma) * (1 - (pressure / p0).to('dimensionless')**(\n",
    "        mpconsts.Rd * gamma / mpconsts.g))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([487]) <Unit('hectopascal')>, 6000 <Unit('meter')>)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_get_bound_pressure_height(p, 6000*units.m, heights=heights, interpolate=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "**Top of 6km layer (from get_layer)**"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==================================\n",
      "Pressure: 486.0 hectopascal\n",
      "Height (from sounding): 6083.279983975857 meter\n",
      "Height (from std atm): 5.780369569654552 kilometer\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/markdown": [
       "**Expected top of 6km layer**"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==========================\n",
      "Pressure bounds from _get_bound_pressure_height: 468.6454879416241 hectopascal\n",
      "Height bounds from _get_bound_pressure_height: 6046 meter\n",
      "Base Height (from std atm): 46 meter\n",
      "6km Pressure (from std atm): 468.6454879416241 hectopascal\n",
      "6km-Pressure Height (from std atm): 6.046000000000002 kilometer\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/markdown": [
       "**Bottom of 6km layer**"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=================\n",
      "Pressure: 995.0 hectopascal\n",
      "Heights (from sounding): 46.0 meter\n",
      "Heights (from std atm): 0.15295996804345655 kilometer\n",
      "\n",
      "U-Shear component: 21.819830420597967 knot\n",
      "V-Shear component: -27.258009986877884 knot\n",
      "Bulk Shear Speed: 34.91567138160724 knot\n"
     ]
    },
    {
     "data": {
      "text/markdown": [
       "\n",
       "**Sounding Location: data_IAD.json**"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Desired layer: 6.0km\n",
      "---------------------------------\n",
      "u-bulk shear: 21.819830420597967 knot\n",
      "v-bulk shear: -27.258009986877884 knot\n",
      "bulk shear speed: 34.91567138160724 knot\n",
      "bulk shear direction: 321.32299241704106 degree\n"
     ]
    }
   ],
   "source": [
    "p_layer, t_layer, u_layer, v_layer, dir_layer, spd_layer, hghts_layer = get_layer(\n",
    "    p, t, u, v, dir, spd, heights, heights=heights, depth=6000*units.meter)\n",
    "\n",
    "printmd(\"**Top of 6km layer (from get_layer)**\")\n",
    "print(\"==================================\")\n",
    "print(\"Pressure: \" + str(p_layer[-1]))\n",
    "print(\"Height (from sounding): \" + str(hghts_layer[-1]))\n",
    "print(\"Height (from std atm): \" + str(pressure_to_height_std(p_layer[-1])))\n",
    "print()\n",
    "\n",
    "printmd(\"**Expected top of 6km layer**\")\n",
    "print(\"==========================\")\n",
    "#baseHght = pressure_to_height_std(p[0])\n",
    "baseHght = heights[0]\n",
    "pressure_bound, height_bound = _get_bound_pressure_height(p, 6000*units.meter+baseHght)\n",
    "stdPres = height_to_pressure_std(6000*units.meter+baseHght)\n",
    "stdHght = pressure_to_height_std(stdPres)\n",
    "print(\"Pressure bounds from _get_bound_pressure_height: \" + str(pressure_bound))\n",
    "print(\"Height bounds from _get_bound_pressure_height: \" + str(height_bound))\n",
    "print(\"Base Height (from std atm): \" + str(baseHght))\n",
    "print(\"6km Pressure (from std atm): \" + str(stdPres))\n",
    "print(\"6km-Pressure Height (from std atm): \" + str(stdHght))\n",
    "print()\n",
    "\n",
    "printmd(\"**Bottom of 6km layer**\")\n",
    "print(\"=================\")\n",
    "print(\"Pressure: \" + str(p_layer[0]))\n",
    "print(\"Heights (from sounding): \" + str(hghts_layer[0]))\n",
    "print(\"Heights (from std atm): \" + str(pressure_to_height_std(p_layer[0])))\n",
    "print()\n",
    "\n",
    "ushr = u_layer[-1] - u_layer[0]\n",
    "vshr = v_layer[-1] - v_layer[0]\n",
    "\n",
    "print(\"U-Shear component: \" + str(ushr))\n",
    "print(\"V-Shear component: \" + str(vshr))\n",
    "print(\"Bulk Shear Speed: \" + str(np.sqrt(ushr**2 + vshr**2)))\n",
    "\n",
    "u_bulk_shear, v_bulk_shear, bulk_shear_speed, bulk_shear_dir = get_bulk_shear(\"data_IAD.json\", 6000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Markdown, display\n",
    "def printmd(string):\n",
    "    display(Markdown(string))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 452,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_BRN(pressure, u, v, temp, cape, heights=None):\n",
    "    r\"\"\"Calculate Bulk Richardson Number and BRN shear.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    pressure : `pint.Quantity`\n",
    "                Atmospheric pressure profile\n",
    "    u        : `pint.Quantity`\n",
    "                U-component of wind.\n",
    "    v        : `pint.Quantity`\n",
    "                V-component of wind.\n",
    "    temp     : `pint.Quantity`\n",
    "                Atmospheric temperature profile (can be either temperature or virtual temperature)\n",
    "    cape     : `pint.Quantity`\n",
    "                CAPE value to use as the numerator of the BRN calculation.\n",
    "    heights  : `pint.Quantity`, optional\n",
    "                Heights in meters from sounding (not AGL)\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    brn     : `pint.Quantity`\n",
    "               Bulk Richadson Number\n",
    "    brn_shr : `pint.Quantity`\n",
    "               Bulk Richardson Shear (the denominator of the BRN)\n",
    "\n",
    "    \"\"\"\n",
    "    \n",
    "    from metpy import constants as mpconsts\n",
    "    \n",
    "    baseZ = heights[0]\n",
    "    sumZ = -baseZ\n",
    "    rho6km = 0\n",
    "    u6km = 0\n",
    "    v6km = 0\n",
    "    u500 = 0\n",
    "    v500 = 0\n",
    "\n",
    "    for i, p in enumerate(pressure):\n",
    "        rho = (p.to('Pa')/(temp[i].to('K')*mpconsts.Rd)).to('kg / m^3')\n",
    "        print(rho)\n",
    "        sumZ += heights[i]\n",
    "        u_weighted = u[i]*rho\n",
    "        v_weighted = v[i]*rho\n",
    "        u_weighted = u_weighted.magnitude\n",
    "        v_weighted = v_weighted.magnitude\n",
    "        \n",
    "        if (sumZ >= 6000 * units.meter):\n",
    "            layerZ=heights[i] - heights[i-1]\n",
    "            fraction = (6000 * units.meter - layerZ)/sumZ\n",
    "            fraction = fraction.to('dimensionless')\n",
    "            u_frac = ((u_weighted - uprev) * fraction) + uprev\n",
    "            v_frac = ((v_weighted - vprev) * fraction) + vprev\n",
    "\n",
    "            rho6km += rho\n",
    "            u_weighted = u_frac*rho\n",
    "            v_weighted = v_frac*rho\n",
    "            u_weighted = u_weighted.magnitude\n",
    "            v_weighted = v_weighted.magnitude\n",
    "            u6km += u_weighted\n",
    "            v6km += v_weighted\n",
    "            break\n",
    "\n",
    "        rho6km += rho\n",
    "        u6km += u_weighted\n",
    "        v6km += v_weighted\n",
    "\n",
    "        if (sumZ < 500 * units.meter):\n",
    "            u500 += u_weighted\n",
    "            v500 += v_weighted\n",
    "            divisor = i+1\n",
    "\n",
    "        uprev = u_weighted\n",
    "        vprev = v_weighted\n",
    "\n",
    "    u6kmAvg = u6km/(i+1)\n",
    "    v6kmAvg = v6km/(i+1)\n",
    "    print(u6kmAvg, v6kmAvg)\n",
    "    u500Avg = u500/divisor\n",
    "    v500Avg = v500/divisor\n",
    "    print(u500Avg, v500Avg)\n",
    "    uDiff = (u6kmAvg-u500Avg)\n",
    "    vDiff = (v6kmAvg-v500Avg)\n",
    "    print(uDiff, vDiff)\n",
    "    mag = np.sqrt(uDiff**2+vDiff**2)\n",
    "    brnshr = (mag**2)*.5\n",
    "    brnshr = brnshr.to('m^2/s^2')\n",
    "    brn = cape/brnu\n",
    "    brn = brn.magnitude\n",
    "\n",
    "    return brn, brnshr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 453,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.1212028197954187 kilogram / meter ** 3\n",
      "1.1124058850537535 kilogram / meter ** 3\n",
      "1.1024174949009757 kilogram / meter ** 3\n",
      "1.0759423235482144 kilogram / meter ** 3\n",
      "1.063655836068516 kilogram / meter ** 3\n",
      "1.0517231890518561 kilogram / meter ** 3\n",
      "1.0494171181571672 kilogram / meter ** 3\n",
      "1.0216280462859253 kilogram / meter ** 3\n",
      "0.9964913714755524 kilogram / meter ** 3\n",
      "0.9686903983105212 kilogram / meter ** 3\n",
      "7.3759664942045475 5.234657939223949\n",
      "5.181513773810572 8.974645116357859\n",
      "2.1944527203939757 -3.7399871771339104\n"
     ]
    },
    {
     "ename": "AttributeError",
     "evalue": "'numpy.float64' object has no attribute 'to'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-453-182fed9617e0>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mcalc_BRN\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m5065\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0munits\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'J/kg'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mheights\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mheights\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-452-e07cb8f9e12f>\u001b[0m in \u001b[0;36mcalc_BRN\u001b[0;34m(pressure, u, v, temp, cape, heights)\u001b[0m\n\u001b[1;32m     84\u001b[0m     \u001b[0mmag\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msqrt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0muDiff\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mvDiff\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     85\u001b[0m     \u001b[0mbrnshr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mmag\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;36m.5\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 86\u001b[0;31m     \u001b[0mbrnshr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbrnshr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'm^2/s^2'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     87\u001b[0m     \u001b[0mbrn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcape\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mbrnu\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     88\u001b[0m     \u001b[0mbrn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbrn\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmagnitude\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'numpy.float64' object has no attribute 'to'"
     ]
    }
   ],
   "source": [
    "calc_BRN(p,u,v,t,5065 * units('J/kg'),heights=heights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 450,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_BRN_2(pressure, u, v, temp, cape, heights=None):\n",
    "    r\"\"\"Calculate Bulk Richardson Number and BRN shear.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    pressure : `pint.Quantity`\n",
    "                Atmospheric pressure profile\n",
    "    u        : `pint.Quantity`\n",
    "                U-component of wind.\n",
    "    v        : `pint.Quantity`\n",
    "                V-component of wind.\n",
    "    temp     : `pint.Quantity`\n",
    "                Atmospheric temperature profile (can be either temperature or virtual temperature)\n",
    "    cape     : `pint.Quantity`\n",
    "                CAPE value to use as the numerator of the BRN calculation.\n",
    "    heights  : `pint.Quantity`, optional\n",
    "                Heights in meters from sounding (not AGL)\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    brn     : `pint.Quantity`\n",
    "               Bulk Richadson Number\n",
    "    brn_shr : `pint.Quantity`\n",
    "               Bulk Richardson Shear (the denominator of the BRN)\n",
    "\n",
    "    \"\"\"\n",
    "    \n",
    "    from metpy import constants as mpconsts\n",
    "    \n",
    "    baseZ = heights[0]\n",
    "    sumZ = -baseZ\n",
    "    rho6km = 0\n",
    "    u6km = 0\n",
    "    v6km = 0\n",
    "    u500 = 0\n",
    "    v500 = 0\n",
    "\n",
    "    rho = (pressure[0].to('Pa')/(temp[0].to('K')*mpconsts.Rd)).to('kg / m^3')\n",
    "    u_weighted = u[0]*rho\n",
    "    v_weighted = v[0]*rho\n",
    "    u_weighted = u_weighted.magnitude\n",
    "    v_weighted = v_weighted.magnitude\n",
    "    ubase_weighted = u_weighted\n",
    "    vbase_weighted = v_weighted\n",
    "\n",
    "    for i, p in enumerate(pressure):\n",
    "        rho = (p.to('Pa')/(temp[i].to('K')*mpconsts.Rd)).to('kg / m^3')\n",
    "        print(rho)\n",
    "        sumZ += heights[i]\n",
    "        u_weighted = u[i]*rho\n",
    "        v_weighted = v[i]*rho\n",
    "        u_weighted = u_weighted.magnitude\n",
    "        v_weighted = v_weighted.magnitude\n",
    "        \n",
    "        if (sumZ >= 6000 * units.meter):\n",
    "            layerZ=heights[i] - heights[i-1]\n",
    "            fraction = (6000 * units.meter - layerZ)/sumZ\n",
    "            fraction = fraction.to('dimensionless')\n",
    "            u_frac = ((u_weighted - uprev) * fraction) + uprev\n",
    "            v_frac = ((v_weighted - vprev) * fraction) + vprev\n",
    "\n",
    "            u_weighted = u_frac*rho\n",
    "            v_weighted = v_frac*rho\n",
    "            u_weighted = u_weighted.magnitude\n",
    "            v_weighted = v_weighted.magnitude\n",
    "            u6km += u_weighted\n",
    "            v6km += v_weighted\n",
    "            \n",
    "            break\n",
    "\n",
    "        if (sumZ < 500 * units.meter):\n",
    "            u500 += u_weighted\n",
    "            v500 += v_weighted\n",
    "            divisor = i+1\n",
    "\n",
    "        uprev = u_weighted\n",
    "        vprev = v_weighted\n",
    "\n",
    "    u6kmAvg = u6km/(i+1)\n",
    "    v6kmAvg = v6km/(i+1)\n",
    "    print(u6kmAvg, v6kmAvg)\n",
    "    u500Avg = u500/divisor\n",
    "    v500Avg = v500/divisor\n",
    "    print(u500Avg, v500Avg)\n",
    "    ushr = (u6kmAvg-u500Avg)\n",
    "    vshr = (v6kmAvg-v500Avg)\n",
    "    print(ushr, vshr)\n",
    "    mag = np.sqrt(ushr**2+vshr**2)\n",
    "    print(mag)\n",
    "#    brnu = (mag)*.5\n",
    "    brnshr = (mag**2)*.5\n",
    "    brnshr = brnshr.to('m^2/s^2')\n",
    "    brn = cape/brnshr\n",
    "    brn = brn.magnitude\n",
    "\n",
    "    return brn, brnshr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 451,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.1212028197954187 kilogram / meter ** 3\n",
      "1.1124058850537535 kilogram / meter ** 3\n",
      "1.1024174949009757 kilogram / meter ** 3\n",
      "1.0759423235482144 kilogram / meter ** 3\n",
      "1.063655836068516 kilogram / meter ** 3\n",
      "1.0517231890518561 kilogram / meter ** 3\n",
      "1.0494171181571672 kilogram / meter ** 3\n",
      "1.0216280462859253 kilogram / meter ** 3\n",
      "0.9964913714755524 kilogram / meter ** 3\n",
      "0.9686903983105212 kilogram / meter ** 3\n",
      "1.3013121276200768 -0.4485221039445009\n",
      "5.181513773810572 8.974645116357859\n",
      "-3.8802016461904953 -9.42316722030236\n",
      "10.190782368242408\n"
     ]
    },
    {
     "ename": "AttributeError",
     "evalue": "'numpy.float64' object has no attribute 'to'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-451-72d774803723>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mcalc_BRN_2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m5065\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0munits\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'J/kg'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mheights\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mheights\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-450-7302b8521aeb>\u001b[0m in \u001b[0;36mcalc_BRN_2\u001b[0;34m(pressure, u, v, temp, cape, heights)\u001b[0m\n\u001b[1;32m     90\u001b[0m \u001b[0;31m#    brnu = (mag)*.5\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     91\u001b[0m     \u001b[0mbrnshr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mmag\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;36m.5\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 92\u001b[0;31m     \u001b[0mbrnshr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbrnshr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'm^2/s^2'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     93\u001b[0m     \u001b[0mbrn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcape\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mbrnshr\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     94\u001b[0m     \u001b[0mbrn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbrn\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmagnitude\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'numpy.float64' object has no attribute 'to'"
     ]
    }
   ],
   "source": [
    "calc_BRN_2(p,u,v,t,5065 * units('J/kg'),heights=heights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 456,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_BRN_3(pressure, u, v, cape, heights=None):\n",
    "    r\"\"\"Calculate Bulk Richardson Number and BRN shear.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    pressure : `pint.Quantity`\n",
    "                Atmospheric pressure profile\n",
    "    u        : `pint.Quantity`\n",
    "                U-component of wind.\n",
    "    v        : `pint.Quantity`\n",
    "                V-component of wind.\n",
    "    cape     : `pint.Quantity`\n",
    "                CAPE value to use as the numerator of the BRN calculation.\n",
    "    heights  : `pint.Quantity`, optional\n",
    "                Heights in meters from sounding (not AGL)\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    brn     : `pint.Quantity`\n",
    "               Bulk Richadson Number\n",
    "    brn_shr : `pint.Quantity`\n",
    "               Bulk Richardson Shear (the denominator of the BRN)\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    p_layer500, u_layer500, v_layer500 = get_layer(p, u, v, heights=None, depth=500*units.meter)\n",
    "    u_layer500 = u_layer500.to('m/s')\n",
    "    v_layer500 = v_layer500.to('m/s')\n",
    "    avgU500 = sum(u_layer500 * p_layer500) / sum(p_layer500)\n",
    "    avgV500 = sum(v_layer500 * p_layer500) / sum(p_layer500)\n",
    "    print(avgU500, avgV500)\n",
    "    \n",
    "    p_layer6km, u_layer6km, v_layer6km = get_layer(p, u, v, heights=None, depth=6000*units.meter)\n",
    "    u_layer6km = u_layer6km.to('m/s')\n",
    "    v_layer6km = v_layer6km.to('m/s')\n",
    "    avgU6km = sum(u_layer6km * p_layer6km) / sum(p_layer6km)\n",
    "    avgV6km = sum(v_layer6km * p_layer6km) / sum(p_layer6km)\n",
    "    print(avgU6km, avgV6km)\n",
    "\n",
    "    ushr = avgU6km-avgU500\n",
    "    vshr = avgV6km-avgV500\n",
    "    print(ushr, vshr)\n",
    "    \n",
    "    mag = np.sqrt(ushr**2+vshr**2)\n",
    "    print(mag)\n",
    "    \n",
    "    brnshr = (mag**2)*.5\n",
    "    brn = cape/brnshr\n",
    "    brn = brn.magnitude\n",
    "\n",
    "    return brn, brnshr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 457,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.6598553905049123 meter / second 4.577493242528277 meter / second\n",
      "8.234877431680705 meter / second -2.8785490523619663 meter / second\n",
      "5.575022041175792 meter / second -7.456042294890244 meter / second\n",
      "9.309857005496275 meter / second\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(116.87548453756854, 43.33671873139403 <Unit('meter ** 2 / second ** 2')>)"
      ]
     },
     "execution_count": 457,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calc_BRN_3(p,u,v,5065 * units('J/kg'),heights=heights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 446,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_BRN_4(pressure, spd, cape, heights=None):\n",
    "    r\"\"\"Calculate Bulk Richardson Number and BRN shear.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    pressure : `pint.Quantity`\n",
    "                Atmospheric pressure profile\n",
    "    spd        : `pint.Quantity`\n",
    "                Wind speed.\n",
    "    cape     : `pint.Quantity`\n",
    "                CAPE value to use as the numerator of the BRN calculation.\n",
    "    heights  : `pint.Quantity`, optional\n",
    "                Heights in meters from sounding (not AGL)\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    brn     : `pint.Quantity`\n",
    "               Bulk Richadson Number\n",
    "    brn_shr : `pint.Quantity`\n",
    "               Bulk Richardson Shear (the denominator of the BRN)\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    p_layer500, spd_layer500 = get_layer(p, spd, heights=None, depth=500*units.meter)\n",
    "    avg500 = sum(spd_layer500 * p_layer500) / sum(p_layer500)\n",
    "    print(avg500)\n",
    "    \n",
    "    p_layer6km, spd_layer6km = get_layer(p, spd, heights=None, depth=6000*units.meter)\n",
    "    avg6km = sum(spd_layer6km * p_layer6km) / sum(p_layer6km)\n",
    "    print(avg6km)\n",
    "\n",
    "    shr = avg6km-avg500\n",
    "    print(shr)\n",
    "\n",
    "    print(spd_layer6km)\n",
    "    \n",
    "    brnshr = (shr**2)*.5\n",
    "    brnshr = brnshr.to('m^2/s^2')\n",
    "    brn = cape/brnshr\n",
    "    brn = brn.magnitude\n",
    "\n",
    "    return brn, brnshr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 447,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.292262720773174 knot\n",
      "19.41731217708562 knot\n",
      "9.125049456312446 knot\n",
      "[7.0 9.0 12.0 12.0 10.0 8.0 8.0 9.0 12.0 15.0 18.0 21.0 21.0 22.0 22.0 24.0 24.0 23.0 22.0 22.0 22.0 23.0 23.0 30.0 30.0 31.0 31.0 32.0 33.0 33.0 34.0 34.0 34.741301661987926] knot\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(459.6867949666405, 11.01837175977084 <Unit('meter ** 2 / second ** 2')>)"
      ]
     },
     "execution_count": 447,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calc_BRN_4(p, spd, 5065 * units('J/kg'), heights=heights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 403,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_BRN_5(pressure, u, v, cape, heights=None):\n",
    "    r\"\"\"Calculate Bulk Richardson Number and BRN shear.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    pressure : `pint.Quantity`\n",
    "                Atmospheric pressure profile\n",
    "    u        : `pint.Quantity`\n",
    "                U-component of wind.\n",
    "    v        : `pint.Quantity`\n",
    "                V-component of wind.\n",
    "    cape     : `pint.Quantity`\n",
    "                CAPE value to use as the numerator of the BRN calculation.\n",
    "    heights  : `pint.Quantity`, optional\n",
    "                Heights in meters from sounding (not AGL)\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    brn     : `pint.Quantity`\n",
    "               Bulk Richadson Number\n",
    "    brn_shr : `pint.Quantity`\n",
    "               Bulk Richardson Shear (the denominator of the BRN)\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    p_layer500, u_layer500, v_layer500 = get_layer(p, u, v, heights=None, depth=500*units.meter)\n",
    "    print(u_layer500, v_layer500)\n",
    "    avgU500 = sum(u_layer500) / p_layer500.size\n",
    "    avgV500 = sum(v_layer500) / p_layer500.size\n",
    "    print(avgU500, avgV500)\n",
    "    \n",
    "    p_layer6km, u_layer6km, v_layer6km = get_layer(p, u, v, heights=None, depth=6000*units.meter)\n",
    "    print(u_layer6km, v_layer6km, p_layer6km)\n",
    "    avgU6km = sum(u_layer6km) / p_layer6km.size\n",
    "    avgV6km = sum(v_layer6km) / p_layer6km.size\n",
    "    print(avgU6km, avgV6km)\n",
    "\n",
    "    ushr = avgU6km-avgU500\n",
    "    vshr = avgV6km-avgV500\n",
    "    print(ushr, vshr)\n",
    "    \n",
    "    mag = np.sqrt(ushr**2+vshr**2)\n",
    "    print(mag)\n",
    "    \n",
    "    brnshr = (mag**2)*.5\n",
    "    brn = cape/brnshr\n",
    "    brn = brn.magnitude\n",
    "\n",
    "    return brn, brnshr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 404,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3.500000000000001 4.500000000000001 6.000000000000002 6.000000000000002 5.955203263159421] knot [6.06217782649107 7.794228634059947 10.392304845413264 10.392304845413264 10.01919862464582] knot\n",
      "5.1910406526318855 knot 8.932042955204674 knot\n",
      "[3.500000000000001 4.500000000000001 6.000000000000002 6.000000000000002 5.735764363510461 5.455986880499987 5.65685424949238 9.0 11.817693036146498 14.095389311788624 17.38666487320323 19.733545036504076 19.733545036504076 19.05255888325765 18.857680615446473 19.659649062935802 20.12809363069018 20.30779463575532 20.801408663184972 21.250368178359505 20.673237657289985 20.845079101842945 20.845079101842945 24.574561328669756 24.574561328669756 25.079526825623375 24.757700811466087 25.216344115415097 25.27946662292628 25.27946662292628 27.153607341607966 28.187277466871414 29.531685895191128] knot [6.06217782649107 7.794228634059947 10.392304845413264 10.392304845413264 8.191520442889917 5.8508296129533655 5.6568542494923815 1.7763568394002505e-15 -2.0837781320031596 -5.130302149885035 -4.658742811845365 -7.18242300983905 -7.18242300983905 -11.0 -11.33083764802119 -13.765834472425105 -13.07133684036064 -10.79784594407548 -7.162499398057439 -5.6940189922554465 -7.524443153164719 -9.720220020036091 -9.720220020036091 -17.20729309053138 -17.20729309053138 -18.22134282106666 -18.656265717713485 -19.701167210421072 -21.211991119655796 -21.211991119655796 -20.46171078716963 -19.012558718005398 -18.28146767383094] knot [995.0 984.0 972.0 939.7 925.0 911.0 908.4 877.6 850.0 818.5 790.1 763.0 762.8 735.8 732.0 700.0 696.0 686.0 667.0 659.7 612.3 568.2 568.0 501.0 500.0 496.0 494.0 492.0 487.0 486.2 477.0 467.0 461.8007636293167] hectopascal\n",
      "17.89910880841279 knot -7.96538746950643 knot\n",
      "12.708068155780905 knot -16.897430424711104 knot\n",
      "21.142803768893213 knot\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(24.759620990876353, 223.50907560496253 <Unit('knot ** 2')>)"
      ]
     },
     "execution_count": 404,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calc_BRN_5(p,u,v,5534 * units('J/kg'),heights=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 458,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_BRN_6(pressure, u, v, temp, cape, heights=None):\n",
    "    r\"\"\"Calculate Bulk Richardson Number and BRN shear.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    pressure : `pint.Quantity`\n",
    "                Atmospheric pressure profile\n",
    "    u        : `pint.Quantity`\n",
    "                U-component of wind.\n",
    "    v        : `pint.Quantity`\n",
    "                V-component of wind.\n",
    "    temp     : `pint.Quantity`\n",
    "                Atmospheric temperature profile\n",
    "    cape     : `pint.Quantity`\n",
    "                CAPE value to use as the numerator of the BRN calculation.\n",
    "    heights  : `pint.Quantity`, optional\n",
    "                Heights in meters from sounding (not AGL)\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    brn     : `pint.Quantity`\n",
    "               Bulk Richadson Number\n",
    "    brn_shr : `pint.Quantity`\n",
    "               Bulk Richardson Shear (the denominator of the BRN)\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    p_layer500, u_layer500, v_layer500, t_layer500 = get_layer(p, u, v, t, heights=None, depth=500*units.meter)\n",
    "    rho_layer500 = (p_layer500.to('Pa')/(t_layer500.to('K')*mpconsts.Rd)).to('kg/m^3')\n",
    "    print(rho_layer500)\n",
    "    \n",
    "    u_layer500 = u_layer500.to('m/s')\n",
    "    v_layer500 = v_layer500.to('m/s')\n",
    "    \n",
    "    avgU500 = sum(u_layer500 * rho_layer500) / sum(rho_layer500)\n",
    "    avgV500 = sum(v_layer500 * rho_layer500) / sum(rho_layer500)\n",
    "    print(avgU500, avgV500)\n",
    "    \n",
    "    p_layer6km, u_layer6km, v_layer6km, t_layer6km = get_layer(p, u, v, t, heights=None, depth=6000*units.meter)\n",
    "    rho_layer6km = (p_layer6km.to('Pa')/(t_layer6km.to('K')*mpconsts.Rd)).to('kg/m^3')\n",
    "    print(rho_layer6km)\n",
    "    \n",
    "    u_layer6km = u_layer6km.to('m/s')\n",
    "    v_layer6km = v_layer6km.to('m/s')\n",
    "    \n",
    "    avgU6km = sum(u_layer6km * rho_layer6km) / sum(rho_layer6km)\n",
    "    avgV6km = sum(v_layer6km * rho_layer6km) / sum(rho_layer6km)\n",
    "    print(avgU6km, avgV6km)\n",
    "\n",
    "    ushr = avgU6km-avgU500\n",
    "    vshr = avgV6km-avgV500\n",
    "    print(ushr, vshr)\n",
    "    \n",
    "    mag = np.sqrt(ushr**2+vshr**2)\n",
    "    print(mag)\n",
    "    \n",
    "    brnshr = (mag**2)*.5\n",
    "    brnshr = brnshr.to(\"m^2/s^2\")\n",
    "    brn = cape/brnshr\n",
    "    brn = brn.magnitude\n",
    "\n",
    "    return brn, brnshr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 464,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.1212028197954187 1.1124058850537535 1.1024174949009757 1.0759423235482144 1.0738480240513477] kilogram / meter ** 3\n",
      "2.66281996716048 meter / second 4.58238029054885 meter / second\n",
      "[1.1212028197954187 1.1124058850537535 1.1024174949009757 1.0759423235482144 1.063655836068516 1.0517231890518561 1.0494171181571672 1.0216280462859253 0.9964913714755524 0.9686903983105212 0.9434123310970598 0.9192458814127077 0.9190049257426126 0.8920289175479302 0.8880401711622702 0.85577526593426 0.8514827581661372 0.8374841380319189 0.8177273526188595 0.8113475273728411 0.7687938166965704 0.7286548692608193 0.7283983909541453 0.6649997777863086 0.6636724329204676 0.6568613700050064 0.6537157077292952 0.6515641008920628 0.6449425145008832 0.6441279296878698 0.6336263559163309 0.623194388962491 0.6180111600294343] kilogram / meter ** 3\n",
      "8.451420318874625 meter / second -3.1497349061140034 meter / second\n",
      "5.788600351714145 meter / second -7.732115196662853 meter / second\n",
      "9.65885601126395 meter / second\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(118.63634728770406, 46.64674972316487 <Unit('meter ** 2 / second ** 2')>)"
      ]
     },
     "execution_count": 464,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calc_BRN_6(p,u,v,t,5534 * units('J/kg'),heights=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 469,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_BRN_7(pressure, u, v, cape, heights=None):\n",
    "    r\"\"\"Calculate Bulk Richardson Number and BRN shear.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    pressure : `pint.Quantity`\n",
    "                Atmospheric pressure profile\n",
    "    u        : `pint.Quantity`\n",
    "                U-component of wind.\n",
    "    v        : `pint.Quantity`\n",
    "                V-component of wind.\n",
    "    temp     : `pint.Quantity`\n",
    "                Atmospheric temperature profile\n",
    "    cape     : `pint.Quantity`\n",
    "                CAPE value to use as the numerator of the BRN calculation.\n",
    "    heights  : `pint.Quantity`, optional\n",
    "                Heights in meters from sounding (not AGL)\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    brn     : `pint.Quantity`\n",
    "               Bulk Richadson Number\n",
    "    brn_shr : `pint.Quantity`\n",
    "               Bulk Richardson Shear (the denominator of the BRN)\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    p_layer500, u_layer500, v_layer500 = get_layer(p, u, v, heights=None, depth=500*units.meter)\n",
    "    u_layer500 = u_layer500.to('m/s')\n",
    "    v_layer500 = v_layer500.to('m/s')\n",
    "    \n",
    "    avgU500 = sum(u_layer500 * p_layer500) / sum(p_layer500)\n",
    "    avgV500 = sum(v_layer500 * p_layer500) / sum(p_layer500)\n",
    "    print(avgU500, avgV500)\n",
    "    \n",
    "    p_layer6km, u_layer6km, v_layer6km = get_layer(p, u, v, heights=None, depth=6000*units.meter)    \n",
    "    u_layer6km = u_layer6km.to('m/s')\n",
    "    v_layer6km = v_layer6km.to('m/s')\n",
    "    \n",
    "    avgU6km = sum(u_layer6km * p_layer6km) / sum(p_layer6km)\n",
    "    avgV6km = sum(v_layer6km * p_layer6km) / sum(p_layer6km)\n",
    "    print(avgU6km, avgV6km)\n",
    "\n",
    "    ushr = avgU6km-avgU500\n",
    "    vshr = avgV6km-avgV500\n",
    "    print(ushr, vshr)\n",
    "    \n",
    "    mag = np.sqrt(ushr**2+vshr**2)\n",
    "    print(mag)\n",
    "    \n",
    "    brnshr = (mag**2)*.5\n",
    "    brn = cape/brnshr\n",
    "    brn = brn.magnitude\n",
    "\n",
    "    return brn, brnshr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 471,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.6598553905049123 meter / second 4.577493242528277 meter / second\n",
      "8.234877431680705 meter / second -2.8785490523619663 meter / second\n",
      "5.575022041175792 meter / second -7.456042294890244 meter / second\n",
      "9.309857005496275 meter / second\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(127.69771597846088, 43.33671873139403 <Unit('meter ** 2 / second ** 2')>)"
      ]
     },
     "execution_count": 471,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calc_BRN_7(p,u,v,5534 * units('J/kg'),heights=None)"
   ]
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}