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