|
1
|
|
|
# Copyright (c) 2008-2017 MetPy Developers. |
|
2
|
|
|
# Distributed under the terms of the BSD 3-Clause License. |
|
3
|
|
|
# SPDX-License-Identifier: BSD-3-Clause |
|
4
|
|
|
""" |
|
5
|
|
|
Isentropic Cross-section |
|
6
|
|
|
======================== |
|
7
|
|
|
|
|
8
|
|
|
Make an isentropic cross section with MetPy |
|
9
|
|
|
|
|
10
|
|
|
With the function `metpy.calc.extract_cross_section`, data can be extracted from 3-dimensional |
|
11
|
|
|
input along a specified line and plotted in a 2-dimensional cross-section. |
|
12
|
|
|
""" |
|
13
|
|
|
|
|
14
|
|
|
################################ |
|
15
|
|
|
import cartopy.crs as ccrs |
|
16
|
|
|
import cartopy.feature as cfeature |
|
17
|
|
|
import matplotlib.pyplot as plt |
|
18
|
|
|
from netCDF4 import Dataset, num2date |
|
19
|
|
|
import numpy as np |
|
20
|
|
|
|
|
21
|
|
|
import metpy.calc as mcalc |
|
22
|
|
|
from metpy.cbook import get_test_data |
|
23
|
|
|
from metpy.units import units |
|
24
|
|
|
|
|
25
|
|
|
####################################### |
|
26
|
|
|
# **Getting the data** |
|
27
|
|
|
# |
|
28
|
|
|
# In this example, NARR reanalysis data for 18 UTC 04 April 1987 from the National Centers |
|
29
|
|
|
# for Environmental Information (https://nomads.ncdc.noaa.gov) will be used. |
|
30
|
|
|
|
|
31
|
|
|
data = Dataset(get_test_data('narr_example.nc', False)) |
|
32
|
|
|
|
|
33
|
|
|
########################## |
|
34
|
|
|
print(list(data.variables)) |
|
35
|
|
|
|
|
36
|
|
|
############################# |
|
37
|
|
|
# We will reduce the dimensionality of the data as it is pulled in to remove an empty time |
|
38
|
|
|
# dimension. Additionally, units are required for input data, so the proper units will also |
|
39
|
|
|
# be attached. |
|
40
|
|
|
|
|
41
|
|
|
|
|
42
|
|
|
# Assign data to variable names |
|
43
|
|
|
dtime = data.variables['Geopotential_height'].dimensions[0] |
|
44
|
|
|
dlev = data.variables['Geopotential_height'].dimensions[1] |
|
45
|
|
|
lat = data.variables['lat'][:] |
|
46
|
|
|
lon = data.variables['lon'][:] |
|
47
|
|
|
lev = data.variables[dlev][:] * units(data.variables[dlev].units) |
|
48
|
|
|
times = data.variables[dtime] |
|
49
|
|
|
vtimes = num2date(times[:], times.units) |
|
50
|
|
|
|
|
51
|
|
|
temps = data.variables['Temperature'] |
|
52
|
|
|
tmp = temps[0, :] * units.kelvin |
|
53
|
|
|
uwnd = data.variables['u_wind'][0, :] * units(data.variables['u_wind'].units) |
|
54
|
|
|
vwnd = data.variables['v_wind'][0, :] * units(data.variables['v_wind'].units) |
|
55
|
|
|
hgt = data.variables['Geopotential_height'][0, :] * units.meter |
|
56
|
|
|
spech = (data.variables['Specific_humidity'][0, :] * |
|
57
|
|
|
units(data.variables['Specific_humidity'].units)) |
|
58
|
|
|
|
|
59
|
|
|
|
|
60
|
|
|
######################### |
|
61
|
|
|
# Create the coordinate system and projection for plotting the inset map |
|
62
|
|
|
crs = ccrs.LambertConformal(central_longitude=-100.0, central_latitude=45.0) |
|
63
|
|
|
tlatlons = crs.transform_points(ccrs.PlateCarree(), lon, lat) |
|
64
|
|
|
tlons = tlatlons[:, :, 0] |
|
65
|
|
|
tlats = tlatlons[:, :, 1] |
|
66
|
|
|
|
|
67
|
|
|
########################## |
|
68
|
|
|
# Get data to plot state and province boundaries |
|
69
|
|
|
states_provinces = cfeature.NaturalEarthFeature(category='cultural', |
|
70
|
|
|
name='admin_1_states_provinces_lakes', |
|
71
|
|
|
scale='50m', |
|
72
|
|
|
facecolor='none') |
|
73
|
|
|
|
|
74
|
|
|
########################## |
|
75
|
|
|
# Define a 2-D lat/lon index function to find the nearest grid point to a specified coordinate |
|
76
|
|
|
|
|
77
|
|
|
|
|
78
|
|
|
def lat_lon_2d_index(y, x, lat1, lon1): |
|
79
|
|
|
"""Calculate the index values of the grid point nearest a given lat/lon point. |
|
80
|
|
|
|
|
81
|
|
|
This function calculates the distance from a desired lat/lon point |
|
82
|
|
|
to each element of a 2D array of lat/lon values, typically from model output, |
|
83
|
|
|
and determines the index value corresponding to the nearest lat/lon grid point. |
|
84
|
|
|
|
|
85
|
|
|
y = latitude array |
|
86
|
|
|
x = longitude array |
|
87
|
|
|
lon1 = longitude point (signle value) |
|
88
|
|
|
lat1 = latitude point (single value) |
|
89
|
|
|
|
|
90
|
|
|
Returns the index value for nearest lat/lon point on grid |
|
91
|
|
|
|
|
92
|
|
|
Equations for variable distiance between longitudes from |
|
93
|
|
|
http://andrew.hedges.name/experiments/haversine/, code by Kevin Goebbert. |
|
94
|
|
|
""" |
|
95
|
|
|
r = 6373. * 1000. # Earth's Radius in meters |
|
96
|
|
|
rad = np.pi / 180. |
|
97
|
|
|
x1 = np.ones(x.shape) * lon1 |
|
98
|
|
|
y1 = np.ones(y.shape) * lat1 |
|
99
|
|
|
dlon = np.abs(x - x1) |
|
100
|
|
|
dlat = np.abs(y - y1) |
|
101
|
|
|
a = (np.sin(rad * dlat / 2.))**2 + (np.cos(rad * y1) * np.cos(rad * y) * |
|
102
|
|
|
(np.sin(rad * dlon / 2.))**2) |
|
103
|
|
|
c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a)) |
|
104
|
|
|
d = r * c |
|
105
|
|
|
return np.unravel_index(d.argmin(), d.shape) |
|
106
|
|
|
|
|
107
|
|
|
|
|
108
|
|
|
############################## |
|
109
|
|
|
# **Insentropic Cross Section** |
|
110
|
|
|
# |
|
111
|
|
|
# With the projection set up, we will define the endpoints, extract the |
|
112
|
|
|
# cross-sectional data, and plot with an inset map. |
|
113
|
|
|
|
|
114
|
|
|
|
|
115
|
|
|
def cross_section(isentlev, num, left_lat, left_lon, right_lat, right_lon): |
|
116
|
|
|
"""Plot an isentropic cross-section.""" |
|
117
|
|
|
# get the coordinates of the endpoints for the cross-section |
|
118
|
|
|
left_coord = np.array((float(left_lat), float(left_lon))) |
|
119
|
|
|
right_coord = np.array((float(right_lat), float(right_lon))) |
|
120
|
|
|
|
|
121
|
|
|
# Calculate data for the inset isentropic map |
|
122
|
|
|
isent_anal = mcalc.isentropic_interpolation(float(isentlev) * units.kelvin, lev, tmp, |
|
123
|
|
|
spech, tmpk_out=True) |
|
124
|
|
|
isentprs = isent_anal[0] |
|
125
|
|
|
isenttmp = isent_anal[1] |
|
126
|
|
|
isentspech = isent_anal[2] |
|
127
|
|
|
isentrh = mcalc.relative_humidity_from_specific_humidity(isentspech, isenttmp, isentprs) |
|
128
|
|
|
|
|
129
|
|
|
# Find index values for the cross section slice |
|
130
|
|
|
iright = lat_lon_2d_index(lat, lon, right_coord[0], right_coord[1]) |
|
131
|
|
|
ileft = lat_lon_2d_index(lat, lon, left_coord[0], left_coord[1]) |
|
132
|
|
|
|
|
133
|
|
|
# Get the cross-section slice data |
|
134
|
|
|
cross_data = mcalc.extract_cross_section(ileft, iright, lat, lon, tmp, uwnd, vwnd, spech, |
|
135
|
|
|
num=num) |
|
136
|
|
|
cross_lat = cross_data[0] |
|
137
|
|
|
cross_lon = cross_data[1] |
|
138
|
|
|
cross_t = cross_data[2] |
|
139
|
|
|
cross_u = cross_data[3] |
|
140
|
|
|
cross_v = cross_data[4] |
|
141
|
|
|
cross_spech = cross_data[5] |
|
142
|
|
|
|
|
143
|
|
|
# Calculate theta and RH on the cross-section |
|
144
|
|
|
cross_theta = mcalc.potential_temperature(lev[:, np.newaxis], cross_t) |
|
145
|
|
|
cross_rh = mcalc.relative_humidity_from_specific_humidity(cross_spech, cross_t, |
|
146
|
|
|
lev[:, np.newaxis]) |
|
147
|
|
|
|
|
148
|
|
|
# Create figure for ploting |
|
149
|
|
|
fig = plt.figure(1, figsize=(17., 12.)) |
|
150
|
|
|
|
|
151
|
|
|
# Plot the cross section |
|
152
|
|
|
ax1 = plt.axes() |
|
153
|
|
|
ax1.set_yscale('symlog') |
|
154
|
|
|
ax1.grid() |
|
155
|
|
|
cint = np.arange(250, 450, 5) |
|
156
|
|
|
|
|
157
|
|
|
# Determine whether to label x-axis with lat or lon values |
|
158
|
|
|
if np.abs(left_lon - right_lon) > np.abs(left_lat - right_lat): |
|
159
|
|
|
cs = ax1.contour(cross_lon, lev[::-1], cross_theta[::-1, :], cint, colors='tab:red') |
|
160
|
|
|
cf = ax1.contourf(cross_lon, lev[::-1], cross_rh[::-1, :], range(10, 106, 5), |
|
161
|
|
|
cmap=plt.cm.gist_earth_r) |
|
162
|
|
|
ax1.barbs(cross_lon[4::4], lev, cross_u[:, 4::4], cross_v[:, 4::4], length=6) |
|
163
|
|
|
plt.xlabel('Longitude (Degrees East)') |
|
164
|
|
|
else: |
|
165
|
|
|
cs = ax1.contour(cross_lat[::-1], lev[::-1], cross_theta[::-1, ::-1], cint, |
|
166
|
|
|
colors='tab:red') |
|
167
|
|
|
cf = ax1.contourf(cross_lat[::-1], lev[::-1], cross_rh[::-1, ::-1], range(10, 106, 5), |
|
168
|
|
|
cmap=plt.cm.gist_earth_r) |
|
169
|
|
|
ax1.barbs(cross_lat[::-4], lev, cross_u[:, ::-4], cross_v[:, ::-4], length=6) |
|
170
|
|
|
plt.xlim(cross_lat[0], cross_lat[-1]) |
|
171
|
|
|
plt.xlabel('Latitude (Degrees North)') |
|
172
|
|
|
|
|
173
|
|
|
# Label the cross section axes |
|
174
|
|
|
plt.clabel(cs, fontsize=10, inline=1, inline_spacing=7, |
|
175
|
|
|
fmt='%i', rightside_up=True, use_clabeltext=True) |
|
176
|
|
|
cb = plt.colorbar(cf, orientation='horizontal', extend=max, aspect=65, shrink=0.75, |
|
177
|
|
|
pad=0.06, extendrect='True') |
|
178
|
|
|
cb.set_label('Relative Humidity', size='x-large') |
|
179
|
|
|
plt.ylabel('Pressure (hPa)') |
|
180
|
|
|
ax1.set_yticklabels(np.arange(1000, 50, -50)) |
|
181
|
|
|
plt.ylim(lev[0], lev[-1]) |
|
182
|
|
|
plt.yticks(np.arange(1000, 50, -50)) |
|
183
|
|
|
|
|
184
|
|
|
# Add a title |
|
185
|
|
|
plt.title(('NARR Isentropic Cross-Section: ' + str(left_coord[0]) + ' N, ' + |
|
186
|
|
|
str(left_coord[1]) + ' E to ' + str(right_coord[0]) + ' N, ' + |
|
187
|
|
|
str(right_coord[1]) + ' E'), loc='left') |
|
188
|
|
|
plt.title('VALID: {:s}'.format(str(vtimes[0])), loc='right') |
|
189
|
|
|
|
|
190
|
|
|
# Add Inset Map |
|
191
|
|
|
ax2 = fig.add_axes([0.125, 0.643, 0.25, 0.25], projection=crs) |
|
192
|
|
|
|
|
193
|
|
|
# Coordinates to limit map area |
|
194
|
|
|
bounds = [(-122., -75., 25., 50.)] |
|
195
|
|
|
|
|
196
|
|
|
# Convert endpoints of cross-section line |
|
197
|
|
|
left = crs.transform_point(left_coord[1], left_coord[0], ccrs.PlateCarree()) |
|
198
|
|
|
right = crs.transform_point(right_coord[1], right_coord[0], ccrs.PlateCarree()) |
|
199
|
|
|
|
|
200
|
|
|
# Limit extent of inset map |
|
201
|
|
|
ax2.set_extent(*bounds, crs=ccrs.PlateCarree()) |
|
202
|
|
|
ax2.coastlines('50m', edgecolor='black', linewidth=0.75) |
|
203
|
|
|
ax2.add_feature(states_provinces, edgecolor='black', linewidth=0.5) |
|
204
|
|
|
|
|
205
|
|
|
# Plot the surface |
|
206
|
|
|
clevisent = np.arange(0, 1000, 25) |
|
207
|
|
|
cs = ax2.contour(tlons, tlats, isentprs[0, :, :], clevisent, |
|
208
|
|
|
colors='k', linewidths=1.0, linestyles='solid') |
|
209
|
|
|
plt.clabel(cs, fontsize=10, inline=1, inline_spacing=7, |
|
210
|
|
|
fmt='%i', rightside_up=True, use_clabeltext=True) |
|
211
|
|
|
|
|
212
|
|
|
# Plot RH |
|
213
|
|
|
cf = ax2.contourf(tlons, tlats, isentrh[0, :, :], range(10, 106, 5), |
|
214
|
|
|
cmap=plt.cm.gist_earth_r) |
|
215
|
|
|
|
|
216
|
|
|
# Plot the cross section line |
|
217
|
|
|
plt.plot([left[0], right[0]], [left[1], right[1]], color='r') |
|
218
|
|
|
plt.show() |
|
219
|
|
|
|
|
220
|
|
|
|
|
221
|
|
|
# Finally, call the plotting function and show the map |
|
222
|
|
|
cross_section(296., 80, 37., -95., 35.5, -65.) |
|
223
|
|
|
|
Unidata /
MetPy
