Meteorology in Python

# Name: narr_850mb_temp_adv_map.py

#

# Purpose: Extracts numerical weather 

# prediction model output from 

# North American Regional Reanalysis

# (NARR). Calculates several key 

# atmospheric properties from 

# output. Plots contours, 

# filled-contours, and barbs from output.

# Here I focus on creating a synoptic scale, 850-mb map.  

#

# Caveats: Run code with Python 3.6 

# or higher. Code requires making 

# modification to one line of MetPy's 

# cartopy_utils.py code to function 

# with later releases of Cartopy 

# (see footnote [1] below).  

# 

# 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 cartopy.crs as ccrs

import cartopy.feature as cfeature

import matplotlib.pyplot as plt

import metpy.calc as mpcalc

import numpy as np

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)

# Pressure level of interest in mb.  

p0 = 850

# 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')

# Use query to obtain NetCDF data. 

data = ncss.get_data(query)

# Extract data and assign units. 

# Temperature in K. 

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

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

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

# Geopotential height in m AMSL.

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

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

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

# Extract coordinate data for plotting (i.e., latitude, longitude, and height).

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

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

p = data.variables['isobaric1'][:]

# Extract time for plotting. Check whether

# time2, time1, or time by printing data: 

# print(data)

time = data.variables['time2']

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

# Calculate dx and dy for calculations.

dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)

# Model pressure level k in (i, j, k) 

# grid index coordinate system. 

k = np.where(p == p0)[0][0]

# Geopotential height at k-th 

# index in m.

h = h[k]

# Temperature at k-th index

# in K. 

T = T[k]

# West-to-east (zonal) wind 

# at k-th index in m/s.  

u = u[k]

# South-to-north (meridional) 

# wind at k-th index in m/s.  

v = v[k]

# Temperature advection in C/s.

T_adv = mpcalc.advection(T,  [u,  v], (dx,  dy), dim_order = 'yx').to('degC/s')

# Data projection. 

dataproj = ccrs.PlateCarree()

# Plot projection.

plotproj = ccrs.LambertConformal(central_longitude = -100.,  central_latitude = 40.,  standard_parallels = [30,  60])

# Define background map. 

def create_map_background():

    fig = plt.figure(figsize = (14, 12))

    ax = plt.subplot(111, projection = plotproj)

    ax.set_extent([-125,  -73,  25,  50],ccrs.PlateCarree())

    ax.coastlines("50m",  linewidth = 0.75)

    ax.add_feature(cfeature.STATES, linewidth = 0.5)

    return fig, ax

# Create figure. 

fig, ax = create_map_background()

# First layer are temperatures 

# in deg. C contours. 

ly1 = ax.contour(lon, lat, T.to('degC'), range(-50, 50, 2), colors = 'grey', linestyles = 'dotted', transform = dataproj)

plt.clabel(ly1, fontsize = 10, inline = 1, inline_spacing = 10, fmt = '%i', rightside_up = True, use_clabeltext = True)

# Second layer are geopotential 

# heights in m AMSL contours.

geo_range = np.arange(0, 4000, 30)

ly2 = ax.contour(lon, lat, h, geo_range, colors = 'k', linewidths = 1.0, linestyles = 'solid', transform = dataproj)

plt.clabel(ly2, fontsize = 10, inline = 1, inline_spacing=10, fmt = '%i', rightside_up = True, use_clabeltext = True)

# Third layer are colorfill temperature 

# advection in deg. C/h.

adv_range = [-3, -2.2, -2, -1.5, -1, -0.5, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0]

cf = ax.contourf(lon, lat, T_adv*3600, adv_range, cmap = 'bwr', extend = 'both', transform = dataproj) 

plt.colorbar(cf, orientation = 'horizontal', pad = 0, aspect = 50, extendrect = True, ticks = adv_range)

# Fourth layer are wind barbs in kt. 

ax.barbs(lon, lat, u.to('kts').m, v.to('kts').m, regrid_shape = 15, transform = dataproj)

# Titles

plt.title('NARR 850-mb geo heights (m), temp (\xb0C), temp adv (\xb0C hr$^{-1}$), and wind barbs (kt)', loc = 'left')

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

# Save figure.

plt.savefig('narr-5-29-2004-12utc-850mb-temp-adv-map.png')