1# Copyright (c) 2016,2017 MetPy Developers.
2# Distributed under the terms of the BSD 3-Clause License.
3# SPDX-License-Identifier: BSD-3-Clause
4"""
5Station Plot with Layout
6========================
7
8Make a station plot, complete with sky cover and weather symbols, using a
9station plot layout built into MetPy.
10
11The station plot itself is straightforward, but there is a bit of code to perform the
12data-wrangling (hopefully that situation will improve in the future). Certainly, if you have
13existing point data in a format you can work with trivially, the station plot will be simple.
14
15The `StationPlotLayout` class is used to standardize the plotting various parameters
16(i.e. temperature), keeping track of the location, formatting, and even the units for use in
17the station plot. This makes it easy (if using standardized names) to re-use a given layout
18of a station plot.
19"""
20import cartopy.crs as ccrs
21import cartopy.feature as cfeature
22import matplotlib.pyplot as plt
23import pandas as pd
24
25from metpy.calc import wind_components
26from metpy.cbook import get_test_data
27from metpy.plots import (add_metpy_logo, simple_layout, StationPlot, StationPlotLayout,
28                         wx_code_map)
29from metpy.units import units
30
31###########################################
32# The setup
33# ---------
34#
35# First read in the data. We use `numpy.loadtxt` to read in the data and use a structured
36# `numpy.dtype` to allow different types for the various columns. This allows us to handle
37# the columns with string data.
38with get_test_data('station_data.txt') as f:
39    data_arr = pd.read_csv(f, header=0, usecols=(1, 2, 3, 4, 5, 6, 7, 17, 18, 19),
40                           names=['stid', 'lat', 'lon', 'slp', 'air_temperature',
41                                  'cloud_fraction', 'dew_point_temperature', 'weather',
42                                  'wind_dir', 'wind_speed'],
43                           na_values=-99999)
44
45    data_arr.set_index('stid', inplace=True)
46
47###########################################
48# This sample data has *way* too many stations to plot all of them. Instead, we just select
49# a few from around the U.S. and pull those out of the data file.
50
51# Pull out these specific stations
52selected = ['OKC', 'ICT', 'GLD', 'MEM', 'BOS', 'MIA', 'MOB', 'ABQ', 'PHX', 'TTF',
53            'ORD', 'BIL', 'BIS', 'CPR', 'LAX', 'ATL', 'MSP', 'SLC', 'DFW', 'NYC', 'PHL',
54            'PIT', 'IND', 'OLY', 'SYR', 'LEX', 'CHS', 'TLH', 'HOU', 'GJT', 'LBB', 'LSV',
55            'GRB', 'CLT', 'LNK', 'DSM', 'BOI', 'FSD', 'RAP', 'RIC', 'JAN', 'HSV', 'CRW',
56            'SAT', 'BUY', '0CO', 'ZPC', 'VIH']
57
58# Loop over all the whitelisted sites, grab the first data, and concatenate them
59data_arr = data_arr.loc[selected]
60
61# Drop rows with missing winds
62data_arr = data_arr.dropna(how='any', subset=['wind_dir', 'wind_speed'])
63
64# First, look at the names of variables that the layout is expecting:
65simple_layout.names()
66
67###########################################
68# Next grab the simple variables out of the data we have (attaching correct units), and
69# put them into a dictionary that we will hand the plotting function later:
70
71# This is our container for the data
72data = {}
73
74# Copy out to stage everything together. In an ideal world, this would happen on
75# the data reading side of things, but we're not there yet.
76data['longitude'] = data_arr['lon'].values
77data['latitude'] = data_arr['lat'].values
78data['air_temperature'] = data_arr['air_temperature'].values * units.degC
79data['dew_point_temperature'] = data_arr['dew_point_temperature'].values * units.degC
80data['air_pressure_at_sea_level'] = data_arr['slp'].values * units('mbar')
81
82###########################################
83# Notice that the names (the keys) in the dictionary are the same as those that the
84# layout is expecting.
85#
86# Now perform a few conversions:
87#
88# - Get wind components from speed and direction
89# - Convert cloud fraction values to integer codes [0 - 8]
90# - Map METAR weather codes to WMO codes for weather symbols
91
92# Get the wind components, converting from m/s to knots as will be appropriate
93# for the station plot
94u, v = wind_components(data_arr['wind_speed'].values * units('m/s'),
95                       data_arr['wind_dir'].values * units.degree)
96data['eastward_wind'], data['northward_wind'] = u, v
97
98# Convert the fraction value into a code of 0-8, which can be used to pull out
99# the appropriate symbol
100data['cloud_coverage'] = (8 * data_arr['cloud_fraction']).fillna(10).values.astype(int)
101
102# Map weather strings to WMO codes, which we can use to convert to symbols
103# Only use the first symbol if there are multiple
104wx_text = data_arr['weather'].fillna('')
105data['current_wx1_symbol'] = [wx_code_map[s.split()[0] if ' ' in s else s] for s in wx_text]
106
107###########################################
108# All the data wrangling is finished, just need to set up plotting and go:
109# Set up the map projection and set up a cartopy feature for state borders
110proj = ccrs.LambertConformal(central_longitude=-95, central_latitude=35,
111                             standard_parallels=[35])
112
113###########################################
114# The payoff
115# ----------
116
117# Change the DPI of the resulting figure. Higher DPI drastically improves the
118# look of the text rendering
119plt.rcParams['savefig.dpi'] = 255
120
121# Create the figure and an axes set to the projection
122fig = plt.figure(figsize=(20, 10))
123add_metpy_logo(fig, 1080, 290, size='large')
124ax = fig.add_subplot(1, 1, 1, projection=proj)
125
126# Add some various map elements to the plot to make it recognizable
127ax.add_feature(cfeature.LAND)
128ax.add_feature(cfeature.OCEAN)
129ax.add_feature(cfeature.LAKES)
130ax.add_feature(cfeature.COASTLINE)
131ax.add_feature(cfeature.STATES)
132ax.add_feature(cfeature.BORDERS, linewidth=2)
133
134# Set plot bounds
135ax.set_extent((-118, -73, 23, 50))
136
137#
138# Here's the actual station plot
139#
140
141# Start the station plot by specifying the axes to draw on, as well as the
142# lon/lat of the stations (with transform). We also the fontsize to 12 pt.
143stationplot = StationPlot(ax, data['longitude'], data['latitude'],
144                          transform=ccrs.PlateCarree(), fontsize=12)
145
146# The layout knows where everything should go, and things are standardized using
147# the names of variables. So the layout pulls arrays out of `data` and plots them
148# using `stationplot`.
149simple_layout.plot(stationplot, data)
150
151plt.show()
152
153###########################################
154# or instead, a custom layout can be used:
155
156# Just winds, temps, and dewpoint, with colors. Dewpoint and temp will be plotted
157# out to Fahrenheit tenths. Extra data will be ignored
158custom_layout = StationPlotLayout()
159custom_layout.add_barb('eastward_wind', 'northward_wind', units='knots')
160custom_layout.add_value('NW', 'air_temperature', fmt='.1f', units='degF', color='darkred')
161custom_layout.add_value('SW', 'dew_point_temperature', fmt='.1f', units='degF',
162                        color='darkgreen')
163
164# Also, we'll add a field that we don't have in our dataset. This will be ignored
165custom_layout.add_value('E', 'precipitation', fmt='0.2f', units='inch', color='blue')
166
167# Create the figure and an axes set to the projection
168fig = plt.figure(figsize=(20, 10))
169add_metpy_logo(fig, 1080, 290, size='large')
170ax = fig.add_subplot(1, 1, 1, projection=proj)
171
172# Add some various map elements to the plot to make it recognizable
173ax.add_feature(cfeature.LAND)
174ax.add_feature(cfeature.OCEAN)
175ax.add_feature(cfeature.LAKES)
176ax.add_feature(cfeature.COASTLINE)
177ax.add_feature(cfeature.STATES)
178ax.add_feature(cfeature.BORDERS, linewidth=2)
179
180# Set plot bounds
181ax.set_extent((-118, -73, 23, 50))
182
183#
184# Here's the actual station plot
185#
186
187# Start the station plot by specifying the axes to draw on, as well as the
188# lon/lat of the stations (with transform). We also the fontsize to 12 pt.
189stationplot = StationPlot(ax, data['longitude'], data['latitude'],
190                          transform=ccrs.PlateCarree(), fontsize=12)
191
192# The layout knows where everything should go, and things are standardized using
193# the names of variables. So the layout pulls arrays out of `data` and plots them
194# using `stationplot`.
195custom_layout.plot(stationplot, data)
196
197plt.show()
198