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

# In[1]:


# Import essential packages and setup inline plotting.
import matplotlib.pyplot as plt
from metpy.units import units
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('u-component_of_wind_height_above_ground', 'v-component_of_wind_height_above_ground', 'Wind_speed_gust_surface', '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']
wind_gust = nc['Wind_speed_gust_surface']
height_temp = nc['height_above_ground4']
height_wind = nc['height_above_ground3']


# In[7]:


print(temp.units,u_wind.units,v_wind.units,wind_gust.units,height_temp.units,height_wind.units)


# In[8]:


temp, u_wind, v_wind, wind_gust, height_temp, height_wind


# In[9]:


height_temp[:], height_wind[:]


# In[47]:


u_wind_surf = u_wind[:,0,:,:] * units('m/s').to('mph')
v_wind_surf = v_wind[:,0,:,:] * units('m/s').to('mph')
t_surface = temp[:,0,:,:] * units('K').to('degF')
surf_wind = np.sqrt((u_wind_surf[:,:,:]*u_wind_surf[:,:,:])+(v_wind_surf[:,:,:]*v_wind_surf[:,:,:]))


# In[48]:


u_wind_surf.shape, v_wind_surf.shape, surf_wind.shape, t_surface.shape


# In[49]:


import metpy.calc
surf_wind_metpy = metpy.calc.get_wind_speed(u_wind_surf,v_wind_surf)


# In[50]:


surf_wind_metpy.shape


# In[51]:


wind_chill_metpy = metpy.calc.windchill(t_surface,surf_wind,mask_undefined=False)


# In[61]:


wind_chill_calc = 35.74* (0.6215*t_surface[:]) - (35.75*(surf_wind[:])**0.16) + (0.4275*t_surface[:]*(surf_wind[:])**0.16)


# In[42]:


import cartopy
import cartopy.crs as ccrs


# In[43]:


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[44]:


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


# In[45]:


# 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. 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[92].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('Wind Chill')



# In[ ]:





# In[ ]:





# In[ ]: