THREDDS-served GRIB Data: GFS (Global Grid) Model: Eclipse Cloud Cover Forecast¶

This notebook reads in realtime GFS model output from a THREDDS server.¶

Overview:¶

Determine latest available model run
Use dictionaries to deal with coordinate dimension names that change from run-to-run
Subset data and map regions
Create a 4-panel plot

Imports¶

import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader
from matplotlib import pyplot as plt
from matplotlib.dates import DateFormatter, AutoDateLocator,HourLocator,DayLocator
from matplotlib.collections import LineCollection
from matplotlib.colors import BoundaryNorm, ListedColormap
import numpy as np
import xarray as xr
import pandas as pd
from datetime import datetime, timedelta
from pyproj import Proj
from metpy.plots import USCOUNTIES
import seaborn as sns
import geopandas as gpd
sns.set()

Access real-time NWP data from Unidata’s THREDDS server¶

Parse the current time and choose a recent hour when we might expect the model run to be complete. There is normally a ~90 minute lag, but we’ll give it a bit more than that.

now = datetime.now()
year = now.year
month = now.month
day = now.day
hour = now.hour
minute = now.minute

# Allow for a six-hour lag between GFS initialization and current time
if (hour >= 18):
    runHour = 12
    hrDelta = hour - runHour
elif (hour >= 12):
    runHour = 6
    hrDelta = hour - runHour
elif (hour >= 6):
    runHour = 0
    hrDelta = hour - runHour
else:
    runHour = 18
    hrDelta = hour + 6
    
runTime = now - timedelta(hours=hrDelta)
runTimeStr = runTime.strftime('%Y%m%d %H00 UTC')
modelDate = runTime.strftime('%Y%m%d')
modelHour = runTime.strftime('%H')
modelMin = runTime.strftime('%M')
modelDay = runTime.strftime('%D')
titleTime = modelDay + ' ' +  modelHour + '00 UTC'
print (modelDay)
print (modelHour)
print (runTimeStr)
print(titleTime)

04/29/24
18
20240429 1800 UTC
04/29/24 1800 UTC

Open the most recent NCEP real-timeGFS from Unidata’s THREDDS server

ds = xr.open_dataset(f'https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/Global_0p25deg/GFS_Global_0p25deg_{modelDate}_{modelHour}00.grib2')

Leverage MetPy’s parse_cf function, which among other things will make it easy to extract the dataset’s map projection.

ds = ds.metpy.parse_cf()

Note that there now exists a metpy_crs coordinate variable, which we’ll soon use to determine and use the Dataset’s map projection.

Select some of the variables from the Dataset; in this case, we’ll choose those relating to cloud cover.

varTotal = ds.Total_cloud_cover_entire_atmosphere
varLow = ds.Low_cloud_cover_low_cloud
varMed = ds.Medium_cloud_cover_middle_cloud
varHigh = ds.High_cloud_cover_high_cloud

Examine one of the grids

GRIB files that are served by THREDDS are automatically converted into NetCDF. A quirk in this translation can result in certain dimensions (typically those that are time or vertically-based) varying from run-to-run. For example, the time coordinate dimension for a given variable might be time in one run, but time1 (or time2, time3, etc.) in another run. This poses a challenge any time we want to subset a DataArray. For example, if a particular dataset had vertical levels in units of hPa, the dimension name for the geopotential height grid might be isobaric. If we wanted to select only the 500 and 700 hPa levels, we might write the following line of code:

Z = ds.sel(isobaric=[500, 700])

This would fail if we pointed to model output from a different time, and the vertical coordinate name for geopotential height was instead isobaric1.

Via strategic use of Python dictionaries, we can programmatically deal with this!

First, let’s use some MetPy accessors that determine the dimension names for times and vertical levels:

timeDim= varTotal.metpy.time.name
timeDim

'time1'

We next construct dictionaries corresponding to each desired subset of the relevant dimensions. Here, we specify a desired date and time:

verifTime = [datetime(2024, 4, 8, 18), datetime(2024, 4, 8, 21)]

timeDict = {timeDim: verifTime}

timeDict

{'time1': [datetime.datetime(2024, 4, 8, 18, 0),
  datetime.datetime(2024, 4, 8, 21, 0)]}

Get the map projection from one of the DataArrays. They will be the same across the Dataset.

proj_data = varTotal.metpy.cartopy_crs
proj_data

<cartopy.crs.PlateCarree object at 0x1468e540ac90>

Create DataArray objects for the horizontal coordinate dimensions. In this case, the model is output on global cylindrical equidistant (PlateCarree) projection, which means the horizontal dimensions are named lon and lat and represent grid point locations in degrees.

x, y = ds.lon, ds.lat

ATM350 Spring 2024

THREDDS-served GRIB Data: GFS (Global Grid) Model: Eclipse Cloud Cover Forecast

Contents

THREDDS-served GRIB Data: GFS (Global Grid) Model: Eclipse Cloud Cover Forecast¶

This notebook reads in realtime GFS model output from a THREDDS server.¶

Overview:¶

Imports¶

Access real-time NWP data from Unidata’s THREDDS server¶

Subset the dataset to a smaller region¶

We can directly select latitude and longitude values to perform regional subset selections; no need to use the Pyproj library.¶