diff --git a/examples/HurricaneTracker.py b/examples/HurricaneTracker.py index f3ceeb7a..8e2c0e2e 100644 --- a/examples/HurricaneTracker.py +++ b/examples/HurricaneTracker.py @@ -1,119 +1,75 @@ """ -Hurricane Tracker with NHC Data -=============================== - +GUI for National Hurricane Center Data +====================================== By: Aodhan Sweeney -This program is a recreation of the 2014 hur_tracker.py -originally written by Unidata Intern Florita Rodriguez. The -2019 version comes with updated interface and functionality, -as well as changing certain dependencies. - +This GUI takes data pulled from the National +Hurricane Center and plots specific storms and +their tracks. The GUI runs using ipywidgets, and +is intended to be used as an Jupyter notebook. """ -import gzip -from io import BytesIO -from io import StringIO - import cartopy.crs as ccrs +from IPython.display import display import ipywidgets as widgets import matplotlib.pyplot as plt import numpy as np -from pandas import DataFrame -import requests - - -def readUrlFile(url): - """readUrlFile is a function created to read a .dat file from a given url - and compile it into a list of strings with given headers.""" - - headers = {'User-agent': 'Unidata Python Client Test'} - response = requests.get(url, headers=headers) - # Store data response in a string buffer - string_buffer = StringIO(response.text) - # Read from the string buffer as if it were a physical file - data = string_buffer.getvalue() - return data.splitlines() - - -def readGZFile(url): - """readGZFile is a function which opens and reads zipped files. In this case it takes in a - .gzfile containing information on each storm and returns a byte buffer split based on - lines.""" - - headers = {'User-agent': 'Unidata Python Client Test'} - response = requests.get(url, headers=headers) - # Store data response in a bytes buffer - bio_buffer = BytesIO(response.content) - # Read from the string buffer as if it were a physical file - gzf = gzip.GzipFile(fileobj=bio_buffer) - data = gzf.read() - return data.splitlines() - - -def split_storm_info(storm_list): - """split_storm_info takes a list of strings and creates a pandas dataframe - for the data set taken off the NHC archive. This function is called in the main to - find all storms.""" +from siphon.simplewebservice import nhc - name, cycloneNum, year, stormType, basin, filename = [], [], [], [], [], [] - for line in storm_list[1:]: - fields = line.split(',') - name.append(fields[0].strip()) - basin.append(fields[1].strip()) - cycloneNum.append(fields[7].strip()) - year.append(fields[8].strip()) - stormType.append(fields[9].strip()) - filename.append(fields[-1].strip().lower()) - storms = DataFrame({'Name': name, 'Basin': basin, 'CycloneNum': np.array(cycloneNum), - 'Year': np.array(year), 'StormType': stormType, - 'Filename': filename}) - return(storms) +class NHC_GUI: + """ + Graphic User Interface designed to allow users to access National Hurricane Center data. - -class Storm_Selection_gui: - """Storm_Selection_gui is a graphic user interface object designed for selecting storms - from the National Hurricane Center's storm list database.""" + This class uses ipython widgets, and the order in which the functions appear in this script + correspond to the order in which the functions and widgets are called/used. + """ def __init__(self): - """__init__ is the initiation function for the Storm_Selection_gui object. This - initiation creates a dataframe for the storm data from the National Hurricane - Center. In addition, this initiation creates widgets to select a specific storm - from the National Hurricane Center database with both a widget for a year - selection and also for a track button which actually retrieves the models - and tracks for a given storm. """ - - # Setting up storm object table - fileLines = readUrlFile('http://ftp.nhc.noaa.gov/atcf/index/storm_list.txt') - self.storm_table = split_storm_info(fileLines) - - # Creation of widgets - # Year Selector Slider (year_slider) + """ + Create object that references NHC.py and thus the National Hurricane Center. + + This initiation creates the National Hurricane Center object and also creates + a widget that allows the user to 1.) select the year in which to find storms, + and 2.) indicate when they have chosen the year and when to continue with parsing + the storm_table. + """ + self.NHCD = nhc.NHCD() + self.storm_table = self.NHCD.storm_table + # Year Slider Widget to select year for which to retrieve storm data. self.year_slider = widgets.IntSlider(min=1851, max=2019, value=2019, description='Storm Year: ') widgets.interact(self.get_storms_slider, year_slider=self.year_slider) - # Button to retrieve storm tracks (track_button) + # Storm Track toggle button to initiate storm track retrieval. self.track_button = widgets.ToggleButton(value=False, description='Get Storm Tracks', disabled=False, button_style='info', tooltip='Description') - widgets.interact(self.get_tracks, track_button=self.track_button) + widgets.interact(self.get_track, track_button=self.track_button) def get_storms_slider(self, year_slider): - """get_storms_slider is a function that is written to take a user defined year from - 1851 to 2019 and then trim the split_storms_info dataframe to just have - storms from that year. This takes in the interactive widget called year_slider.""" - - self.one_year_table = self.storm_table[self.storm_table.Year == str(year_slider)] + """ + Take a year chosen by the user, and create a list of all storms during that year. + + Parameters + ---------- + year_slider: ipywidget + tells value of year chosen by the user + """ + self.year = str(year_slider) + self.one_year_table = self.storm_table[self.storm_table.Year == year_slider] self.storm_names = widgets.Dropdown(options=self.one_year_table['Name'], description='Storm Names: ') widgets.interact(self.get_name_dropdown, storm_names=self.storm_names) def get_name_dropdown(self, storm_names): - """get_name_dropdown is a function that allows for selection of a specific storm - based on its name. This function returns the filename of a given storm. This - function interacts with the widget called forecast_select.""" - + """ + Take names created by previous function and allow selection of which models to plot. + + Parameters + ---------- + storm_names: list + all storms in a given year + """ name = self.storm_names.value one_storm_row = self.one_year_table[self.one_year_table.Name == name] self.filename = one_storm_row.Filename @@ -121,182 +77,103 @@ def get_name_dropdown(self, storm_names): if self.filename.empty is False: self.filename = file_name[0] elif self.filename.empty is True: - print('No file name data for this year.') - - def get_models(self, models): - """get_models is a function that is linked to the selection widget which allows for - multiple models to be selected to plot for a specific hurricane track.""" - - self.model_selection_latlon(models) - - def time_slider(self, time_slide): - """time_slider is a function which changes the time for which to plot the models for - a given hurricane track.""" - - time = self.date_times[time_slide] - return(time) - - def get_tracks(self, track_button): - """get_tracks is a function that will create the url and pull the data for either - the forecast track or best track for a given storm. The Url is made by using both - the year and the filename. This function will then read the data and create a data - frame for both the forecast and best tracks and compile these data frames into a - dictionary. This function returns this dictionary of forecast and best tracks. """ - - year = str(self.year_slider.value) - filename = self.filename - data_dictionary = {} - # Current year data is stored in a different location - if year == '2019': - urlf = 'http://ftp.nhc.noaa.gov/atcf/aid_public/a{}.dat.gz'.format(filename) - urlb = 'http://ftp.nhc.noaa.gov/atcf/btk/b{}.dat'.format(filename) - else: - urlf = 'http://ftp.nhc.noaa.gov/atcf/archive/{}/a{}.dat.gz'.format(year, filename) - urlb = 'http://ftp.nhc.noaa.gov/atcf/archive/{}/b{}.dat.gz'.format(year, filename) - - url_links = [urlf, urlb] - url_count = 0 - if track_button is True: - for url in url_links: - # Checking if url is valid, if status_code is 200 then website is active - if requests.get(url).status_code == 200: - if url.endswith('.dat'): - lines = readUrlFile(url) - else: - lines = readGZFile(url) - - # Splitting the method for which we will create the dataframe - lat, lon, basin, cycloneNum = [], [], [], [] - warnDT, model, forecast_hour = [], [], [] - for line in lines: - line = str(line) - line = line[2:] - fields = line.split(',') - # Joins together lattitude and longitude strings without - # directional letters. - # Includes type conversion in order to divide by 10 to - # get the correct coordinate. - latSingle = int(fields[6][:-1])/10.0 - lonSingle = -(int(fields[7][:-1])/10.0) - lat.append(latSingle) - lon.append(lonSingle) - basin.append(fields[0]) - forecast_hour.append(fields[5]) - cycloneNum.append(fields[1].strip()) - warnDT.append(fields[2].strip()) - model.append(fields[4].strip()) - - # Combining data from file into a Pandas Dataframe. - storm_data_frame = DataFrame({'Basin': basin, - 'CycloneNum': np.array(cycloneNum), - 'WarnDT': np.array(warnDT), - 'Model': model, 'Lat': np.array(lat), - 'Lon': np.array(lon), - 'forecast_hour': - np.array(forecast_hour)}) - # Adding this newly created DataFrame to a dictionary - if url_count == 0: - data_dictionary['forecast'] = storm_data_frame - else: - data_dictionary['best_track'] = storm_data_frame - - else: - print('url {} was not valid, select different storm.'.format(url)) - track_button = False - - url_count += 1 - - self.storm_dictionary = data_dictionary - forecast = data_dictionary.get('forecast') - - unique_models, unique_index = list(np.unique(forecast['Model'].values, - return_index=True)) - - # Selection tool to pick from available models to plot (model_select) + raise ValueError('No data for specific file name of the chosen.') + + def get_track(self, track_button): + """ + Query whether track button has been toggled, and create select model widget. + + Parameters + ---------- + track_button: ipywidget + button that when toggled indicates that user is ready to select the model + tracks for the chosen storm. + """ + if self.track_button.value is True: + unique_models = self.NHCD.get_tracks(self.year, self.filename) self.model_select = widgets.SelectMultiple(options=unique_models, value=[unique_models[0]], description='Models: ', disabled=False) widgets.interact(self.get_models, models=self.model_select) - def model_selection_latlon(self, models): - """model_selection_latlon is a function that allows the user to select a model for - a given storm and whether the tracks are forecast or best tracks. The parameters - for this are a string stating whether the user wants forecast or best tracks and - also all model outputs for all forecasts and best tracks compiled into a python - dictionary. The latlon part of this function comes from taking the users selected - model and getting the latitudes and longitudes of all positions of the storm for - this forecast. This function then returns these lats and lons as a pandas.Series""" - - if self.model_select.disabled is False: - # We will always plot best track, and thus must save the coordinates for plotting - best_track = self.storm_dictionary.get('best_track') - self.date_times = best_track['WarnDT'] - lats = best_track['Lat'] - lons = best_track['Lon'] - self.best_track_coordinates = [lats, lons] - - model_tracks = self.storm_dictionary.get('forecast') - - self.model_table = [] - for model in models: - one_model_table = model_tracks[model_tracks['Model'] == model] - self.model_table.append(one_model_table) - - # Slider to choose time frame for which to plot models (plot_slider) - self.plot_slider = widgets.IntSlider(min=0, max=(len(self.date_times)-1), - value=0, description='Tracks Time', - disabled=False) - widgets.interact(self.plotting, plot_slider=self.plot_slider) + def get_models(self, models): + """ + Select models from a list of all models for a given stormself. + + Parameters + ---------- + models: list + list of ran for a given storm + """ + self.NHCD.model_selection_latlon(models) + self.date_times = self.NHCD.date_times.tolist() + self.plot_slider = widgets.IntSlider(min=0, max=(len(self.date_times)-1), + value=0, description='Tracks Time', + disabled=False) + self.play = widgets.Play(interval=800, min=0, max=(len(self.date_times)-1), + value=0, description='Tracks Time') + widgets.jslink((self.plot_slider, 'value'), (self.play, 'value')) + self.box = widgets.HBox([self.plot_slider, self.play]) + display(self.plot_slider) + widgets.interact(self.plotting, plot_slider=self.play) def plotting(self, plot_slider): - """plotting is a function that that plots the models tracks and best tracks for all - selected models. These tracks are then plotted on the Blue Marble projection. """ - + """ + Plot selected model tracks and best track for given storm. + + Parameters + ---------- + plot_slider: ipywidget + Widget where assigned value is plot/play widget value + """ if self.plot_slider.disabled is False: + # Identifying the time associated with the models for time text box year = self.date_times[plot_slider][0: 4] month = self.date_times[plot_slider][4: 6] day = self.date_times[plot_slider][6: 8] hour = self.date_times[plot_slider][8: 10] - time_string = 'Date: {0}/{1}/{2} \n Hour: {3}'.format(month, day, year, hour) + time_string = 'Date: {0}/{1}/{2} \nHour: {3}'.format(month, day, year, hour) # Finding data for best track, and extremes for which to base axis extent on - self.best_lats = np.array(self.best_track_coordinates[0]) - self.best_lons = np.array(self.best_track_coordinates[1]) + self.best_lats = np.array(self.NHCD.best_track_coordinates[0]) + self.best_lons = np.array(self.NHCD.best_track_coordinates[1]) min_best_lat = min(self.best_lats) max_best_lat = max(self.best_lats) min_best_lon = min(self.best_lons) max_best_lon = max(self.best_lons) - # Plotting the track on a cartopy stock image + # Plotting the tracks on top of a cartopy stock image projection self.fig = plt.figure(figsize=(14, 11)) self.ax = self.fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) self.ax.stock_img() self.data_projection = ccrs.PlateCarree() self.ax.plot(self.best_lons, self.best_lats, marker='o', color='white', label='Best Track', transform=self.data_projection) - self.ax.set_extent([(min_best_lon - 30), (max_best_lon + 30), - (min_best_lat - 30), (max_best_lat + 30)]) + + self.ax.set_extent([(min_best_lon - 40), (max_best_lon + 40), + (min_best_lat - 40), (max_best_lat + 40)]) + jet = plt.get_cmap('jet') colors = iter(jet(np.linspace(0.2, 1, (len(self.model_select.value)+1)))) left = .1 bottom = .1 self.ax.text(left, bottom, time_string, transform=self.ax.transAxes, fontsize=14, color='black') - for model_type in self.model_table: - one_model_time = model_type[model_type['WarnDT'] == - self.date_times[plot_slider]] - lats = one_model_time['Lat'].tolist() - lons = one_model_time['Lon'].tolist() + + for model_type in self.NHCD.model_table: + one_time = model_type[model_type['WarnDT'] == self.date_times[plot_slider]] + lats = one_time['Lat'].tolist() + lons = one_time['Lon'].tolist() if len(lats) != 0: model_list = model_type['Model'].tolist() self.ax.plot(lons, lats, marker='o', color=next(colors), label=model_list[0]) plt.title('Storm Name: {0} Year: {1}'.format(self.storm_names.value, - str(self.year_slider.value))) + self.year)) plt.legend() ###################################################################### -storm_selection = Storm_Selection_gui() +NHC_GUI() diff --git a/examples/SPCTracker.py b/examples/SPCTracker.py new file mode 100644 index 00000000..6facc9de --- /dev/null +++ b/examples/SPCTracker.py @@ -0,0 +1,160 @@ +""" +GUI for Storm Prediction Center +=============================== + +This GUI takes in SPC data and plots +all events from a day. It is meant +to show the ease of data access and +creation of a useful plot using the +Siphon Simple Web Service for the SPC. +""" + + +import cartopy.crs as ccrs +import cartopy.feature as cfeature +from IPython.display import display +import ipywidgets as widgets +import matplotlib.patches as mpatches +import matplotlib.pyplot as plt +import pandas as pd +from siphon.simplewebservice import spc + + +class SPC_GUI: + """ + Graphic User Interface designed to allow users to access Storm Prediction Center data. + + This class uses ipython widgets. NOTE: date chosen must be from 2011 to present because + methods of data parsing are different and will not work for dates earlier. + """ + def __init__(self): + """ + Create object that references SPC.py and thus the Storm Prediction Center. + + This initiation creates the SPC object and also creates a widget that allows + the user to select a date for which to animate SPC reports from the years of + 2011 and onward. + """ + self.datepicker = widgets.DatePicker(description='Pick a Date:', disabled=False) + widgets.interact(self.format_date, datepicker=self.datepicker) + + def format_date(self, datepicker): + """ + Allow user to chose a date for which to plot the SPC events. + + Parameters + ========== + datepicker: ipywidget + allows chosing of date + """ + if datepicker is not None: + year_string = datepicker.strftime('%Y') + month_string = datepicker.strftime('%m') + day_string = datepicker.strftime('%d') + self.date_string = year_string + month_string + day_string + self.stormpicker = widgets.SelectMultiple(options=['tornado', 'hail', 'wind'], + description='Event Type: ') + widgets.interact(self.fetch_spc, stormpicker=self.stormpicker) + + def fetch_spc(self, stormpicker): + """ + Use a date chosen by the user and create widgets that allow you to iterate through + the pandas dataframe holding the data. + + Parameters + ========== + stormpicker: ipywidget + allows choice of the storm type to be plotted + """ + self.spc_data_table = [] + for storm_type in stormpicker: + one_storm_type = spc.SPC(storm_type, self.date_string) + # Storm must be after year 2011 + if int(self.date_string[:4]) < 2011: + raise ValueError('SPC GUI does not support events before 2012.') + self.spc_data_table.append(one_storm_type.day_table) + + if self.spc_data_table != []: + self.plot_button = widgets.ToggleButton(value=False, description='Plot SPC Events', + disabled=False, button_style='danger', + tooltip='Description') + + self.all_events = pd.concat(self.spc_data_table, sort=False) + self.all_events = self.all_events.sort_index() + + if self.all_events.empty: + raise ValueError('No storm reports of any type for this date.') + + self.plot_slider = widgets.IntSlider(min=0, max=(len(self.all_events)-1), + value=0, description='Event #', + disabled=False) + self.play = widgets.Play(interval=900, min=0, max=(len(self.all_events)-1), + value=0, description='Event #') + widgets.jslink((self.plot_slider, 'value'), (self.play, 'value')) + self.plot_button = widgets.ToggleButton(value=False, + description='Plot SPC Events', + disabled=False, button_style='danger', + tooltip='description') + self.box = widgets.HBox([self.plot_slider, self.play]) + display(self.plot_slider) + widgets.interact(self.plot_events, plot_slider=self.play, + plot_button=self.plot_button) + + def plot_events(self, plot_slider, plot_button): + """ + Plot the storm events chosen for the date picked by the user and the type of storms + selected. + + Parameters + ========== + plot_sider: ipywidget + pases a iteration of the dataframe holding all storm data to be plotted + + plot_button: ipywidget + toggle button that when activated allows for plotting + """ + spc_event = self.all_events.iloc[plot_slider] + if self.plot_button.value is True: + trimmed_event = spc_event.dropna() + # Plotting the tracks on top of a cartopy stock image projection + self.fig = plt.figure(figsize=(14, 11)) + self.ax = self.fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) + self.ax.stock_img() + state_lines = 'admin_1_states_provinces_lines' + states_provinces = cfeature.NaturalEarthFeature(category='cultural', + name=state_lines, + scale='50m', + facecolor='none') + self.ax.add_feature(cfeature.BORDERS) + self.ax.add_feature(cfeature.COASTLINE) + self.ax.add_feature(cfeature.LAKES) + + self.ax.add_feature(states_provinces, edgecolor='gray') + self.ax.set_title('SPC Events for {}/{}/{}'.format(self.date_string[4:6], + self.date_string[6:8], + self.date_string[0:4])) + self.data_projection = ccrs.PlateCarree() + if 'Size (in)' in trimmed_event: + self.ax.plot(spc_event['Lon'], spc_event['Lat'], marker='o', color='blue', + label='Hail', transform=self.data_projection, markersize=7) + if 'Speed (kt)' in trimmed_event: + self.ax.plot(spc_event['Lon'], spc_event['Lat'], marker='o', color='green', + label='Wind', transform=self.data_projection, markersize=7) + if 'F-Scale' in trimmed_event: + self.ax.plot(spc_event['Lon'], spc_event['Lat'], marker='o', color='red', + label='Tornado', transform=self.data_projection, markersize=7) + self.ax.set_extent([-60, -130, 23, 50]) + hail_event = mpatches.Patch(color='blue', label='Hail') + wind_event = mpatches.Patch(color='green', label='Wind') + torn_event = mpatches.Patch(color='red', label='Tornado') + self.ax.legend(handles=[hail_event, wind_event, torn_event]) + comment = spc_event['Comment'] + if len(comment) > 70: + comment = comment[:70] + '\n' + comment[70:] + if len(comment) > 140: + comment = comment[:140] + '\n' + comment[140:] + self.fig.text(0.15, 0.2, comment, fontsize=15) + + +###################################################################### +SPC_GUI() diff --git a/examples/Upperair_Obs.py b/examples/Upperair_Obs.py new file mode 100644 index 00000000..02a5f818 --- /dev/null +++ b/examples/Upperair_Obs.py @@ -0,0 +1,393 @@ +""" +DIFAX Replication +================= + +This example replicates the traditional DIFAX images for upper-level +observations. + +By: Kevin Goebbert + +Observation data comes from Iowa State Archive, accessed through the +Siphon package. Contour data comes from the GFS 0.5 degree analysis. +Classic upper-level data of Geopotential Height and Temperature are +plotted. + +""" + +from datetime import datetime, timedelta +import urllib.request + +import cartopy.crs as ccrs +import cartopy.feature as cfeature +import matplotlib.pyplot as plt +import metpy.calc as mpcalc +from metpy.plots import StationPlot +from metpy.units import units +import numpy as np +from siphon.simplewebservice.iastate import IAStateUpperAir +import xarray as xr + + +###################################################################### +# Plotting High/Low Symbols +# ------------------------- +# +# A helper function to plot a text symbol (e.g., H, L) for relative +# maximum/minimum for a given field (e.g., geopotential height). +# + + +def plot_maxmin_points(lon, lat, data, extrema, nsize, symbol, color='k', + plotValue=True, transform=None): + """ + This function will find and plot relative maximum and minimum for a 2D grid. The function + can be used to plot an H for maximum values (e.g., High pressure) and an L for minimum + values (e.g., low pressue). It is best to used filetered data to obtain a synoptic scale + max/min value. The symbol text can be set to a string value and optionally the color of the + symbol and any plotted value can be set with the parameter color. + + Parameters + ---------- + lon : 2D array + Plotting longitude values + lat : 2D array + Plotting latitude values + data : 2D array + Data that you wish to plot the max/min symbol placement + extrema : str + Either a value of max for Maximum Values or min for Minimum Values + nsize : int + Size of the grid box to filter the max and min values to plot a reasonable number + symbol : str + Text to be placed at location of max/min value + color : str + Name of matplotlib colorname to plot the symbol (and numerical value, if plotted) + plot_value : Boolean (True/False) + Whether to plot the numeric value of max/min point + + Return + ------ + The max/min symbol will be plotted on the current axes within the bounding frame + (e.g., clip_on=True) + """ + from scipy.ndimage.filters import maximum_filter, minimum_filter + + if (extrema == 'max'): + data_ext = maximum_filter(data, nsize, mode='nearest') + elif (extrema == 'min'): + data_ext = minimum_filter(data, nsize, mode='nearest') + else: + raise ValueError('Value for hilo must be either max or min') + + if lon.ndim == 1: + lon, lat = np.meshgrid(lon, lat) + + mxx, mxy = np.where(data_ext == data) + + for i in range(len(mxy)): + ax.text(lon[mxx[i], mxy[i]], lat[mxx[i], mxy[i]], symbol, color=color, size=36, + clip_on=True, horizontalalignment='center', verticalalignment='center', + transform=transform) + ax.text(lon[mxx[i], mxy[i]], lat[mxx[i], mxy[i]], + '\n' + str(np.int(data[mxx[i], mxy[i]])), + color=color, size=12, clip_on=True, fontweight='bold', + horizontalalignment='center', verticalalignment='top', transform=transform) + ax.plot(lon[mxx[i], mxy[i]], lat[mxx[i], mxy[i]], marker='o', markeredgecolor='black', + markerfacecolor='white', transform=transform) + ax.plot(lon[mxx[i], mxy[i]], lat[mxx[i], mxy[i]], + marker='x', color='black', transform=transform) + + +###################################################################### +# Station Information +# ------------------- +# +# A helper function for obtaining radiosonde station information (e.g., +# latitude/longitude) requried to plot data obtained from each station. +# Original code by github user sgdecker. +# + + +def station_info(stid): + r"""Provide information about weather stations. + + Parameters + ---------- + stid: str or iterable object containing strs + The ICAO or IATA code(s) for which station information is requested. + with_units: bool + Whether to include units for values that have them. Default True. + + Returns + ------- + info: dict + Information about the station(s) within a dictionary with these keys: + 'state': Two-character ID of the state/province where the station is located, + if applicable + 'name': The name of the station + 'lat': The latitude of the station [deg] + 'lon': The longitude of the station [deg] + 'elevation': The elevation of the station [m] + 'country': Two-character ID of the country where the station is located + + Modified code from Steven Decker, Rutgers University + + """ + # Provide a helper function for later usage + def str2latlon(s): + deg = float(s[:3]) + mn = float(s[-3:-1]) + if s[-1] == 'S' or s[-1] == 'W': + deg = -deg + mn = -mn + return deg + mn / 60. + + # Various constants describing the underlying data + url = 'https://www.aviationweather.gov/docs/metar/stations.txt' + # file = 'stations.txt' + state_bnds = slice(0, 2) + name_bnds = slice(3, 19) + icao_bnds = slice(20, 24) + iata_bnds = slice(26, 29) + lat_bnds = slice(39, 45) + lon_bnds = slice(47, 54) + z_bnds = slice(55, 59) + cntry_bnds = slice(81, 83) + + # Generalize to any number of IDs + if isinstance(stid, str): + stid = [stid] + + # Get the station dataset + infile = urllib.request.urlopen(url) + data = infile.readlines() +# infile = open(file, 'rb') +# data = infile.readlines() + + state = [] + name = [] + lat = [] + lon = [] + z = [] + cntry = [] + + for s in stid: + s = s.upper() + for line_bytes in data: + line = line_bytes.decode('UTF-8') + icao = line[icao_bnds] + iata = line[iata_bnds] + if len(s) == 3 and s in iata or len(s) == 4 and s in icao: + state.append(line[state_bnds].strip()) + name.append(line[name_bnds].strip()) + lat.append(str2latlon(line[lat_bnds])) + lon.append(str2latlon(line[lon_bnds])) + z.append(float(line[z_bnds])) + cntry.append(line[cntry_bnds]) + + break + else: + state.append('NA') + name.append('NA') + lat.append(np.nan) + lon.append(np.nan) + z.append(np.nan) + cntry.append('NA') + + infile.close() + + return {'state': np.array(state), 'name': np.array(name), 'lat': np.array(lat), + 'lon': np.array(lon), 'elevation': np.array(z), 'country': np.array(cntry), + 'units': {'lat': 'deg', 'lon': 'deg', 'z': 'm'}} + + +###################################################################### +# Observation Data +# ---------------- +# +# Set a date and time for upper-air observations (should only be 00 or 12 +# UTC for the hour). +# +# Request all data from Iowa State using the Siphon package. The result is +# a pandas DataFrame containing all of the sounding data from all +# available stations. +# + +# Set date for desired UPA data +today = datetime.utcnow() + +# Go back one day to ensure data availability +date = datetime(today.year, today.month, today.day, 0) - timedelta(days=1) + +# Request data using Siphon request for data from Iowa State Archive +data = IAStateUpperAir.request_all_data(date) + + +###################################################################### +# Subset Observational Data +# ------------------------- +# +# From the request above will give all levels from all radisonde sites +# available through the service. For plotting a pressure surface map there +# is only need to have the data from that level. Below the data is subset +# and a few parameters set based on the level chosen. Additionally, the +# station information is obtained and latitude and longitude data is added +# to the DataFrame. +# + +level = 500 + +if (level == 925) | (level == 850) | (level == 700): + cint = 30 + def hght_format(v): return format(v, '.0f')[1:] +elif level == 500: + cint = 60 + def hght_format(v): return format(v, '.0f')[:3] +elif level == 300: + cint = 120 + def hght_format(v): return format(v, '.0f')[:3] +elif level < 300: + cint = 120 + def hght_format(v): return format(v, '.0f')[1:4] + +# Create subset of all data for a given level +data_subset = data.pressure == level +df = data[data_subset] + +# Get station lat/lon from look-up file; add to Dataframe +stn_info = station_info(list(df.station.values)) +df.insert(10, 'latitude', stn_info['lat']) +df.insert(11, 'longitude', stn_info['lon']) + + +###################################################################### +# Gridded Data +# ------------ +# +# Obtain GFS gridded output for contour plotting. Specifically, +# geopotential height and temperature data for the given level and subset +# for over North America. Data are smoothed for aesthetic reasons. +# + +# Get GFS data and subset to North America for Geopotential Height and Temperature +ds = xr.open_dataset('https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/Global_0p5deg_ana/' + 'GFS_Global_0p5deg_ana_{0:%Y%m%d}_{0:%H}00.grib2'.format( + date)).metpy.parse_cf() + +# Geopotential height and smooth +hght = ds.Geopotential_height_isobaric.metpy.sel( + vertical=level*units.hPa, time=date, lat=slice(70, 15), lon=slice(360-145, 360-50)) +smooth_hght = mpcalc.smooth_n_point(hght, 9, 10) + +# Temperature, smooth, and convert to Celsius +tmpk = ds.Temperature_isobaric.metpy.sel( + vertical=level*units.hPa, time=date, lat=slice(70, 15), lon=slice(360-145, 360-50)) +smooth_tmpc = (mpcalc.smooth_n_point(tmpk, 9, 10)).to('degC') + + +###################################################################### +# Create DIFAX Replication +# ------------------------ +# +# Plot the observational data and contours on a Lambert Conformal map and +# add features that resemble the historic DIFAX maps. +# + +# Set up map coordinate reference system +mapcrs = ccrs.LambertConformal( + central_latitude=45, central_longitude=-100, standard_parallels=(30, 60)) + +# Set up station locations for plotting observations +point_locs = mapcrs.transform_points( + ccrs.PlateCarree(), df['longitude'].values, df['latitude'].values) + +# Start figure and set graphics extent +fig = plt.figure(1, figsize=(17, 15)) +ax = plt.subplot(111, projection=mapcrs) +ax.set_extent([-125, -70, 20, 55]) + +# Add map features for geographic reference +ax.add_feature(cfeature.COASTLINE.with_scale('50m'), edgecolor='grey') +ax.add_feature(cfeature.LAND.with_scale('50m'), facecolor='white') +ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='grey') + +# Plot plus signs every degree lat/lon +plus_lat = [] +plus_lon = [] +other_lat = [] +other_lon = [] + +for x in hght.lon.values[::2]: + for y in hght.lat.values[::2]: + if (x % 5 == 0) | (y % 5 == 0): + plus_lon.append(x) + plus_lat.append(y) + else: + other_lon.append(x) + other_lat.append(y) +ax.scatter(other_lon, other_lat, s=5, marker='o', + transform=ccrs.PlateCarree(), color='lightgrey', zorder=-1) +ax.scatter(plus_lon, plus_lat, s=30, marker='+', + transform=ccrs.PlateCarree(), color='lightgrey', zorder=-1) + +# Add gridlines for every 5 degree lat/lon +ax.gridlines(linestyle='solid', ylocs=range(15, 71, 5), xlocs=range(-150, -49, 5)) + +# Start the station plot by specifying the axes to draw on, as well as the +# lon/lat of the stations (with transform). We also the fontsize to 10 pt. +stationplot = StationPlot(ax, df['longitude'].values, df['latitude'].values, clip_on=True, + transform=ccrs.PlateCarree(), fontsize=10) + +# Plot the temperature and dew point to the upper and lower left, respectively, of +# the center point. +stationplot.plot_parameter('NW', df['temperature'], color='black') +stationplot.plot_parameter('SW', df['dewpoint'], color='black') + +# A more complex example uses a custom formatter to control how the geopotential height +# values are plotted. This is set in an earlier if-statement to work appropriate for +# different levels. +stationplot.plot_parameter('NE', df['height'], formatter=hght_format) + +# Add wind barbs +stationplot.plot_barb(df['u_wind'], df['v_wind'], length=7, pivot='tip') + +# Plot Solid Contours of Geopotential Height +cs = ax.contour(hght.lon, hght.lat, smooth_hght, + range(0, 20000, cint), colors='black', transform=ccrs.PlateCarree()) +clabels = plt.clabel(cs, fmt='%d', colors='white', inline_spacing=5, use_clabeltext=True) + +# Contour labels with black boxes and white text +for t in clabels: + t.set_bbox({'facecolor': 'black', 'pad': 4}) + t.set_fontweight('heavy') + +# Plot Dashed Contours of Temperature +cs2 = ax.contour(hght.lon, hght.lat, smooth_tmpc, range(-60, 51, 5), + colors='black', transform=ccrs.PlateCarree()) +clabels = plt.clabel(cs2, fmt='%d', colors='white', inline_spacing=5, use_clabeltext=True) + +# Set longer dashes than default +for c in cs2.collections: + c.set_dashes([(0, (5.0, 3.0))]) + +# Contour labels with black boxes and white text +for t in clabels: + t.set_bbox({'facecolor': 'black', 'pad': 4}) + t.set_fontweight('heavy') + +# Plot filled circles for Radiosonde Obs +ax.scatter(df['longitude'].values, df['latitude'].values, s=12, + marker='o', color='black', transform=ccrs.PlateCarree()) + +# Use definition to plot H/L symbols +plot_maxmin_points(hght.lon, hght.lat, smooth_hght.m, 'max', 50, + symbol='H', color='black', transform=ccrs.PlateCarree()) +plot_maxmin_points(hght.lon, hght.lat, smooth_hght.m, 'min', 25, + symbol='L', color='black', transform=ccrs.PlateCarree()) + +# Add titles +plt.title('Upper-air Observations at {}-hPa Analysis Heights/Temperature'.format(level), + loc='left') +plt.title('Valid: {}'.format(date), loc='right') + +plt.show()