Meteorology in Python

# Name: narr_500mb_vor_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, 500-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.gridspec as gridspec

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 = 500

# 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', 'u-component_of_wind_isobaric', 'v-component_of_wind_isobaric')

# Use the query to obtain NetCDF data.

data = ncss.get_data(query)

# Extract data and assign units

# Geopotential height in m AMSL.

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

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

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.

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

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

p = 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)

# Calcualte 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]

# 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]

# Relative vorticity (curl) in hertz / 1/s.

vor = mpcalc.vorticity(u,  v,  dx,  dy,  dim_order = 'yx')

# Planetary vorticity (curl) in hertz / 1/s.

f = mpcalc.coriolis_parameter(np.deg2rad(lat)).to('1/s')

# Absolute vorticity (curl) in hertz / 1/s.

avor = vor + f

# Absolute vorticity advection in hertz / 1/s.

absvor_adv = mpcalc.advection(avor,  [u,  v],  (dx,  dy),  dim_order = 'yx') * 1e9

# Data projection.

dataproj = ccrs.PlateCarree()

# Plot projection.

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

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

# Cell content replaced by load magic replacement.

fig, ax = create_map_background()

# First layer are geopotential 

# height contours in m AMSL. 

geo_range = np.arange(5100, 6061, 60)

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

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

# Second layer are absolute 

# vorticity contours in 10^5

# Hz / 1/s. 

avor_range = np.arange(-9, 50, 5)

ly2 = ax.contour(lon, lat, avor*10**5, avor_range, colors = 'grey', linewidths = 1.25, linestyles = 'dashed', transform = dataproj)

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

# Third layer are colorfill 

# vorticity advections in 10^8

# Hz^2 / 1/s^2.

avoradv_range = np.arange(-30, 31, 2)

ly3 = ax.contourf(lon, lat, absvor_adv, avoradv_range[avoradv_range != 0], extend = 'both', cmap = 'bwr', transform = dataproj)

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

# Fourth layer are wind barbs 

# in kt. 

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

# Titles

plt.title('NARR 500-mb geo heights (m), avor$*10^5$ ($s^{-1}$), avor adv$*10^8$ ($s^{-2}$), and wind barbs (kt)', loc = 'left')

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

# Save figure.

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