1# Copyright (c) 2017,2018 MetPy Developers.
2# Distributed under the terms of the BSD 3-Clause License.
3# SPDX-License-Identifier: BSD-3-Clause
4"""
5===============================
6Sigma to Pressure Interpolation
7===============================
8
9By using `metpy.calc.log_interp`, data with sigma as the vertical coordinate can be
10interpolated to isobaric coordinates.
11"""
12
13import cartopy.crs as ccrs
14import cartopy.feature as cfeature
15import matplotlib.pyplot as plt
16from netCDF4 import Dataset, num2date
17
18from metpy.cbook import get_test_data
19from metpy.interpolate import log_interpolate_1d
20from metpy.plots import add_metpy_logo, add_timestamp
21from metpy.units import units
22
23######################################
24# **Data**
25#
26# The data for this example comes from the outer domain of a WRF-ARW model forecast
27# initialized at 1200 UTC on 03 June 1980. Model data courtesy Matthew Wilson, Valparaiso
28# University Department of Geography and Meteorology.
29
30
31data = Dataset(get_test_data('wrf_example.nc', False))
32lat = data.variables['lat'][:]
33lon = data.variables['lon'][:]
34time = data.variables['time']
35vtimes = num2date(time[:], time.units)
36temperature = units.Quantity(data.variables['temperature'][:], 'degC')
37pressure = units.Quantity(data.variables['pressure'][:], 'Pa')
38hgt = units.Quantity(data.variables['height'][:], 'meter')
39
40####################################
41# Array of desired pressure levels
42plevs = [700.] * units.hPa
43
44#####################################
45# **Interpolate The Data**
46#
47# Now that the data is ready, we can interpolate to the new isobaric levels. The data is
48# interpolated from the irregular pressure values for each sigma level to the new input
49# mandatory isobaric levels. `mpcalc.log_interp` will interpolate over a specified dimension
50# with the `axis` argument. In this case, `axis=1` will correspond to interpolation on the
51# vertical axis. The interpolated data is output in a list, so we will pull out each
52# variable for plotting.
53
54height, temp = log_interpolate_1d(plevs, pressure, hgt, temperature, axis=1)
55
56####################################
57# **Plotting the Data for 700 hPa.**
58
59# Set up our projection
60crs = ccrs.LambertConformal(central_longitude=-100.0, central_latitude=45.0)
61
62# Set the forecast hour
63FH = 1
64
65# Create the figure and grid for subplots
66fig = plt.figure(figsize=(17, 12))
67add_metpy_logo(fig, 470, 320, size='large')
68
69# Plot 700 hPa
70ax = plt.subplot(111, projection=crs)
71ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.75)
72ax.add_feature(cfeature.STATES, linewidth=0.5)
73
74# Plot the heights
75cs = ax.contour(lon, lat, height[FH, 0, :, :], transform=ccrs.PlateCarree(),
76                colors='k', linewidths=1.0, linestyles='solid')
77cs.clabel(fontsize=10, inline=1, inline_spacing=7, fmt='%i', rightside_up=True,
78          use_clabeltext=True)
79
80# Contour the temperature
81cf = ax.contourf(lon, lat, temp[FH, 0, :, :], range(-20, 20, 1), cmap=plt.cm.RdBu_r,
82                 transform=ccrs.PlateCarree())
83cb = fig.colorbar(cf, orientation='horizontal', aspect=65, shrink=0.5, pad=0.05,
84                  extendrect='True')
85cb.set_label('Celsius', size='x-large')
86
87ax.set_extent([-106.5, -90.4, 34.5, 46.75], crs=ccrs.PlateCarree())
88
89# Make the axis title
90ax.set_title(f'{plevs[0]:~.0f} Heights (m) and Temperature (C)', loc='center', fontsize=10)
91
92# Set the figure title
93fig.suptitle(f'WRF-ARW Forecast VALID: {vtimes[FH]} UTC', fontsize=14)
94add_timestamp(ax, vtimes[FH], y=0.02, high_contrast=True)
95
96plt.show()
97