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"""
5Wind and Sea Level Pressure Interpolation
6=========================================
7
8Interpolate sea level pressure, as well as wind component data,
9to make a consistent looking analysis, featuring contours of pressure and wind barbs.
10"""
11import cartopy.crs as ccrs
12import cartopy.feature as cfeature
13from matplotlib.colors import BoundaryNorm
14import matplotlib.pyplot as plt
15import numpy as np
16import pandas as pd
17
18from metpy.calc import wind_components
19from metpy.cbook import get_test_data
20from metpy.interpolate import interpolate_to_grid, remove_nan_observations
21from metpy.plots import add_metpy_logo
22from metpy.units import units
23
24to_proj = ccrs.AlbersEqualArea(central_longitude=-97., central_latitude=38.)
25
26###########################################
27# Read in data
28with get_test_data('station_data.txt') as f:
29    data = pd.read_csv(f, header=0, usecols=(2, 3, 4, 5, 18, 19),
30                       names=['latitude', 'longitude', 'slp', 'temperature', 'wind_dir',
31                              'wind_speed'],
32                       na_values=-99999)
33
34###########################################
35# Project the lon/lat locations to our final projection
36lon = data['longitude'].values
37lat = data['latitude'].values
38xp, yp, _ = to_proj.transform_points(ccrs.Geodetic(), lon, lat).T
39
40###########################################
41# Remove all missing data from pressure
42x_masked, y_masked, pressure = remove_nan_observations(xp, yp, data['slp'].values)
43
44###########################################
45# Interpolate pressure using Cressman interpolation
46slpgridx, slpgridy, slp = interpolate_to_grid(x_masked, y_masked, pressure,
47                                              interp_type='cressman', minimum_neighbors=1,
48                                              search_radius=400000, hres=100000)
49
50##########################################
51# Get wind information and mask where either speed or direction is unavailable
52wind_speed = (data['wind_speed'].values * units('m/s')).to('knots')
53wind_dir = data['wind_dir'].values * units.degree
54
55good_indices = np.where((~np.isnan(wind_dir)) & (~np.isnan(wind_speed)))
56
57x_masked = xp[good_indices]
58y_masked = yp[good_indices]
59wind_speed = wind_speed[good_indices]
60wind_dir = wind_dir[good_indices]
61
62###########################################
63# Calculate u and v components of wind and then interpolate both.
64#
65# Both will have the same underlying grid so throw away grid returned from v interpolation.
66u, v = wind_components(wind_speed, wind_dir)
67
68windgridx, windgridy, uwind = interpolate_to_grid(x_masked, y_masked, np.array(u),
69                                                  interp_type='cressman', search_radius=400000,
70                                                  hres=100000)
71
72_, _, vwind = interpolate_to_grid(x_masked, y_masked, np.array(v), interp_type='cressman',
73                                  search_radius=400000, hres=100000)
74
75###########################################
76# Get temperature information
77x_masked, y_masked, t = remove_nan_observations(xp, yp, data['temperature'].values)
78tempx, tempy, temp = interpolate_to_grid(x_masked, y_masked, t, interp_type='cressman',
79                                         minimum_neighbors=3, search_radius=400000, hres=35000)
80
81temp = np.ma.masked_where(np.isnan(temp), temp)
82
83###########################################
84# Set up the map and plot the interpolated grids appropriately.
85levels = list(range(-20, 20, 1))
86cmap = plt.get_cmap('viridis')
87norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)
88
89fig = plt.figure(figsize=(20, 10))
90add_metpy_logo(fig, 360, 120, size='large')
91view = fig.add_subplot(1, 1, 1, projection=to_proj)
92
93view.set_extent([-120, -70, 20, 50])
94view.add_feature(cfeature.STATES.with_scale('50m'))
95view.add_feature(cfeature.OCEAN)
96view.add_feature(cfeature.COASTLINE.with_scale('50m'))
97view.add_feature(cfeature.BORDERS, linestyle=':')
98
99cs = view.contour(slpgridx, slpgridy, slp, colors='k', levels=list(range(990, 1034, 4)))
100view.clabel(cs, inline=1, fontsize=12, fmt='%i')
101
102mmb = view.pcolormesh(tempx, tempy, temp, cmap=cmap, norm=norm)
103fig.colorbar(mmb, shrink=.4, pad=0.02, boundaries=levels)
104
105view.barbs(windgridx, windgridy, uwind, vwind, alpha=.4, length=5)
106
107view.set_title('Surface Temperature (shaded), SLP, and Wind.')
108
109plt.show()
110