talks/MetPy Exercise.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from datetime import datetime, timedelta\n",
"from siphon.catalog import TDSCatalog\n",
"from siphon.ncss import NCSS\n",
"import numpy as np\n",
"import metpy.calc as mpcalc\n",
"from metpy.units import units, concatenate\n",
"from metpy.plots import SkewT"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Just some helper code to make things easier. `metpy_units_handler` plugins into siphon to automatically add units to variables. `post_process_data` is used to clean up some oddities from the NCSS point feature collection."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"units.define('degrees_north = 1 degree')\n",
"units.define('degrees_east = 1 degree')\n",
"unit_remap = dict(inches='inHg', Celsius='celsius')\n",
"def metpy_units_handler(vals, unit):\n",
" arr = np.array(vals)\n",
" if unit:\n",
" unit = unit_remap.get(unit, unit)\n",
" arr = arr * units(unit)\n",
" return arr\n",
"\n",
"# Fix dates and sorting\n",
"def sort_list(list1, list2):\n",
" return [l1 for (l1, l2) in sorted(zip(list1, list2), key=lambda i: i[1])]\n",
"\n",
"def post_process_data(data):\n",
" data['time'] = [datetime.strptime(d.decode('ascii'), '%Y-%m-%d %H:%M:%SZ') for d in data['time']]\n",
" ret = dict()\n",
" for key,val in data.items():\n",
" try:\n",
" val = units.Quantity(sort_list(val.magnitude.tolist(), data['time']), val.units)\n",
" except AttributeError:\n",
" val = sort_list(val, data['time'])\n",
" ret[key] = val\n",
" return ret"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# METAR Meteogram"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First we need to grab the catalog for the METAR feature collection data from http://thredds.ucar.edu/thredds/catalog.html"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"cat = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/nws/metar/ncdecoded/catalog.xml?dataset=nws/metar/ncdecoded/Metar_Station_Data_fc.cdmr')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Set up NCSS access to the dataset"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ds = list(cat.datasets.values())[0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"ncss = NCSS(ds.access_urls['NetcdfSubset'])\n",
"ncss.unit_handler = metpy_units_handler"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What variables do we have available?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ncss.variables"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create a query for the last 7 days of data for a specific lon/lat point. We should ask for: air temperature, dewpoint temperature, wind speed, and wind direction."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"now = datetime.utcnow()\n",
"query = ncss.query().accept('csv')\n",
"query.lonlat_point(-97, 35.25).time_range(now - timedelta(days=7), now)\n",
"query.variables('air_temperature', 'dew_point_temperature', 'wind_speed', 'wind_from_direction')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Get the data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"scrolled": true
},
"outputs": [],
"source": [
"data = ncss.get_data(query)\n",
"data = post_process_data(data) # Fixes for NCSS point"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Heat Index"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, we need relative humidity:\n",
" $$RH = e / e_s$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"e = mpcalc.saturation_vapor_pressure(data['dew_point_temperature'])\n",
"e_s = mpcalc.saturation_vapor_pressure(data['air_temperature'])\n",
"rh = e / e_s"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculate heat index:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# RH should be [0, 100]\n",
"hi = mpcalc.heat_index(data['air_temperature'], rh * 100)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot the temperature, dewpoint, and heat index. Bonus points to also plot wind speed and direction."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"times = data['time']\n",
"fig, axes = plt.subplots(2, 1, figsize=(9, 9))\n",
"axes[0].plot(times, data['air_temperature'].to('degF'), 'r', linewidth=2)\n",
"axes[0].plot(times, data['dew_point_temperature'].to('degF'), 'g', linewidth=2)\n",
"axes[0].plot(times, hi, color='darkred', linestyle='--', linewidth=2)\n",
"axes[0].grid(True)\n",
"axes[1].plot(times, data['wind_speed'].to('mph'), 'b')\n",
"twin = plt.twinx(axes[1])\n",
"twin.plot(times, data['wind_from_direction'], 'kx')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#Sounding"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First grab the catalog for the Best dataset from the GSD HRRR from http://thredds.ucar.edu/thredds/catalog.html"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"cat = TDSCatalog('http://thredds-jumbo.unidata.ucar.edu/thredds/catalog/grib/HRRR/CONUS_3km/wrfprs/catalog.xml?dataset=grib/HRRR/CONUS_3km/wrfprs/Best')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Set up NCSS access to the dataset"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"best_ds = list(cat.datasets.values())[0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ncss = NCSS(best_ds.access_urls['NetcdfSubset'])\n",
"ncss.unit_handler = metpy_units_handler"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What variables do we have?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ncss.variables"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Set up a query for the most recent set of data from a point. We should request temperature, dewpoint, and U and V."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"query = ncss.query().accept('csv')\n",
"query.lonlat_point(-105, 40).time(datetime.utcnow())\n",
"query.variables('Temperature_isobaric', 'Dewpoint_temperature_isobaric',\n",
" 'u-component_of_wind_isobaric', 'v-component_of_wind_isobaric')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Get the data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"data = ncss.get_data(query)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"T = data['Temperature_isobaric'].to('degC')\n",
"Td = data['Dewpoint_temperature_isobaric'].to('degC')\n",
"p = data['vertCoord'].to('mbar')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot a sounding of the data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"scrolled": false
},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(9, 9))\n",
"skew = SkewT(fig=fig)\n",
"skew.plot(p, T, 'r')\n",
"skew.plot(p, Td, 'g')\n",
"skew.ax.set_ylim(1050, 100)\n",
"skew.plot_mixing_lines()\n",
"skew.plot_dry_adiabats()\n",
"skew.plot_moist_adiabats()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Also calculate the parcel profile and add that to the plot"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"prof = mpcalc.parcel_profile(p[::-1], T[-1], Td[-1])\n",
"skew.plot(p[::-1], prof.to('degC'), 'k', linewidth=2)\n",
"fig"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's also plot the location of the LCL and the 0 isotherm as well:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"lcl = mpcalc.lcl(p[-1], T[-1], Td[-1])\n",
"lcl_temp = mpcalc.dry_lapse(concatenate((p[-1], lcl)), T[-1])[-1].to('degC')\n",
"skew.plot(lcl, lcl_temp, 'bo')\n",
"skew.ax.axvline(0, color='blue', linestyle='--', linewidth=2)\n",
"fig"
]
}
],
"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.4.3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}