Working with Surface Observations in Siphon and MetPy

In [ ]:

    

from siphon.catalog import TDSCatalog

# copied from the browser url box
metar_cat_url = ('http://thredds.ucar.edu/thredds/catalog/'
                 'irma/metar/catalog.xml?dataset=irma/metar/Metar_Station_Data_-_Irma_fc.cdmr')

# Parse the xml
catalog = TDSCatalog(metar_cat_url)

# what datasets are here?
print(list(catalog.datasets))


    

In [ ]:

    

metar_dataset = catalog.datasets['Feature Collection']


    

In [ ]:

    

# Can safely ignore the warnings
ncss = metar_dataset.subset()


    

In [ ]:

    

ncss.variables


    

In [ ]:

    

from datetime import datetime

query = ncss.query()
query.lonlat_box(north=34, south=24, east=-80, west=-90)
query.time(datetime(2017, 9, 10, 12))
query.variables('temperature', 'dewpoint', 'altimeter_setting',
                'wind_speed', 'wind_direction', 'sky_coverage')
query.accept('csv')


    

In [ ]:

    

# Get the data
data = ncss.get_data(query)
data


    

In [ ]:

    

import numpy as np

import metpy.calc as mpcalc
from metpy.units import units

# Since we used the CSV data, this is just a dictionary of arrays
lats = data['latitude']
lons = data['longitude']
tair = data['temperature']
dewp = data['dewpoint']
alt = data['altimeter_setting']

# Convert wind to components
u, v = mpcalc.wind_components(data['wind_speed'], data['wind_direction'] * units.degree)

# Need to handle missing (NaN) and convert to proper code
cloud_cover = 8 * data['sky_coverage'] / 100.
cloud_cover[np.isnan(cloud_cover)] = 10
cloud_cover = cloud_cover.astype(np.int)

# For some reason these come back as bytes instead of strings
stid = np.array([s.tostring().decode() for s in data['station']])


    

In [ ]:

    

%matplotlib inline

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

from metpy.plots import StationPlot, sky_cover

# Set up a plot with map features
fig = plt.figure(figsize=(12, 12))
proj = ccrs.Stereographic(central_longitude=-95, central_latitude=35)
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.add_feature(cfeature.STATES, edgecolor='black')
ax.coastlines(resolution='50m')
ax.gridlines()

# Create a station plot pointing to an Axes to draw on as well as the location of points
stationplot = StationPlot(ax, lons, lats, transform=ccrs.PlateCarree(),
                          fontsize=12)
stationplot.plot_parameter('NW', tair, color='red')

# Add wind barbs
stationplot.plot_barb(u, v)

# Plot the sky cover symbols in the center. We give it the integer code values that
# should be plotted, as well as a mapping class that can convert the integer values
# to the appropriate font glyph.
stationplot.plot_symbol('C', cloud_cover, sky_cover)


    

In [ ]:

    

# Project points so that we're filtering based on the way the stations are laid out on the map
proj = ccrs.Stereographic(central_longitude=-95, central_latitude=35)
xy = proj.transform_points(ccrs.PlateCarree(), lons, lats)

# Reduce point density so that there's only one point within a 200km circle
mask = mpcalc.reduce_point_density(xy, 200000)


    

In [ ]:

    

# Set up a plot with map features
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.add_feature(cfeature.STATES, edgecolor='black')
ax.coastlines(resolution='50m')
ax.gridlines()

# Create a station plot pointing to an Axes to draw on as well as the location of points
stationplot = StationPlot(ax, lons[mask], lats[mask], transform=ccrs.PlateCarree(),
                          fontsize=12)
stationplot.plot_parameter('NW', tair[mask], color='red')
stationplot.plot_barb(u[mask], v[mask])
stationplot.plot_symbol('C', cloud_cover[mask], sky_cover)


    

In [ ]:

    

# Use reduce_point_density

# Set up a plot with map features
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.add_feature(cfeature.STATES, edgecolor='black')
ax.coastlines(resolution='50m')
ax.gridlines()

# Create a station plot pointing to an Axes to draw on as well as the location of points

# Plot dewpoint

# Plot altimeter setting--formatter can take a function that formats values

# Plot station id


    

In [ ]:

    

from datetime import timedelta

# define the time range we are interested in
end_time = datetime(2017, 9, 12, 0)
start_time = end_time - timedelta(days=2)

# build the query
query = ncss.query()
query.lonlat_point(-80.25, 25.8)
query.time_range(start_time, end_time)
query.variables('altimeter_setting', 'temperature', 'dewpoint',
                'wind_direction', 'wind_speed')
query.accept('csv')


    

In [ ]:

    

data = ncss.get_data(query)


    

In [ ]:

    

print(list(data.keys()))


    

In [ ]:

    

station_id = data['station'][0].tostring()
print(station_id)


    

In [ ]:

    

station_id = station_id.decode('ascii')
print(station_id)


    

In [ ]:

    

data['time']


    

In [ ]:

    

time = [datetime.strptime(s.decode('ascii'), '%Y-%m-%dT%H:%M:%SZ') for s in data['time']]


    

In [ ]:

    

from matplotlib.dates import DateFormatter, AutoDateLocator

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(time, data['wind_speed'], color='tab:blue')

ax.set_title(f'Site: {station_id}      Date: {time[0]:%Y/%m/%d}')
ax.set_xlabel('Hour of day')
ax.set_ylabel('Wind Speed')
ax.grid(True)

# Improve on the default ticking
locator = AutoDateLocator()
hoursFmt = DateFormatter('%H')
ax.xaxis.set_major_locator(locator)
ax.xaxis.set_major_formatter(hoursFmt)


    

In [ ]:

    

# Your code goes here