#!/usr/bin/env python
# coding: utf-8

# # Instructor Help (IH) - Units

# More on units can be found at the MetPy help page at: http://metpy.readthedocs.io/en/latest/ under 'Unit Support'.

# In[1]:


# Import essential packages and setup inline plotting.
import matplotlib.pyplot as plt
import numpy as np
from metpy.units import units
get_ipython().run_line_magic('matplotlib', 'inline')


# In[2]:


# Start setting up thredds (Thematic Real-time Environmental Distributed Data Services) access 
# using siphon ncss (netCDF subset service) method.

from siphon.catalog import get_latest_access_url
gfs_catalog = "http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/Global_0p25deg/catalog.html"
latest_gfs_ncss = get_latest_access_url(gfs_catalog, "NetcdfSubset")

# Set up access via NCSS
from siphon.ncss import NCSS
ncss = NCSS(latest_gfs_ncss)

# Create a query to ask for all times in netcdf4 format for
# the Temperature_surface variable, with a bounding box centered
# on lat,lon with a height and width speficied in degrees.
query = ncss.query()
lat = 42
lon = -92.5
width = 25
height = 20


# In[3]:


# Let's take a look at the available variables.
ncss.variables


# In[4]:


# Build our query time and domain.
query.all_times().accept('netcdf4').variables('Pressure_reduced_to_MSL_msl')
query.lonlat_box(north=lat+height/2., south=lat - height/2., east=lon + width/2., west=lon - width/2.)

# Get the data!
nc = ncss.get_data(query)


# In[5]:


# Let's see what this netCDF variable looks like.
nc


# In[6]:


# Grab a variable.  I am going to grad mean sea-level pressure and mess around with the units.
var='Pressure_reduced_to_MSL_msl'
ncvar = nc[var]
ncvar


# In[7]:


# What I didn't come up with in class (embarrassing!!!!) is that all these varialbes already have units.  Of course they do, 
# they are netCDF objects.  Let's see what they are.

print(ncvar.units)


# In[8]:


# What do they look like?

ncvar[:]


# In[9]:


# Pascals are great, but let's change them to something more useful, like millibars.

ncvar = (ncvar[:] * units('Pa')).to('mbar')
print(ncvar.units)


# In[10]:


# Let's take a look at the units again.  They should be in millibars.

ncvar[:]