bryan_bulk_shear_calc_SkewT.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate Bulk Wind Shear for JSON Sounding Files\n",
    "* Have specific data structure, please see https://github.com/bryanguarente/d3-skewt-and-indices/blob/master/data_DNR.json "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Problem with Data:\n",
    "* Cointains wind speed and direction, but MetPy's ```metpy.calc.bulk_shear``` requires u and v-components \n",
    "\n",
    "Solution:\n",
    "* Use MetPy's ```metpy.calc.wind_components``` to get u and v from speed and direction.\n",
    "\n",
    "Issues:\n",
    "* Potential for incomplete data ie \n",
    "    * empty entries -> causes calculation errors\n",
    "    * others??"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 236,
   "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": "markdown",
   "metadata": {},
   "source": [
    "Copied data from https://github.com/bryanguarente/d3-skewt-and-indices/blob/master/data_DNR.json to file locally called ```dnr.json```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "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": 153,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_json(\"data_IAD.json\")\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['pres', 'hght', 'tmpc', 'dwpt', 'wdir', 'wspd'], dtype='object')"
      ]
     },
     "execution_count": 154,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "pandas.core.series.Series"
      ]
     },
     "execution_count": 155,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(df['pres'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "list"
      ]
     },
     "execution_count": 156,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(df['pres'][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "list"
      ]
     },
     "execution_count": 157,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(df['pres'][0][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "int"
      ]
     },
     "execution_count": 158,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(df['pres'][0][0][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 282,
   "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": 278,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "99\n",
      "99\n",
      "99\n",
      "99\n",
      "99\n"
     ]
    }
   ],
   "source": [
    "print(p.size)\n",
    "print(heights.size)\n",
    "print(dir.size)\n",
    "print(t.size)\n",
    "print(td.size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "487.0 hectopascal\n",
      "6004 meter\n"
     ]
    }
   ],
   "source": [
    "print(p[28])\n",
    "print(heights[28])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "pint.quantity.build_quantity_class.<locals>.Quantity"
      ]
     },
     "execution_count": 166,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Need to address any entries that are empty strings and not integers or floats\n",
    "* in the ```data_DNR.json``` example, there are some in wind speed and wind direction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[353, 20]"
      ]
     },
     "execution_count": 167,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df[\"wdir\"].values[0][0][-2:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "int"
      ]
     },
     "execution_count": 168,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(df[\"wdir\"].values[0][0][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "int"
      ]
     },
     "execution_count": 169,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(df[\"wdir\"].values[0][0][-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pandas ```replace``` should work, but isn't in this case\n",
    "* might be because the data values in each column are in lists to begin with..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.replace('', 0, regex=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Didn't seem to work..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[353, 20]"
      ]
     },
     "execution_count": 171,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df[\"wdir\"].values[0][0][-2:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get it done the old fashion way with direct element replacement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [],
   "source": [
    "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 # can choose anything, ie np.nan, -9999, etc. \n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run any/all columns through to capture all instances of empty string"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in df.columns:\n",
    "    replace_empty_str(i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Double check that the strings got replaced by 0 (or whatever else is desired)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([353, 20], [14, 17])"
      ]
     },
     "execution_count": 174,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df[\"wdir\"][0][0][-2:],df[\"wspd\"][0][0][-2:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Change speed and direction to u and v components with MetPy\n",
    "* ```metpy.calc.bulk_shear``` requires u and v components\n",
    "* Use: ```metpy.calc.wind_components``` (requires speed and direction and produces u and v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "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": [
    "### Get 6k bulk shear u and v components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_shear_6km,v_shear_6km = mpcalc.bulk_shear(p,u,v,depth=6000 * units.meter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "u-6km shear: 26.0 knot\n",
      "v-6km shear: -24.0 knot\n"
     ]
    }
   ],
   "source": [
    "print(f\"u-6km shear: {np.round(u_shear_6km)}\\nv-6km shear: {np.round(v_shear_6km)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Grab speed and direction if desired fo r 6km bulk shear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [],
   "source": [
    "bulk_shear_6km = np.round(np.sqrt(u_shear_6km**2 + v_shear_6km**2))\n",
    "bulk_shear_6km_dir = np.round(mpcalc.wind_direction(u_shear_6km,v_shear_6km))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6km shear speed: 36.0 knot\n",
      "6km shear direction: 313.0 degree\n"
     ]
    }
   ],
   "source": [
    "print(f\"6km shear speed: {bulk_shear_6km}\\n6km shear direction: {bulk_shear_6km_dir}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quick function to get bulk shear at desired layer (meters) for JSON files\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 273,
   "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",
    "    print(f\"\\nSounding 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",
    "    \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,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": 200,
   "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": 318,
   "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": 324,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([487]) <Unit('hectopascal')>, 6000 <Unit('meter')>)"
      ]
     },
     "execution_count": 324,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_get_bound_pressure_height(p, 6000*units.meter, heights=heights, interpolate=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 330,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "461.8007636293167 hectopascal\n",
      "6362.090869242044 meter\n",
      "29.531685895191128 knot\n",
      "-18.28146767383094 knot\n",
      "301.77609501403623 degree\n",
      "\n",
      "995.0 hectopascal\n",
      "46.0 meter\n",
      "3.500000000000001 knot\n",
      "6.06217782649107 knot\n",
      "210.0 degree\n",
      "\n",
      "0.2870579959589441 joule / gram / kelvin\n",
      "9.80665 meter / second ** 2\n",
      "461.8007636293167 hectopascal 6152.959968043457 meter\n",
      "0.15295996804345655 kilometer\n",
      "471.614328198986 hectopascal\n",
      "0.2870579959589441 joule / gram / kelvin\n",
      "9.80665 meter / second ** 2\n",
      "5.999999999999998 kilometer\n",
      "\n",
      "26.031685895191128 knot\n",
      "-24.34364550032201 knot\n",
      "35.6407315692487 knot\n"
     ]
    }
   ],
   "source": [
    "#_, u_layer, v_layer = get_layer(p, u, v, \n",
    "#    heights=heights, depth=6000*units.meter) # Using heights=heights makes this calculation wrong.\n",
    "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=None, depth=6000*units.meter)\n",
    "\n",
    "#print(p)\n",
    "#print(p_layer)\n",
    "#print(hghts_layer)\n",
    "#print(dir_layer)\n",
    "#print(spd_layer)\n",
    "#print(u_layer)\n",
    "#print(v_layer)\n",
    "\n",
    "print(p_layer[-1])\n",
    "print(hghts_layer[-1])\n",
    "print(u_layer[-1])\n",
    "print(v_layer[-1])\n",
    "print(dir_layer[-1])\n",
    "print()\n",
    "print(p_layer[0])\n",
    "print(hghts_layer[0])\n",
    "print(u_layer[0])\n",
    "print(v_layer[0])\n",
    "print(dir_layer[0])\n",
    "print()\n",
    "ushr = u_layer[-1] - u_layer[0]\n",
    "vshr = v_layer[-1] - v_layer[0]\n",
    "\n",
    "baseHght = pressure_to_height_std(p[0])\n",
    "pressure_bound, height_bound = _get_bound_pressure_height(p, 6000*units.meter+baseHght)\n",
    "\n",
    "print(pressure_bound, height_bound)\n",
    "\n",
    "\n",
    "print(baseHght)\n",
    "\n",
    "stdPres = height_to_pressure_std(6000*units.meter)\n",
    "print(stdPres)\n",
    "stdHght = pressure_to_height_std(stdPres)\n",
    "print(stdHght)\n",
    "print()\n",
    "\n",
    "#for i in p_layer:\n",
    "#    Z = pressure_to_height_std(i)\n",
    "#    print(Z*1000*units.meter/units.kilometer)\n",
    "    \n",
    "#print()\n",
    "\n",
    "#for i in p_layer:\n",
    "#    print(i)\n",
    "\n",
    "print(ushr)\n",
    "print(vshr)\n",
    "print(np.sqrt(ushr**2 + vshr**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 329,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Sounding Location: data_IAD.json\n",
      "Desired layer: 6.0km\n",
      "---------------------------------\n",
      "u-bulk shear: 26.031685895191128 knot\n",
      "v-bulk shear: -24.34364550032201 knot\n",
      "bulk shear speed: 35.6407315692487 knot\n",
      "bulk shear direction: 313.08077859722147 degree\n"
     ]
    }
   ],
   "source": [
    "u_bulk_shear, v_bulk_shear, bulk_shear_speed, bulk_shear_dir = get_bulk_shear(\"data_IAD.json\", 6000)"
   ]
  },
  {
   "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
}