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