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

# # Instructor Help (IH) - Time Series

# In[1]:


# Import essential packages and setup inline plotting.
import matplotlib.pyplot as plt
import numpy as np
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.xml"
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]:


ncss.variables


# In[4]:


query.all_times().accept('netcdf4').variables('Temperature_height_above_ground')
#query.lonlat_box(north=lat+height/2., south=lat - height/2., east=lon + width/2., west=lon - width/2.)
query.lonlat_point(lon,lat)
# get the data!
nc = ncss.get_data(query)


# In[5]:


nc


# In[6]:


import numpy
from metpy.units import units
from netCDF4 import num2date

var='Temperature_height_above_ground'
ncvar = nc[var]

# ncvar currently is a two-dimensional array with one of the dimensions of length 1.  
# That is a waste and a minor headache to keep remembering and working with.  The squeeze method
# remove the dimension of length 1 making a simple one-dimensional array along the variable.

ncvar = ncvar[:].squeeze()

# This picks out the two-meter temperature.

ncvar = ncvar[:,0]

# This section gets the time variable and converts the units to something that matplotlib 
# recognizes and will plot.  num2date takes the time variable and the units and does the
# conversion. Notice the time variable is also two-dimensional with the first dimension of
# length one (one point, one station, etc.  Let's also squeeze it to get the dimensions of our
# data and time varialbes to match.

timevar='time'
time = num2date(nc[timevar][:],nc[timevar].units)
time = time[:].squeeze()


# In[7]:


# Now they both have the same shape......
print(ncvar.shape)
print(time.shape)


# In[8]:


# This hack is needed due to a known bug.  Hopefully it will be resolved, but for now, it is needed (11/28/2018)

import datetime as dt
newtimes = []
for tm in time:
    newtimes.append(dt.datetime(tm.year, tm.month, tm.day, tm.hour, tm.minute))


# In[9]:


#....and we can plot them. 
plt.plot(newtimes,ncvar)


# In[10]:


# Let's put a little bit of effort into making the plot look better.  Of course, to make this even better, 
# we would incorporate MetPy units and change the units to something we better understand.

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(1, 1, 1)
ax.grid()
ax.set_title('2-m Air Temperature - Ames, IA - Latest GFS Run', fontdict={'size':22})
ax.set_ylabel('K', fontsize=15)
ax.set_xlabel('Date', fontsize=15)
fig.autofmt_xdate()
plt.plot(newtimes,ncvar)


# In[11]:


#Of course, to make this even better, we would incorporate MetPy units and change the units 
#to something we better understand.  This has been somewhat complicated by the fact that numpy
#silently drops the units attribute on the masked array.  We can fix this by dropping the mask
#and making it a simple numpy array.

ncvar=(np.array(ncvar[:])*units('degK')).to('degF')
print(ncvar.units)


# In[12]:


# Let's plot it again.......

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(1, 1, 1)
ax.grid()
ax.set_title('2-m Air Temperature - Ames, IA - Latest GFS Run', fontdict={'size':22})
ax.set_ylabel('$^\circ$ F', fontsize=15)
ax.set_xlabel('Date', fontsize=15)
fig.autofmt_xdate()
plt.plot(newtimes,ncvar)


# In[ ]: