|
1
|
|
|
# Copyright (c) 2008-2016 MetPy Developers. |
|
2
|
|
|
# Distributed under the terms of the BSD 3-Clause License. |
|
3
|
|
|
# SPDX-License-Identifier: BSD-3-Clause |
|
4
|
|
|
""" |
|
5
|
|
|
Four Panel Map |
|
6
|
|
|
=============== |
|
7
|
|
|
|
|
8
|
|
|
By reading model output data from a netCDF file, we can create a four panel plot showing: |
|
9
|
|
|
|
|
10
|
|
|
* 300 hPa heights and winds |
|
11
|
|
|
* 500 hPa heights and absolute vorticity |
|
12
|
|
|
* Surface temperatures |
|
13
|
|
|
* Precipitable water |
|
14
|
|
|
""" |
|
15
|
|
|
|
|
16
|
|
|
########################################### |
|
17
|
|
|
|
|
18
|
|
|
import netCDF4 |
|
19
|
|
|
import numpy as np |
|
20
|
|
|
import cartopy.crs as ccrs |
|
21
|
|
|
import cartopy.feature as cfeature |
|
22
|
|
|
import matplotlib.gridspec as gridspec |
|
23
|
|
|
import matplotlib.pyplot as plt |
|
24
|
|
|
import scipy.ndimage as ndimage |
|
25
|
|
|
from metpy.cbook import get_test_data |
|
26
|
|
|
from metpy.units import units |
|
27
|
|
|
|
|
28
|
|
|
########################################### |
|
29
|
|
|
|
|
30
|
|
|
|
|
31
|
|
|
# Function used to create the map subplots |
|
32
|
|
|
def plot_background(ax): |
|
33
|
|
|
ax.set_extent([235., 290., 20., 55.]) |
|
34
|
|
|
ax.coastlines('50m', edgecolor='black', linewidth=0.5) |
|
35
|
|
|
ax.add_feature(states_provinces, edgecolor='black', linewidth=0.5) |
|
36
|
|
|
ax.add_feature(country_borders, edgecolor='black', linewidth=0.5) |
|
37
|
|
|
return ax |
|
38
|
|
|
|
|
39
|
|
|
|
|
40
|
|
|
########################################### |
|
41
|
|
|
|
|
42
|
|
|
# Open the example netCDF data |
|
43
|
|
|
ds = netCDF4.Dataset(get_test_data('gfs_output.nc', False)) |
|
44
|
|
|
print(ds) |
|
45
|
|
|
|
|
46
|
|
|
########################################### |
|
47
|
|
|
|
|
48
|
|
|
# Convert number of hours since the reference time into an actual date |
|
49
|
|
|
time_vals = netCDF4.num2date(ds.variables['time'][:].squeeze(), ds.variables['time'].units) |
|
50
|
|
|
|
|
51
|
|
|
########################################### |
|
52
|
|
|
|
|
53
|
|
|
# Combine 1D latitude and longitudes into a 2D grid of locations |
|
54
|
|
|
lon_2d, lat_2d = np.meshgrid(ds.variables['lon'][:], ds.variables['lat'][:]) |
|
55
|
|
|
|
|
56
|
|
|
########################################### |
|
57
|
|
|
|
|
58
|
|
|
# Assign units |
|
59
|
|
|
vort_500 = ds.variables['vort_500'][0] * units('1/s') # To plot a e-5 s^-1 |
|
60
|
|
|
surface_temp = ds.variables['temp'][0] * units.kelvin |
|
61
|
|
|
precip_water = ds.variables['precip_water'][0] * units.mm |
|
62
|
|
|
winds_300 = ds.variables['winds_300'][0] * units('m/s') |
|
63
|
|
|
|
|
64
|
|
|
########################################### |
|
65
|
|
|
|
|
66
|
|
|
# Do unit conversions to what we wish to plot |
|
67
|
|
|
vort_500 = vort_500 * 1e5 |
|
68
|
|
|
surface_temp = surface_temp.to('degF') |
|
69
|
|
|
precip_water = precip_water.to('inches') |
|
70
|
|
|
winds_300 = winds_300.to('knots') |
|
71
|
|
|
|
|
72
|
|
|
########################################### |
|
73
|
|
|
|
|
74
|
|
|
# Smooth the height data |
|
75
|
|
|
heights_300 = ndimage.gaussian_filter(ds.variables['heights_300'][0], sigma=1.5, order=0) |
|
76
|
|
|
heights_500 = ndimage.gaussian_filter(ds.variables['heights_500'][0], sigma=1.5, order=0) |
|
77
|
|
|
|
|
78
|
|
|
########################################### |
|
79
|
|
|
|
|
80
|
|
|
# Add state boundaries to plot |
|
81
|
|
|
states_provinces = cfeature.NaturalEarthFeature(category='cultural', |
|
82
|
|
|
name='admin_1_states_provinces_lines', |
|
83
|
|
|
scale='50m', facecolor='none') |
|
84
|
|
|
|
|
85
|
|
|
# Add country borders to plot |
|
86
|
|
|
country_borders = cfeature.NaturalEarthFeature(category='cultural', |
|
87
|
|
|
name='admin_0_countries', |
|
88
|
|
|
scale='50m', facecolor='none') |
|
89
|
|
|
|
|
90
|
|
|
crs = ccrs.LambertConformal(central_longitude=-100.0, central_latitude=45.0) |
|
91
|
|
|
|
|
92
|
|
|
########################################### |
|
93
|
|
|
|
|
94
|
|
|
# Create the figure |
|
95
|
|
|
fig = plt.figure(figsize=(20, 15)) |
|
96
|
|
|
gs = gridspec.GridSpec(5, 2, height_ratios=[1, .05, 1, .05, 0], bottom=.05, top=.95, wspace=.1) |
|
97
|
|
|
|
|
98
|
|
|
|
|
99
|
|
|
# Upper left plot - 300-hPa winds and geopotential heights |
|
100
|
|
|
ax1 = plt.subplot(gs[0, 0], projection=crs) |
|
101
|
|
|
plot_background(ax1) |
|
102
|
|
|
cf1 = ax1.contourf(lon_2d, lat_2d, winds_300, cmap='cool', transform=ccrs.PlateCarree()) |
|
103
|
|
|
c1 = ax1.contour(lon_2d, lat_2d, heights_300, colors='black', linewidth=2, |
|
104
|
|
|
transform=ccrs.PlateCarree()) |
|
105
|
|
|
plt.clabel(c1, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True) |
|
106
|
|
|
|
|
107
|
|
|
ax2 = plt.subplot(gs[1, 0]) |
|
108
|
|
|
cb1 = plt.colorbar(cf1, cax=ax2, orientation='horizontal') |
|
109
|
|
|
cb1.set_label('knots', size='x-large') |
|
110
|
|
|
ax1.set_title('300-hPa Wind Speeds and Heights', fontsize=16) |
|
111
|
|
|
|
|
112
|
|
|
# Upper right plot - 500mb absolute vorticity and geopotential heights |
|
113
|
|
|
ax3 = plt.subplot(gs[0, 1], projection=crs) |
|
114
|
|
|
plot_background(ax3) |
|
115
|
|
|
cf2 = ax3.contourf(lon_2d, lat_2d, vort_500, cmap='BrBG', transform=ccrs.PlateCarree(), |
|
116
|
|
|
zorder=0, norm=plt.Normalize(-32, 32), latlon=True) |
|
117
|
|
|
c2 = ax3.contour(lon_2d, lat_2d, heights_500, colors='k', lw=2, transform=ccrs.PlateCarree()) |
|
118
|
|
|
plt.clabel(c2, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True) |
|
119
|
|
|
|
|
120
|
|
|
ax4 = plt.subplot(gs[1, 1]) |
|
121
|
|
|
cb2 = plt.colorbar(cf2, cax=ax4, orientation='horizontal') |
|
122
|
|
|
cb2.set_label(r'$10^{-5}$ s$^{-1}$', size='x-large') |
|
123
|
|
|
ax3.set_title('500-hPa Absolute Vorticity and Heights', fontsize=16) |
|
124
|
|
|
|
|
125
|
|
|
# Lower left plot - surface temperatures |
|
126
|
|
|
ax5 = plt.subplot(gs[2, 0], projection=crs) |
|
127
|
|
|
plot_background(ax5) |
|
128
|
|
|
cf3 = ax5.contourf(lon_2d, lat_2d, surface_temp, cmap='YlOrRd', |
|
129
|
|
|
transform=ccrs.PlateCarree(), zorder=0) |
|
130
|
|
|
|
|
131
|
|
|
ax6 = plt.subplot(gs[3, 0]) |
|
132
|
|
|
cb3 = plt.colorbar(cf3, cax=ax6, orientation='horizontal') |
|
133
|
|
|
cb3.set_label(r'$^o$F', size='x-large') |
|
134
|
|
|
ax5.set_title('Surface Temperatures', fontsize=16) |
|
135
|
|
|
|
|
136
|
|
|
# Lower right plot - precipitable water entire atmosphere |
|
137
|
|
|
ax7 = plt.subplot(gs[2, 1], projection=crs) |
|
138
|
|
|
plot_background(ax7) |
|
139
|
|
|
cf4 = plt.contourf(lon_2d, lat_2d, precip_water, cmap='Greens', |
|
140
|
|
|
transform=ccrs.PlateCarree(), zorder=0) |
|
141
|
|
|
|
|
142
|
|
|
ax8 = plt.subplot(gs[3, 1]) |
|
143
|
|
|
cb4 = plt.colorbar(cf4, cax=ax8, orientation='horizontal') |
|
144
|
|
|
cb4.set_label('in.', size='x-large') |
|
145
|
|
|
ax7.set_title('Precipitable Water', fontsize=16) |
|
146
|
|
|
|
|
147
|
|
|
fig.suptitle('{0:%d %B %Y %H:%MZ}'.format(time_vals), fontsize=24) |
|
148
|
|
|
|
|
149
|
|
|
# Display the plot |
|
150
|
|
|
plt.show() |
|
151
|
|
|
|
Unidata /
MetPy
