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

# 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.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]:


ncss.variables


# In[4]:


query.all_times().accept('netcdf4').variables('Snow_depth_surface')
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]:


nc


# In[6]:


var='Snow_depth_surface'
ncvar = nc[var]
ncvar


# In[7]:


import cartopy
import cartopy.crs as ccrs


# In[8]:


from netCDF4 import num2date
# find the correct time dimension name
for d in ncvar.dimensions:
    if "time" in d: 
        timevar = d

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


# In[9]:


# Open and read netCDF variables
lat = nc['lat']
lon = nc['lon']
data = ncvar[60, :, :]


# In[10]:


# Set up a standard lat/lon projection
proj = ccrs.PlateCarree()

# Construct figure
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(1, 1, 1, projection=proj)

# Setting the plot size and text
plt.figtext(.3, .15,
            'Figure 1. This is not what I am concerned about',
            fontsize=12, ha='center')

# define color map
cmap=plt.get_cmap("Blues")
#cmap = plt.cm.RdBu_r

# Nice high-level, human-readable abstractions for dealing with maps.
# add some common geographic features

import cartopy.feature as cfeature
states_provinces = cfeature.NaturalEarthFeature(
        category='cultural',
        name='admin_1_states_provinces_lines',
        scale='50m',
        facecolor='none')
ax.add_feature(states_provinces, edgecolor='black', zorder=2)
ax.add_feature(cartopy.feature.BORDERS,edgecolor='black',zorder=2)
ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.coastlines(zorder=3)
ax.gridlines()
ax.set_title(time[60].isoformat());

# Color-filled contour plot
cs = ax.contourf(lon[:], lat[:], data[:], 10, cmap=cmap,transform=ccrs.PlateCarree(), zorder=2)

# Color bar
cbar = plt.colorbar(cs)
cbar.set_label('Depth (m)')