# Level 3 example with multiple products
import numpy as np
import matplotlib.pyplot as plt
from numpy import ma

from metpy.cbook import get_test_data
from metpy.io.nexrad import Level3File
from metpy.plots import ctables

# Helper code for making sense of these products. This is hidden from the slideshow
# and eventually, in some form, will make its way into MetPy proper.
def print_tab_pages(prod):
    print(('\n' + '-'*80 + '\n').join(prod.tab_pages))

def print_graph_pages(prod):
    colors = {0:'white', 3:'red', 4:'cyan'}
    for page in prod.graph_pages:
        fig, ax = plt.subplots(1, 1, figsize=(10,10))
        ax.axesPatch.set_facecolor('black')
        for line in page:
            if 'color' in line:
                c = colors[line['color']]
                if 'text' in line:
                    ax.text(line['x'], line['y'], line['text'], color=c,
                            transform=ax.transData, verticalalignment='top',
                            horizontalalignment='left', fontdict={'family':'monospace'},
                            fontsize=8)
                else:
                    vecs = np.array(line['vectors'])
                    ax.plot(vecs[:, ::2], vecs[:, 1::2], color=c)
        ax.set_xlim(0, 639)
        ax.set_ylim(511, 0)
        ax.set_aspect('equal', 'box')
        ax.xaxis.set_major_formatter(plt.NullFormatter())
        ax.xaxis.set_major_locator(plt.NullLocator())
        ax.yaxis.set_major_formatter(plt.NullFormatter())
        ax.yaxis.set_major_locator(plt.NullLocator())
        for s in ax.spines: ax.spines[s].set_color('none')

def plot_prod(prod, cmap, norm, ax=None):
    if ax is None:
        ax = plt.gca()

    data_block = prod.sym_block[0][0]
    data = np.array(data_block['data'])
    data = prod.map_data(data)
    data = np.ma.array(data, mask=np.isnan(data))
    if 'start_az' in data_block:
        az = np.array(data_block['start_az'] + [data_block['end_az'][-1]])
        rng = np.linspace(0, prod.max_range, data.shape[-1] + 1)
        x = rng * np.sin(np.deg2rad(az[:, None]))
        y = rng * np.cos(np.deg2rad(az[:, None]))
    else:
        x = np.linspace(-prod.max_range, prod.max_range, data.shape[1] + 1)
        y = np.linspace(-prod.max_range, prod.max_range, data.shape[0] + 1)
        data = data[::-1]
    pc = ax.pcolormesh(x, y, data, cmap=cmap, norm=norm)
    plt.colorbar(pc, extend='both')
    ax.set_aspect('equal', 'datalim')
    ax.set_xlim(-100, 100)
    ax.set_ylim(-100, 100)
    return pc, data

def plot_points(prod, ax=None):
    if ax is None:
        ax = plt.gca()

    data_block = prod.sym_block[0]
    styles = {'MDA': dict(marker='o', markerfacecolor='None', markeredgewidth=2, size='radius'),
              'MDA (Elev.)': dict(marker='s', markerfacecolor='None', markeredgewidth=2, size='radius'),
              'TVS': dict(marker='v', markerfacecolor='red', markersize=10),
              'Storm ID': dict(text='id'),
              'HDA': dict(marker='o', markersize=10, markerfacecolor='blue', alpha=0.5)}
    artists = []
    for point in data_block:
        if 'type' in point:
            info = styles.get(point['type'], {}).copy()
            x,y = point['x'], point['y']
            text_key = info.pop('text', None)
            if text_key:
                artists.append(ax.text(x, y, point[text_key], transform=ax.transData, clip_box=ax.bbox, **info))
                artists[-1].set_clip_on(True)
            else:
                size_key = info.pop('size', None)
                if size_key:
                    info['markersize'] = np.pi * point[size_key]**2
                artists.append(ax.plot(x, y, **info))

def plot_tracks(prod, ax=None):
    if ax is None:
        ax = plt.gca()
        
    data_block = prod.sym_block[0]

    for track in data_block:
        if 'marker' in track:
            pass
        if 'track' in track:
            x,y = np.array(track['track']).T
            ax.plot(x, y, color='k')

# Read in a bunch of NIDS products
tvs = Level3File(get_test_data('nids/KOUN_SDUS64_NTVTLX_201305202016'))
nmd = Level3File(get_test_data('nids/KOUN_SDUS34_NMDTLX_201305202016'))
nhi = Level3File(get_test_data('nids/KOUN_SDUS64_NHITLX_201305202016'))
n0q = Level3File(get_test_data('nids/KOUN_SDUS54_N0QTLX_201305202016'))
nst = Level3File(get_test_data('nids/KOUN_SDUS34_NSTTLX_201305202016'))

