Meteorology in Python

# Name: narr_skew_t_log_p.py

#

# Purpose: Extracts numerical weather

# prediction model output from

# North American Regional Reanalysis

# (NARR). Generates skew-T log-p

# diagram above location of interest.

#

# Caveats: Run code with Python 3.6

# or higher.

#

# Authors: Source code originally

# developed by University

# Corporation for Atmospheric Research.

# Code modified by Lawrence Burkett

# (2020) for academic use.

# Import the following modules:

import matplotlib

matplotlib.use('Agg')

from matplotlib import pyplot as plt

from datetime import datetime

import matplotlib.pyplot as plt

import metpy.calc as mpcalc

import numpy as np

from metpy.plots import SkewT

from metpy.plots import StationPlot

from metpy.units import units

from netCDF4 import Dataset, num2date

from scipy.ndimage import gaussian_filter

from siphon.catalog import TDSCatalog

# Year, month, day, and hour of interest.

year = 2004

month = 5

day = 29

hour = 12

dt = datetime(year, month, day, hour)

# Location of interest with

# respect to longitude and

# latitude in deg.

loc_lon, loc_lat = -99.84, 35.07

# Read NARR data from THREDDS

# server.

base_url = 'https://www.ncei.noaa.gov/thredds/catalog/model-narr-a-files/'

# Generate URL to year, month,

# day, and hour of interest.

cat = TDSCatalog(f'{base_url}{dt:%Y%m}/{dt:%Y%m%d}/catalog.xml')

# Have Siphon find appropriate dataset.

ds = cat.datasets.filter_time_nearest(dt)

# Interface with data through NetCDF

# Subset Service (NCSS).

ncss = ds.subset()

# Create NCSS query with desired

# specifications.

query = ncss.query()

query.lonlat_box(north = 60, south = 18, east = 300, west = 225)

query.all_times()

query.add_lonlat()

query.accept('netcdf')

query.variables('Geopotential_height_isobaric', 'Temperature_isobaric', 'u-component_of_wind_isobaric', 'v-component_of_wind_isobaric', 'Specific_humidity_isobaric')

# Use query to obtain NetCDF data.

data = ncss.get_data(query)

# Extract data and assign units.

# Geopotential height in m AMSL.

# Make this data an array.

h = np.array(gaussian_filter(data.variables['Geopotential_height_isobaric'][0], sigma = 1.0) * units.meter)

# Specific humidity in kg/kg.

# Make this data an array.

q = np.array(gaussian_filter(data.variables['Specific_humidity_isobaric'][0], sigma = 1.0) * units('kg/kg'))

# Temperature in K. Make

# this data an array.

T = np.array(gaussian_filter(data.variables['Temperature_isobaric'][0], sigma = 1.0) * units.K)

# West-to-east (zonal) wind in m/s.

# Make this data an array.

u = np.array(gaussian_filter(data.variables['u-component_of_wind_isobaric'][0], sigma = 1.0) * units('m/s'))

# South-to-north (meridional) wind in m/s.

# Make this data an array.

v = np.array(gaussian_filter(data.variables['v-component_of_wind_isobaric'][0], sigma = 1.0) * units('m/s'))

# Extract coordinate data. Make

# pressure level data an array.

lat = data.variables['lat'][:]

lon = data.variables['lon'][:]

p = np.array(data.variables['isobaric1'][:])

# Extract time. Check whether

# time2, time1, or time by printing data:

time = data.variables['time2']

vtime = num2date(time[0], units = time.units)

# Find NARR grid cell

# containing location of

# interest.

abs_lon = abs(lon - loc_lon)

abs_lat = abs(lat - loc_lat)

maxi = np.maximum(abs_lon, abs_lat)

loc_index = np.argmin(maxi)

# Get the 2D equivalent

# of index in 1D array.

ij = np.unravel_index(loc_index, (h.shape[1], h.shape[2]))

i_loc_index = ij[0]

j_loc_index = ij[1]

# Geopential heights in

# column above location.

hk = []

for k in range(len(p)):

tmp = (h[k, i_loc_index, j_loc_index])

hk.append(tmp)

# Specific humidity in

# column above location.

qk = []

for k in range(len(p)):

tmp = (q[k, i_loc_index, j_loc_index])

qk.append(tmp)

# Temperatures in column

# above location.

Tk = []

for k in range(len(p)):

tmp = (T[k, i_loc_index, j_loc_index])

Tk.append(tmp)

# Zonal wind in column

# above location.

uk = []

for k in range(len(p)):

tmp = (u[k, i_loc_index, j_loc_index])

uk.append(tmp)

# Meridional wind in

# column above location.

vk = []

for k in range(len(p)):

tmp = (v[k, i_loc_index, j_loc_index])

vk.append(tmp)

# Calculate dewpoint temperature in deg C.

Td = mpcalc.dewpoint_from_specific_humidity(qk * units('kg/kg'), Tk * units.K, p * units.hPa)

# Make lists arrays; and

# convert temperatures to

# deg C, and winds to kt.

Tk = np.array(Tk)

Tk = Tk - 273.15

uk = np.array(uk)

uk = uk * 1.94384

vk = np.array(vk)

vk = vk * 1.94384

# Create figure.

fig = plt.figure(figsize = (11, 11))

# Skew-T plot from MetPy.

skew = SkewT(fig, rotation = 45)

# Plot data.

skew.plot(p, Tk, 'r')

skew.plot(p, Td, 'g')

skew.plot_barbs(p[::3], uk[::3], vk[::3], y_clip_radius = 0.03)

# Set axis limits.

skew.ax.set_xlim(-30, 40)

skew.ax.set_ylim(1020, 100)

# Add key lines.

skew.plot_dry_adiabats(t0 = np.arange(233, 533, 10) * units.K, alpha = 0.25, color = 'orangered')

skew.plot_moist_adiabats(t0 = np.arange(233, 400, 5) * units.K, alpha = 0.25, color = 'tab:green')

skew.plot_mixing_lines(p = np.arange(1000, 99, -20) * units.hPa, linestyle = 'dotted', color = 'tab:blue')

# Add labels.

skew.ax.set_xlabel('temperature (\N{DEGREE CELSIUS})')

skew.ax.set_ylabel('pressure (mb)')

# Title

plt.title('NARR -99.84\xb0, 35.07\xb0: temp (red), dwpt (green), and wind barbs (kt)', loc = 'left')

plt.title(f'valid: {vtime} utc', loc = 'right')

# Save figure.

plt.savefig('narr-5-29-2004-12utc-sounding.png')