|
1.Unidata Python Gallery 基礎庫安裝到位，運行腳本即可出圖 腳本可點擊文末閱讀原文下載 頁面鏈接：https://unidata./python-gallery/examples/index.html 
2.Useful Python Tools
This is a list of useful and/or new Python tools that the Unidata Python Team and community are keeping an eye on or using. Unidata ProjectsMetPy - General meteorological toolkit siphon - Remote data access netCDF4-python - netCDF4 API
Meteorology SpecificData WranglingPlottingCore Python/InterfaceiPython - Interactive Python shell Jupyter - Notebooks and the new Jupyter Lab pathlib - Easy file path manipulation
EducationPerformance3.實例 
1from datetime import datetime
2
3import matplotlib.pyplot as plt
4import metpy.calc as mpcalc
5from metpy.units import units
6import numpy as np
7from pyproj import Geod
8from scipy.interpolate import griddata
9from scipy.ndimage import gaussian_filter
10from siphon.simplewebservice.wyoming import WyomingUpperAir
11
12def vertical_interpolate(vcoord_data, interp_var, interp_levels):
13 '''A function to interpolate sounding data from each station to
14 every millibar. Assumes a log-linear relationship.
15
16 Input
17 -----
18 vcoord_data : A 1D array of vertical level values (e.g., pressure from a radiosonde)
19 interp_var : A 1D array of the variable to be interpolated to all pressure levels
20 vcoord_interp_levels : A 1D array containing veritcal levels to interpolate to
21
22 Return
23 ------
24 interp_data : A 1D array that contains the interpolated variable on the interp_levels
25 '''
26
27 # Make veritcal coordinate data and grid level log variables
28 lnp = np.log(vcoord_data)
29 lnp_intervals = np.log(interp_levels)
30
31 # Use numpy to interpolate from observed levels to grid levels
32 interp_data = np.interp(lnp_intervals[::-1], lnp[::-1], interp_var[::-1])[::-1]
33
34 # Mask for missing data (generally only near the surface)
35 mask_low = interp_levels > vcoord_data[0]
36 mask_high = interp_levels < vcoord_data[-1]
37 interp_data[mask_low] = interp_var[0]
38 interp_data[mask_high] = interp_var[-1]
39
40 return interp_data
41
42def radisonde_cross_section(stns, data, start=1000, end=100, step=10):
43 '''This function takes a list of radiosonde observation sites with a
44 dictionary of Pandas Dataframes with the requesite data for each station.
45
46 Input
47 -----
48 stns : List of statition three-letter identifiers
49 data : A dictionary of Pandas Dataframes containing the radiosonde observations
50 for the stations
51 start : interpolation start value, optional (default = 1000 hPa)
52 end : Interpolation end value, optional (default = 100 hPa)
53 step : Interpolation interval, option (default = 10 hPa)
54
55 Return
56 ------
57 cross_section : A dictionary that contains the following variables
58
59 grid_data : An interpolated grid with 100 points between the first and last station,
60 with the corresponding number of vertical points based on start, end, and interval
61 (default is 90)
62 obs_distance : An array of distances between each radiosonde observation location
63 x_grid : A 2D array of horizontal direction grid points
64 p_grid : A 2D array of vertical pressure levels
65 ground_elevation : A representation of the terrain between radiosonde observation sites
66 based on the elevation of each station converted to pressure using the standard
67 atmosphere
68
69 '''
70 # Set up vertical grid, largest value first (high pressure nearest surface)
71 vertical_levels = np.arange(start, end-1, -step)
72
73 # Number of vertical levels and stations
74 plevs = len(vertical_levels)
75 nstns = len(stns)
76
77 # Create dictionary of interpolated values and include neccsary attribute data
78 # including lat, lon, and elevation of each station
79 lats = []
80 lons = []
81 elev = []
82 keys = data[stns[0]].keys()[:8]
83 tmp_grid = dict.fromkeys(keys)
84
85 # Interpolate all variables for each radiosonde observation
86 # Temperature, Dewpoint, U-wind, V-wind
87 for key in tmp_grid.keys():
88 tmp_grid[key] = np.empty((nstns, plevs))
89 for station, loc in zip(stns, range(nstns)):
90 if key == 'pressure':
91 lats.append(data[station].latitude[0])
92 lons.append(data[station].longitude[0])
93 elev.append(data[station].elevation[0])
94 tmp_grid[key][loc, :] = vertical_levels
95 else:
96 tmp_grid[key][loc, :] = vertical_interpolate(
97 data[station]['pressure'].values, data[station][key].values,
98 vertical_levels)
99
100 # Compute distance between each station using Pyproj
101 g = Geod(ellps='sphere')
102 _, _, dist = g.inv(nstns*[lons[0]], nstns*[lats[0]], lons[:], lats[:])
103
104 # Compute sudo ground elevation in pressure from standard atmsophere and the elevation
105 # of each station
106 ground_elevation = mpcalc.height_to_pressure_std(np.array(elev) * units('meters'))
107
108 # Set up grid for 2D interpolation
109 grid = dict.fromkeys(keys)
110 x = np.linspace(dist[0], dist[-1], 100)
111 nx = len(x)
112
113 pp, xx = np.meshgrid(vertical_levels, x)
114 pdist, ddist = np.meshgrid(vertical_levels, dist)
115
116 # Interpolate to 2D grid using scipy.interpolate.griddata
117 for key in grid.keys():
118 grid[key] = np.empty((nx, plevs))
119 grid[key][:] = griddata((ddist.flatten(), pdist.flatten()),
120 tmp_grid[key][:].flatten(),
121 (xx, pp),
122 method='cubic')
123
124 # Gather needed data in dictionary for return
125 cross_section = {'grid_data': grid, 'obs_distance': dist,
126 'x_grid': xx, 'p_grid': pp, 'elevation': ground_elevation}
127 return cross_section
128# A roughly east-west cross section
129stn_list = ['DNR', 'LBF', 'OAX', 'DVN', 'DTX', 'BUF']
130
131# Set a date and hour of your choosing
132date = datetime(2019, 12, 20, 0)
133
134df = {}
135
136# Loop over stations to get data and put into dictionary
137for station in stn_list:
138 df[station] = WyomingUpperAir.request_data(date, station)
139xsect = radisonde_cross_section(stn_list, df)
140
141potemp = mpcalc.potential_temperature(
142 xsect['p_grid'] * units('hPa'), xsect['grid_data']['temperature'] * units('degC'))
143
144relhum = mpcalc.relative_humidity_from_dewpoint(
145 xsect['grid_data']['temperature'] * units('degC'),
146 xsect['grid_data']['dewpoint'] * units('degC'))
147
148mixrat = mpcalc.mixing_ratio_from_relative_humidity(relhum,
149 xsect['grid_data']['temperature'] *
150 units('degC'),
151 xsect['p_grid'] * units('hPa'))
152
153fig = plt.figure(figsize=(17, 11))
154
155# Specify plotting axis (single panel)
156ax = plt.subplot(111)
157
158# Set y-scale to be log since pressure decreases exponentially with height
159ax.set_yscale('log')
160
161# Set limits, tickmarks, and ticklabels for y-axis
162ax.set_ylim([1030, 101])
163ax.set_yticks(range(1000, 101, -100))
164ax.set_yticklabels(range(1000, 101, -100))
165
166# Invert the y-axis since pressure decreases with increasing height
167ax.invert_yaxis()
168
169# Plot the sudo elevation on the cross section
170ax.fill_between(xsect['obs_distance'], xsect['elevation'].m, 1030,
171 where=xsect['elevation'].m <= 1030, facecolor='lightgrey',
172 interpolate=True, zorder=10)
173# Don't plot xticks
174plt.xticks([], [])
175
176# Plot wind barbs for each sounding location
177for stn, stn_name in zip(range(len(stn_list)), stn_list):
178 ax.axvline(xsect['obs_distance'][stn], ymin=0, ymax=1,
179 linewidth=2, color='blue', zorder=11)
180 ax.text(xsect['obs_distance'][stn], 1100, stn_name, ha='center', color='blue')
181 ax.barbs(xsect['obs_distance'][stn], df[stn_name]['pressure'][::2],
182 df[stn_name]['u_wind'][::2, None],
183 df[stn_name]['v_wind'][::2, None], zorder=15)
184
185# Plot smoothed potential temperature grid (K)
186cs = ax.contour(xsect['x_grid'], xsect['p_grid'], gaussian_filter(
187 potemp, sigma=1.0), range(0, 500, 5), colors='red')
188ax.clabel(cs, fmt='%i')
189
190# Plot smoothed mixing ratio grid (g/kg)
191cs = ax.contour(xsect['x_grid'], xsect['p_grid'], gaussian_filter(
192 mixrat*1000, sigma=2.0), range(0, 41, 2), colors='tab:green')
193ax.clabel(cs, fmt='%i')
194
195# Add some informative titles
196plt.title('Cross-Section from DNR to BUF Potential Temp. '
197 '(K; red) and Mix. Rat. (g/kg; green)', loc='left')
198plt.title(date, loc='right')
199
200plt.show()
|