# What happens when we print one out
tvs

# Can print tabular (ASCII) information in the product
print_tab_pages(tvs)

"""
Explanation: <h1 align="center">MetPy: An Open Source Python Toolkit for Meteorology<h1>
<h3 align="center">Ryan May<br>UCAR/Unidata<br>22 June 2015<br></h3>

## Who am I?

- Software Engineer at Unidata
- Been here for a year and a half
- Work on THREDDS
- Pythonista
- I like long walks on the beach and watching Longhorns lose football games
- Twitter/GitHub: @dopplershift

#What is MetPy?

- Python toolkit for Meteorology
    - Reading in weather data formats
    - Re-usable calculation routines
    - Weather plots
- Intended for use by anyone in the field: Education, Research, etc.
- Goal is to become a community resource to find useful bits and pieces

Design philosophy:
- Fit well with the rest of the scientific Python ecosystem (NumPy, Matplotlib, etc.)
- Simple to use with your own data
- Unit correctness built-in (using `pint`)
- Well-tested and documented
- Development driven by use cases

#What can MetPy do?

## Reading Data
Pure python implementations for NEXRAD files:
- Level 2
- Level 3

## Example of using NIDS decoding
End of explanation
"""


fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(1, 1, 1)
norm, cmap = ctables.registry.get_with_boundaries('NWSReflectivity', np.arange(0, 85, 5))
pc, data = plot_prod(n0q, cmap, norm, ax)
plot_points(tvs)
plot_points(nmd)
plot_points(nhi)
plot_tracks(nst)
ax.set_ylim(-20, 40)
ax.set_xlim(-50, 20)

"""
Explanation: What about the "data" content of the products?
End of explanation
"""


import metpy.calc as mpcalc
from metpy.units import units
temp = 86 * units.degF
press = 860. * units.mbar
humidity = 40 / 100.

dewpt = mpcalc.dewpoint_rh(temp, humidity).to('degF')
dewpt

"""
Explanation: Calculations

Basic calculations (e.g. heat index, wind components)
Kinematic calculations (e.g. vertical vorticity)
Thermodynamic calculations (e.g. dewpoint, LCL)

Turbulence (e.g. turbulence kinetic energy)


Calculations simple and independent to promote re-use

Only requirement is to use unit support

Example
It's 86F outside here (pressure 860 mb), with a relative humidity of 40%. What's the dewpoint?
End of explanation
"""


mpcalc.lcl(press, temp, dewpt)

"""
Explanation: What does the LCL look like for that?
End of explanation
"""


import numpy as np
pressure_levels = np.array([860., 850., 700., 500., 300.]) * units.mbar
mpcalc.parcel_profile(pressure_levels, temp, dewpt)

"""
Explanation: Given those conditions, what does the profile of a parcel look like?
End of explanation
"""


import matplotlib.pyplot as plt
import numpy as np

from metpy.cbook import get_test_data
from metpy.calc import get_wind_components
from metpy.plots import SkewT

# Parse the data
p, T, Td, direc, spd = np.loadtxt(get_test_data('sounding_data.txt'),
        usecols=(0, 2, 3, 6, 7), skiprows=4, unpack=True)
u, v = get_wind_components(spd, np.deg2rad(direc))

# Create a skewT using matplotlib's default figure size
fig = plt.figure(figsize=(8, 8))
skew = SkewT(fig)

# Plot the data using normal plotting functions, in this case using
# log scaling in Y, as dictated by the typical meteorological plot
skew.plot(p, T, 'r')
skew.plot(p, Td, 'g')
skew.plot_barbs(p, u, v)

# Add the relevant special lines
skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()
skew.ax.set_ylim(1000, 100)
fig

"""
Explanation: Plots

Colortables
Skew-T

Example
Plot some sounding data
End of explanation
"""


from siphon.catalog import TDSCatalog
best_gfs = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/Global_0p5deg/catalog.xml?dataset=grib/NCEP/GFS/Global_0p5deg/Best')
best_ds = list(best_gfs.datasets.values())[0]

best_ds.access_urls

"""
Explanation: Siphon: Sucking data down the pipe
THREDDS?

Unidata technology for serving many different data formats
Creates "feature collections" on top of data files to provide access across multiple files
Provides a variety of access services on top of these datasets
Demonstration server: http://thredds.ucar.edu/thredds/catalog.html

Siphon is Unidata's Python library for talking to THREDDS:
- Making sense of catalogs
- Talking to the NetCDF Subset Service (NCSS)
- Talking to the radar query service (radar server)
Example
Let's use Siphon to get the forecast from the GFS. Start by finding it on THREDDS.
End of explanation
"""


