#!/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
import metpy as mp
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]:


ncss.variables


# In[4]:


query.all_times().accept('netcdf4').variables('u-component_of_wind_height_above_ground', 'v-component_of_wind_height_above_ground','Temperature_height_above_ground')
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]:


temp = nc['Temperature_height_above_ground']
u_wind = nc['u-component_of_wind_height_above_ground']
v_wind = nc['v-component_of_wind_height_above_ground']
height_temp = nc['height_above_ground2']
height_wind = nc['height_above_ground4']


# In[34]:


t_surface = temp[:,0,:,:]
u_wind_surf = u_wind[:,0,:,:]
v_wind_surf = v_wind[:,0,:,:]
t_surface = units.Quantity(t_surface, 'K')
t_surface_F = t_surface.to('degF')
u_wind_knots = units.Quantity(u_wind_surf, 'm/s').to('knots')
v_wind_knots = units.Quantity(v_wind_surf, 'm/s').to('knots')


# In[8]:


import cartopy
import cartopy.crs as ccrs


# In[9]:


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

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


# In[45]:


# Open and read netCDF variables
lat_var = nc['lat']
lon_var = nc['lon']
lat = lat_var[:].squeeze()
lon = lon_var[:].squeeze()
lon_2d, lat_2d = np.meshgrid(lon, lat)
lon_2d[lon_2d > 180] = lon_2d[lon_2d > 180] - 360
time_sel = 8
temp_plot = t_surface_F[time_sel,:,:]
u_plot = u_wind_knots[time_sel,:,:]
v_plot = v_wind_knots[time_sel,:,:]
lon_2d.shape, lat_2d.shape, u_plot.shape, v_plot.shape


# In[46]:


# Set up a 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. Just work, would you!',
            fontsize=12, ha='center')

# define color map
cmap=plt.get_cmap("turbo")
#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')
ax.add_feature(cartopy.feature.BORDERS,edgecolor='black')
ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.coastlines()
ax.gridlines()
ax.set_title(time[time_sel].isoformat());

# Color-filled contour plot
cs = ax.contourf(lon_2d, lat_2d, temp_plot, 10, cmap=cmap,transform=proj)

# Color bar
cbar = plt.colorbar(cs)
cbar.set_label('$\circ$ F')

ax.barbs(lon_2d,lat_2d, u_plot, v_plot, length=6, pivot='middle', color='black', regrid_shape = 10, transform=proj)

plt.show()



# In[ ]:





# In[ ]:





# In[ ]: