Meteorology in Python

# Name: narr_grid_cell.py

#

# Purpose: Extracts numerical weather

# prediction model output from

# North American Regional Reanalysis

# (NARR) from a specific location at

# constant pressure. Here I focus

# on output from 850-mb.

#

# 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

# (2019) for academic use.

# Import the following modules:

from datetime import datetime

import numpy as np

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

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

# Temperature in K.

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

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

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

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

# grid index coordinate system.

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

# Geopotential height in m

# at location of interest

# at pressure k.

hk = h[k]

hk_flat = hk.flatten()

loc_h = hk_flat[loc_index]

# Temperature in K at

# location of interest

# at pressure k.

Tk = T[k]

Tk_flat = Tk.flatten()

loc_T = Tk_flat[loc_index]

# Zonal wind in m/s at

# location of interest.

# at pressure k.

uk = u[k]

uk_flat = uk.flatten()

loc_u = uk_flat[loc_index]

# Meridional wind in m/s

# at location of interest.

# at pressure k.

vk = v[k]

vk_flat = vk.flatten()

loc_v = vk_flat[loc_index]

# Print output at location

# of interest at pressure k.

print("###Requested output###")

print("geopotential height:")

print(loc_h)

print("temperature:")

print(loc_T)

print("zonal wind:")

print(loc_u)

print("meridional wind:")

print(loc_v)

###Requested output###

geopotential height:

1437.5032958984375 meter

temperature:

294.10748291015625 kelvin

zonal wind:

6.219133377075195 meter / second

meridional wind:

19.543563842773438 meter / second