# Set up a class to access
from siphon.ncss import NCSS
ncss = NCSS(best_ds.access_urls['NetcdfSubset'])

# Get today's date
from datetime import datetime, timedelta
now = datetime.utcnow()

# Get a query object and set to get temperature for Boulder for the next 7 days
query = ncss.query()
query.lonlat_point(-105, 40).vertical_level(100000).time_range(now, now + timedelta(days=7))
query.variables('Temperature_isobaric').accept('netcdf4')

# Get the Data
data = ncss.get_data(query)
list(data.variables)

# Pull the variables we want from the NetCDF file
temp = data.variables['Temperature_isobaric']
time = data.variables['time']

# Convert time values from numbers to datetime
from netCDF4 import num2date
time_vals = num2date(time[:].squeeze(), time.units)

# Plot
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 1, figsize=(9, 8))
ax.plot(time_vals, temp[:].squeeze(), 'r', linewidth=2)
ax.set_ylabel(temp.standard_name + ' (%s)' % temp.units)
ax.set_xlabel('Forecast Time (UTC)')
ax.grid(True)

"""
Explanation: Ok, let's get data using NCSS
End of explanation
"""


# Re-use the NCSS access from earlier, but make a new query
query = ncss.query()
query.lonlat_point(-105, 40).time(now + timedelta(hours=12)).accept('csv')
query.variables('Temperature_isobaric', 'Relative_humidity_isobaric',
                'u-component_of_wind_isobaric', 'v-component_of_wind_isobaric')
data = ncss.get_data(query)

# Pull out data with some units
p = (data['vertCoord'] * units('Pa')).to('mbar')
T = data['Temperature_isobaric'] * units('kelvin')
Td = mpcalc.dewpoint_rh(T, data['Relative_humidity_isobaric'] / 100.)
u = data['ucomponent_of_wind_isobaric'] * units('m/s')
v = data['vcomponent_of_wind_isobaric'] * units('m/s')

# Create a skewT using matplotlib's default figure size
fig = plt.figure(figsize=(8, 8))
skew = SkewT(fig)

# Plot the data using normal plotting functions, in this case using
# log scaling in Y, as dictated by the typical meteorological plot
skew.plot(p, T.to('degC'), 'r')
skew.plot(p, Td.to('degC'), 'g')

skew.plot_barbs(p[p>=100 * units.mbar], u.to('knots')[p>=100 * units.mbar],
                v.to('knots')[p>=100 * units.mbar])
skew.ax.set_ylim(1000, 100)
skew.ax.set_title(data['date'][0]);

"""
Explanation: "With your powers combined..."

Siphon is a complementary technology to Metpy
Provides an easy way to get data for use in MetPy

Example
Let's plot a forecast sounding from GFS for Boulder 12 hours from now
End of explanation
"""


# This code isn't actually in a MetPy release yet.....
import importlib
import metpy.plots.station_plot
importlib.reload(metpy.plots.station_plot)
station_plot = metpy.plots.station_plot.station_plot

# Set up query for point data
ncss = NCSS('http://thredds.ucar.edu/thredds/ncss/nws/metar/ncdecoded/Metar_Station_Data_fc.cdmr/dataset.xml')
query = ncss.query()
query.lonlat_box(44, 37, -98, -108).time(now - timedelta(days=2)).accept('csv')
query.variables('air_temperature', 'dew_point_temperature',
                'wind_from_direction', 'wind_speed')
data = ncss.get_data(query)

# Some unit conversions
speed = data['wind_speed'] * units('m/s')
data['u'],data['v'] = mpcalc.get_wind_components(speed.to('knots'),
                                                 data['wind_from_direction'] * units('degree'))

# Plot using basemap for now
from mpl_toolkits.basemap import Basemap
fig = plt.figure(figsize=(9, 9))
ax = fig.add_subplot(1, 1, 1)
m = Basemap(lon_0=-105, lat_0=40, lat_ts=40, resolution='i',
            projection='stere', urcrnrlat=44, urcrnrlon=-98, llcrnrlat=37,
            llcrnrlon=-108, ax=ax)
m.bluemarble()

# Just an early prototype...
station_plot(data, proj=m, ax=ax, layout={'NW': 'air_temperature', 'SW': 'dew_point_temperature'},
             styles={'air_temperature': dict(color='r'), 'dew_point_temperature': dict(color='lightgreen')},
             zorder=1);
fig

"""
Explanation: Ideas for the future
A sampling from our issue tracker:
- Station plots (with weather symbols)
- Hodograph
- Pulling in sounding data from Wyoming archive
End of explanation
"""


