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