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

# <div style="width:1000 px">
# 
# <div style="float:right; width:98 px; height:98px;"><img src="https://pbs.twimg.com/profile_images/1187259618/unidata_logo_rgb_sm_400x400.png" alt="Unidata Logo" style="height: 98px;"></div>
# 
# <h1>Using Siphon to get NEXRAD Level 3 data from a TDS</h1>
# <h3>Unidata Python Workshop</h3>
# 
# <div style="clear:both"></div>
# </div>
# 
# <hr style="height:2px;">
# 
# <div style="float:right; width:250 px"><img src="https://upload.wikimedia.org/wikipedia/commons/4/4d/Siphoning.JPG" alt="Siphoning" style="height: 300px;"></div>
# 
# ### Objectives
# 1. Learn more about Siphon
# 2. Use the RadarServer class to retrieve radar data from a TDS
# 3. Plot this data using numpy arrays and matplotlib
# 
# In this example, we'll focus on interacting with the Radar Query Service to retrieve radar data.
# 
# **But first!**
# Bookmark these resources for when you want to use Siphon later!
# + [latest Siphon documentation](http://siphon.readthedocs.org/en/latest/)
# + [Siphon github repo](https://github.com/Unidata/siphon)
# + [TDS documentation](http://www.unidata.ucar.edu/software/thredds/v4.6/tds/TDS.html)
# 

# ## Querying the server

# First, we point at the top level of the Radar Query Service (the "Radar Server") to see what radar collections are available:

# In[1]:


from siphon.catalog import TDSCatalog
cat = TDSCatalog('http://thredds.ucar.edu/thredds/radarServer/catalog.xml')
list(cat.catalog_refs)


# Next we create an instance of the `RadarServer` object to point at one of these collections. This downloads some top level metadata and sets things up so we can easily query the server.

# In[2]:


from siphon.radarserver import RadarServer
rs = RadarServer(cat.catalog_refs['NEXRAD Level III Radar from IDD'].href)


# We can use rs.variables to see a list of radar products available to view from this access URL.

# In[3]:


print(sorted(rs.variables))


# If you're not a NEXRAD radar expert, there is more information available within the metadata downloaded from the server. (**NOTE:** Only the codes above are valid for queries.)

# In[4]:


sorted(rs.metadata['variables'])


# We can also see a list of the stations. Each station has associated location information.

# In[5]:


print(sorted(rs.stations))


# In[6]:


rs.stations['TLX']


# Next, we'll create a new query object to help request the data. Using the chaining methods, let's ask for reflectivity data at the lowest tilt (NOQ) from radar TLX (Oklahoma City) for the current time. We see that when the query is represented as a string, it shows the encoded URL.

# In[7]:


from datetime import datetime
query = rs.query()
query.stations('TLX').time(datetime.utcnow()).variables('N0Q')


# The query also supports time range queries, queries for closest to a lon/lat point, or getting all radars within a lon/lat box.
# 
# We can use the RadarServer instance to check our query, to make sure we have required parameters and that we have chosen valid station(s) and variable(s).

# In[8]:


rs.validate_query(query)


# Make the request, which returns an instance of TDSCatalog. This handles parsing the catalog

# In[9]:


catalog = rs.get_catalog(query)


# We can look at the datasets on the catalog to see what data we found by the query. We find one NIDS file in the return

# In[10]:


catalog.datasets


# ### Exercise: Querying the radar server
# 
# We'll work through doing some more queries on the radar server. Some useful links:
# - RadarQuery [documentation](https://siphon.readthedocs.org/en/latest/api/radarserver.html#siphon.radarserver.RadarQuery)
# - Documentation on Python's [datetime.timedelta](https://docs.python.org/3.5/library/datetime.html#timedelta-objects)
# 
# See if you can write Python code for the following queries:
# 
# Get ZDR (differential reflectivity) for 3 days ago from the radar nearest to Hays, KS (lon -99.324403, lat 38.874929). **No map necessary!**

# In[ ]:





# Get base reflectivity for the last two hours from all of the radars in Wyoming (call it the bounding box with lower left corner 41.008717, -111.056360 and upper right corner 44.981008, -104.042719)

# In[ ]:





# ## Pulling out the data

# We can pull that dataset out of the dictionary and look at the available access URLs. We see URLs for OPeNDAP, CDMRemote, and HTTPServer (direct download).

# In[11]:


ds = list(catalog.datasets.values())[0]
ds.access_urls


# We'll use the CDMRemote reader in Siphon and pass it the appropriate access URL. (This will all behave identically to using the 'OPENDAP' access, if we replace the `Dataset` from Siphon with that from `netCDF4`).

# In[12]:


from siphon.cdmr import Dataset
data = Dataset(ds.access_urls['CdmRemote'])


# The CDMRemote reader provides an interface that is almost identical to the usual python NetCDF interface. 

# In[13]:


list(data.variables)


# We pull out the variables we need for azimuth and range, as well as the data itself.

# In[14]:


rng = data.variables['gate'][:]
az = data.variables['azimuth'][:]
ref = data.variables['BaseReflectivityDR'][:]


# Then convert the polar coordinates to Cartesian using numpy

# In[15]:


import numpy as np
x = rng * np.sin(np.deg2rad(az))[:, None]
y = rng * np.cos(np.deg2rad(az))[:, None]
ref = np.ma.array(ref, mask=np.isnan(ref))


# Finally, we plot them up using matplotlib and cartopy.

# In[16]:


get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt
import cartopy
from metpy.plots import ctables  # For NWS colortable

# Create projection centered on the radar. This allows us to use x
# and y relative to the radar.
proj = cartopy.crs.LambertConformal(central_longitude=data.RadarLongitude,
                                    central_latitude=data.RadarLatitude)

# New figure with specified projection
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(1, 1, 1, projection=proj)

# Grab state borders
state_borders = cartopy.feature.NaturalEarthFeature(
    category='cultural', name='admin_1_states_provinces_lakes',
    scale='50m', facecolor='none')
ax.add_feature(state_borders, edgecolor='black', linewidth=2, zorder=2)

# Counties
counties = cartopy.io.shapereader.Reader('../../data/counties.shp')
ax.add_geometries(counties.geometries(), cartopy.crs.PlateCarree(),
                  facecolor='None', edgecolor='grey', zorder=1)

# Set limits in lat/lon space
ax.set_extent([data.RadarLongitude - 2.5, data.RadarLongitude + 2.5,
               data.RadarLatitude - 2.5, data.RadarLatitude + 2.5])

# Get the NWS typical reflectivity color table, along with an appropriate norm that
# starts at 5 dBz and has steps in 5 dBz increments
norm, cmap = ctables.registry.get_with_steps('NWSReflectivity', 5, 5)
mesh = ax.pcolormesh(x, y, ref, cmap=cmap, norm=norm, zorder=0)


# ### Exercise: Your turn to plot!
# 
# Try making your own plot of radar data. Various options here, but this is pretty open-ended. Some options to inspire you:
# - Try working with Level II Data (All variables are three-dimensional and in a single file)
# - Try plotting Storm Total Precipitation or Digital Accumulation for an area that recently had a bunch of rain
# - Maybe plot Doppler velocity or some dual-pol data (ZDR or correlation coefficient) for recent severe weather. Bonus points for making a multi-panel plot.

# In[ ]: