Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

ERA5 SoS Template

This notebook will create all the necessary files and directories for a visualization on NOAA’s Science-on-a-Sphere.

Overview:

  1. Set up output directories for SoS

  2. Specify the date/time range for your case

  3. Use Cartopy’s add_cyclic_point method to avoid a blank seam at the dateline

  4. Create small and large thumbnail graphics.

  5. Create a standalone colorbar

  6. Create a set of labels for each plot

  7. Create several days’ worth of Science-on-a-Sphere-ready plots

  8. Create an SoS playlist file

Imports

import xarray as xr
import numpy as np
import pandas as pd
import cartopy.feature as cfeature
import cartopy.crs as ccrs
import cartopy.util as cutil
from datetime import datetime as dt
import warnings
import metpy.calc as mpcalc
from metpy.units import units
import matplotlib.pyplot as plt

# Do not output warning messages

warnings.simplefilter("ignore")

Set up output directories for SoS

The software that serves output for the Science-on-a-Sphere expects a directory structure as follows:

  • Top level directory: choose a name that is consistent with your case, e.g.: SS93

  • 2nd level directory: choose a name that is consistent with the graphic, e.g.: SLP_500Z

  • 3rd level directories:

    • 2048: Contains the graphics (resolution: 2048x1024) that this notebook generates

    • labels: Contains one or two files:

      1. (required) A text file, labels.txt, which has as many lines as there are graphics in the 2048 file. Each line functions as the title of each graphic.

      2. (optional) A PNG file, colorbar.png, which colorbar which will get overlaid on your map graphic.

    • media: Contains large and small thumbnails (thumbnail_small.jpg, thumbnail_large.jpg) that serve as icons on the SoS iPad and SoS computer apps

    • playlist: A text file, playlist.sos, which tells the SoS how to display your product

As an example, here is how the directory structure on our SoS looks for the products generated by this notebook. Our SoS computer stores locally-produced content in the /shared/sos/media/site-custom directory (note: the SoS directories are not network-accessible, so you won’t be able to cd into them). The ERA5 visualizations go in the ERA5 subfolder. Your top-level directory sits within the ERA5 folder.

sos@sos1:/shared/sos/media/site-custom/ERA5/cle_superbomb/SLP_500Z$ ls -R
.:
2048  labels  media  playlist

./2048:
ERA5_1978012600-fs8.png  ERA5_1978012606-fs8.png  ERA5_1978012612-fs8.png  ERA5_1978012618-fs8.png
ERA5_1978012601-fs8.png  ERA5_1978012607-fs8.png  ERA5_1978012613-fs8.png  ERA5_1978012619-fs8.png
ERA5_1978012602-fs8.png  ERA5_1978012608-fs8.png  ERA5_1978012614-fs8.png  ERA5_1978012620-fs8.png
ERA5_1978012603-fs8.png  ERA5_1978012609-fs8.png  ERA5_1978012615-fs8.png  ERA5_1978012621-fs8.png
ERA5_1978012604-fs8.png  ERA5_1978012610-fs8.png  ERA5_1978012616-fs8.png  ERA5_1978012622-fs8.png
ERA5_1978012605-fs8.png  ERA5_1978012611-fs8.png  ERA5_1978012617-fs8.png  ERA5_1978012623-fs8.png

./labels:
colorbar.png  labels.txt

./media:
thumbnail_big.jpg  thumbnail_small.jpg

./playlist:
playlist.sos

Define the 1st and 2nd-level directories. The 3rd-level directories follow on from them. Then create the directories. Remove previously-created PNGs.

# You define these
caseDir = 'changethis'
prodDir = 'changethis'

# These remain as is
graphicsDir = caseDir + '/' + prodDir + '/2048/'
labelsDir = caseDir + '/' + prodDir + '/labels/'
mediaDir = caseDir + '/' + prodDir + '/media/'
playlistDir = caseDir + '/' + prodDir + '/playlist/'

! mkdir -p {graphicsDir} {labelsDir} {mediaDir} {playlistDir}

# Remove all PNGs in the graphics directory - COMMENT OUT if you don't want to do this!
! rm -f  {graphicsDir}*.png

Specify a starting and ending date/time, regional extent, vertical levels, and access/subset the ERA5

# Set the bounds of the map ... for the SoS, we will cover the entire globe.
## Do not change these values ##
latN = 90
latS = -90
lonW = 0
lonE = 360

latRange = np.arange(latS, latN + 0.1,.5) # ensure +90 is included
lonRange = np.arange(lonW, lonE, .5) 
## 

# Select the map projection for your figure.
proj_map = ccrs.PlateCarree() # For SoS maps

# Select the projection on which the dataset is based
proj_data = ccrs.PlateCarree()

# Select the resolution of the Cartopy cartographic features: coarsest Natural Earth, and no counties

res = '110m' # Least detailed, best for large/global maps

# Specify desired pressure levels: in this case, a list of one or more
plevel_list = [500, 700, 850]

startYear = 2026
startMonth = 2
startDay = 20
startHour = 0
startMinute = 0
startDateTime = dt(startYear,startMonth,startDay, startHour, startMinute)

endYear = 2026
endMonth = 2
endDay = 24
endHour = 18
endMinute = 0
endDateTime = dt(endYear,endMonth,endDay, endHour, endMinute)

delta_time = endDateTime - startDateTime
time_range_max = 7*86400

if (delta_time.total_seconds() > time_range_max):
    raise RuntimeError("Your time range must not exceed 7 days. Go back and try again.")

# Create a list of date and times based on what we specified for the initial and final times, 
#  using Pandas’ date_range function
dateRange = pd.date_range(startDateTime, endDateTime,freq="6h")

Access either the cloud-served or local ERA5 repository. Then subset it according to your choices above.

endDate = dt(2023,1,10)
if (endDateTime <= endDate): # Use WeatherBench archive

    cloud_source = True
    ds = xr.open_dataset(
     'gs://weatherbench2/datasets/era5/1959-2023_01_10-wb13-6h-1440x721.zarr', 
     chunks={'time': 48},
     consolidated=True,
     engine='zarr'
    )
    
    # Attach units to the pressure coordinate
    ds.coords['level'].attrs['units'] = 'hPa'
    
    # Rename the variable names in this dataset so they use their corresponding short_name attributes
    # Construct a dictionary whose keys are the original variable names, and whose values are their 
    # short names
    
    rename_data_vars = {}
    seen_short_names = set()
    
    for var_name in ds.data_vars:
        short = ds[var_name].attrs.get("short_name")
        # Only rename if:
        # 1. short_name exists
        # 2. We haven't already used it
        if short and short not in seen_short_names:
            rename_data_vars[var_name] = short
            seen_short_names.add(short)
    # Apply renaming
    ds = ds.rename(rename_data_vars)

else: # Use local archive
    import glob, os
    cloud_source = False
    input_directory = '/free/ktyle/era5'
    files = glob.glob(os.path.join(input_directory,'*_era5.nc'))
    ds = xr.open_mfdataset(files,coords='minimal',compat='override')
    # Rename two of the coordinate variables so they match what is in the WeatherBench archive
    ds = ds.rename({'valid_time': 'time', 'pressure_level': 'level'})

# Perform the initial temporal and vertical level subset
ds = ds.sel(time = dateRange, longitude = lonRange, latitude = latRange, level = plevel_list)

Examine the subsetted Dataset

ds

To avoid contouring problems across the dateline, create a function so longitudes run from -180 to 180 instead of 0 to 360.

def reorder_longitudes(da):
    """
    This function reorders the longitudes in a global grid from 0 --> 360 to -180 --> 180
    """
    da = da.assign_coords(longitude=(((da.longitude + 180 ) % 360 ) - 180))
    return da.sortby(da['longitude'])

Create DataArray objects for the data variables we are interested in. Reorder their longitude dimension so they begin at 180 instead of 0.

# Select only those variables that you need for your product. 
# Examples below include Z and SLP
Z = reorder_longitudes(ds['z'])
SLP = reorder_longitudes(ds['msl'])

Select just a single isobaric surface from our initial list. Subset, and then load in the data into RAM via the compute function.

p_level = 850 
%time Zp = Z.sel(level=p_level).compute()
# 3-d variables do not have the isobaric level dimension
%time SLPp = SLP.compute()

Create objects for the relevant coordinate arrays; in this case, longitude, latitude, and time.

Remember to use one of the variables that you selected above ... in this example, we’ll use temperature on the selected isobaric surface (Tp)

lons, lats, times = Tp.longitude, Tp.latitude, Tp.time
lons

Perform any necessary diagnostics, rescaling, smoothing and/or unit conversions.

As written, this will work for the following set of overlays:

Adjust to suit your desired map!

  1. Filled contours of 850 hPa temperature

  2. Contour lines of 850 hPa geopotential height

  3. 850 hPa wind barbs

Zdam = mpcalc.geopotential_to_height(Zp).metpy.convert_units('dam')
Ukts = Up.metpy.convert_units('kts')
Vkts = Vp.metpy.convert_units('kts')
Tc = Tp.metpy.convert_units('degC')

Set contour levels; inspect min and max values for filled contours.

Tc.min().values,Tc.max().values
tContours = np.arange(-45,39,3)

# If plotting geopotential height contours, set intervals as follows
if (p_level >= 700):
    incZ = 3
elif (p_level >= 400):
    incZ = 6
elif (p_level >= 150):
    incZ = 9
else: 
    incZ = 12

zContours = np.arange(0,3000,incZ) # use for contour line
#zContours = np.arange(xxx, yyy, incZ) # use for contour fill

Set colors for map backgrounds and wind barbs

coast_border_col = 'dimgrey'
state_col = 'lightgrey'
wind_col = 'whitesmoke'

Create a single plot for the first time in the time range (isel=0), just to ensure all looks good. Add the cyclic point to avoid the seam at the eastern edge of the domain.

# Edit to correspond to your variables

T0 = Tc.isel(time=0)
U0 = Ukts.isel(time=0)
V0 = Vkts.isel(time=0)
Z0 = Zdam.isel(time=0)

# add cyclic points to data and longitudes
# latitudes are unchanged in 1-dimension
Z_cyc, lons_cyc = cutil.add_cyclic_point(Z0, lons)
T_cyc, lons_cyc = cutil.add_cyclic_point(T0, lons)
U_cyc, lons_cyc = cutil.add_cyclic_point(U0, lons)
V_cyc, lons_cyc = cutil.add_cyclic_point(V0, lons)

dpi = 100
fig = plt.figure(figsize=(2048/dpi, 1024/dpi))
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree(), frameon=False)

ax.set_global()
ax.add_feature(cfeature.COASTLINE.with_scale(res), edgecolor=coast_border_col)
ax.add_feature(cfeature.BORDERS.with_scale(res), edgecolor=coast_border_col)
ax.add_feature(cfeature.STATES.with_scale(res), edgecolor=state_col, linewidth=0.75)

# Temperature (T) contour fills

# Note we don't need the transform argument since the map/data projections are the same, but we'll leave it in
CF = ax.contourf(lons_cyc,lats,T_cyc,levels=tContours,cmap=plt.get_cmap('coolwarm'), extend='both', transform=proj_data) 

# Height (Z) contour lines
CL = ax.contour(lons_cyc,lats,Z_cyc,zContours,linewidths=1.25,colors='yellow', transform=proj_data)
ax.clabel(CL, inline_spacing=0.2, fontsize=8, fmt='%.0f')
fig.tight_layout(pad=.01)

# Plotting wind barbs uses the ax.barbs method. Here, you can't pass in the DataArray directly; you can only pass in the array's values.
# Also need to sample (skip) a selected # of points to keep the plot readable.
skip = 8
# Default length of wind barbs is 7. Shorten them a bit.
ax.barbs(lons_cyc[::skip],lats[::skip],U_cyc[::skip,::skip], V_cyc[::skip,::skip], color=wind_col, length = 5, transform=proj_data)
# Save this figure to your current directory

### Choose your own name for this test figure
fig.savefig('test_ERA5_SOS.png')

If you’d like, go to https://tcpoole.com/WorldForge3dFree/ , upload your PNG, and view how your graphic might look on the sphere.

Create thumbnails, colorbar and text labels.

res = '110m'
dpi = 100
for size in (128, 800):
    if (size == 128):
        sizeStr = 'small'
    else:
        sizeStr ='big'
        
    fig = plt.figure(figsize=(size/dpi, size/dpi))
    ax = plt.subplot(1,1,1,projection=ccrs.Orthographic(central_longitude=-90), frameon=False)
    ax.set_global()
    ax.add_feature(cfeature.COASTLINE.with_scale(res))
    tl1 = caseDir
    tl2 = prodDir
    ax.set_title(f"{tl1}\n{tl2}", color='purple', fontsize=8)

    # Contour fills

    CF = ax.contourf(lons_cyc,lats,T_cyc, levels=tContours,cmap=plt.get_cmap('coolwarm'), extend='both', transform=proj_data) 

    fig.tight_layout(pad=.01)
    fig.savefig(f'{mediaDir}thumbnail_{sizeStr}.jpg')

# draw a new figure and replot the colorbar there
fig,ax = plt.subplots(figsize=(14,2), dpi=125)
# set tick and ticklabel color
tick_color='black'
label_color='orange'
cbar = fig.colorbar(CF, ax=ax, orientation='horizontal')
cbar.ax.xaxis.set_tick_params(color=tick_color, labelcolor=label_color, labelsize=8)
# Remove the Axes object ... essentially the contours and cartographic features ... from the figure.
ax.remove()
# All that remains is the colorbar ... save it to disk. Make the background transparent.
fig.savefig(f'{labelsDir}colorbar.png',transparent=True)

# create labels
labelFname = f'{labelsDir}labels.txt'

Define a function to create the plot for each hour.

Edit the contour and wind barb lines to correspond to the variables you are overlaying!

def make_sos_map (timeIdx):
    
    res = '110m'
    dpi = 100
    fig = plt.figure(figsize=(2048/dpi, 1024/dpi))
    ax = plt.subplot(1,1,1,projection=ccrs.PlateCarree(), frameon=False)
    ax.set_global()
    ax.add_feature(cfeature.COASTLINE.with_scale(res), edgecolor=coast_border_col)
    ax.add_feature(cfeature.BORDERS.with_scale(res), edgecolor=coast_border_col)
    ax.add_feature(cfeature.STATES.with_scale(res), edgecolor=state_col, linewidth=0.75)

#   add cyclic points to data and longitudes
#   latitudes are unchanged in 1-dimension
    Z_cyc, lons_cyc= cutil.add_cyclic_point(Zdam.isel(time=timeIdx), lons)
    T_cyc, lons_cyc = cutil.add_cyclic_point(Tc.isel(time=timeIdx), lons)
    U_cyc, lons_cyc = cutil.add_cyclic_point(Ukts.isel(time=timeIdx), lons)
    V_cyc, lons_cyc = cutil.add_cyclic_point(Vkts.isel(time=timeIdx), lons)

#   Temperature (T) contour fills

#   Note we don't need the transform argument since the map/data projections are the same, but we'll leave it in
    CF = ax.contourf(lons_cyc,lats,T_cyc,levels=tContours,cmap=plt.get_cmap('coolwarm'), extend='both', transform=proj_data) 

#   Height (Z) contour lines
    CL = ax.contour(lons_cyc,lats,Z_cyc,zContours,linewidths=1.25,colors='yellow', transform=proj_data)
    ax.clabel(CL, inline_spacing=0.2, fontsize=8, fmt='%.0f')
    fig.tight_layout(pad=.01)

#   Plotting wind barbs uses the ax.barbs method. Here, you can't pass in the DataArray directly; you can only pass in the array's values.
#   Also need to sample (skip) a selected # of points to keep the plot readable.
    skip = 8
#   Default length of wind barbs is 7. Shorten them a bit.
    ax.barbs(lons_cyc[::skip],lats[::skip],U_cyc[::skip,::skip], V_cyc[::skip,::skip], color=wind_col, length = 5, transform=proj_data)
    frameNum = f'{timeIdx}'.zfill(2)
    figName = f'{graphicsDir}ERA5_{caseDir}_{prodDir}_{frameNum}.png'
    fig.savefig(figName)
    
    # Reduce the size of the PNG image via the Linux pngquant utility. The -f option overwites the resulting file if it already exists.
    # The output file will end in "-fs8.png"
    ! pngquant -f {figName}
    
    # Remove the original PNG
    ! rm -f {figName}
    
    # Do not show the graphic in the notebook
    plt.close()

Loop through the times to produce the individual frames.

nTimes = len(dateRange)

labelFileObject = open(labelFname, 'w')

nFrames = 2
## Uncomment this next line once you are ready to generate all the plots for your time range
#nFrames = nTimes 

for timeIdx in range(0, nFrames, 1):
    
    make_sos_map(timeIdx)
    
    # Construct the title string and write it to the file
    valTime = pd.to_datetime(times[timeIdx].values).strftime("%m/%d/%Y %H UTC")
    tl1 = f'ERA5 {p_level} hPa Z/T/Wind {valTime} \n' # \n is the newline character
    labelFileObject.write(tl1)
    print(tl1)
    
# Close the text file
labelFileObject.close()
print("Finished with SoS graphic loop")

Create the SoS playlist file

# Set the author (creator) name
creaName = "Rovin Lazyle" # Put your name here!

playlistFname = f'{playlistDir}playlist.sos'
plFileObject = open(playlistFname, 'w')
subCat = "ERA5"
cRet = "\n" # New line character code
lBr = "{"
rBr = "}"

nameStr = f'ERA5 {caseDir}{prodDir}'
descriptionStr = f'{lBr}{caseDir}{prodDir}{rBr}'

plFileObject.write(f"name = {nameStr}{cRet}")
plFileObject.write(f"description = {descriptionStr}{cRet}")

plFileObject.write(f"pip = ../labels/colorbar.png{cRet}")
plFileObject.write(f"pipheight = 10{cRet}")
plFileObject.write(f"pipvertical = -35{cRet}")

plFileObject.write(f"label = ../labels/labels.txt{cRet}")
plFileObject.write(f"layer = Grids{cRet}")
plFileObject.write(f"layerdata = ../2048{cRet}")

plFileObject.write(f"firstdwell = 2000{cRet}")
plFileObject.write(f"lastdwell = 3000{cRet}")
plFileObject.write(f"fps = 8{cRet}")

plFileObject.write(f"zrotationenable = 1{cRet}")
plFileObject.write(f"zfps = 30{cRet}")

plFileObject.write(f"subcategory = {subCat}{cRet}")
plFileObject.write(f"Creator = {creaName}{cRet}")
                   
plFileObject.close()

Summary

This notebook presents a streamlined version of the steps necessary to produce a SoS visualization.

Things to try

  1. Copy this template notebook and rename it so it fits the particulars of your SoS product. E.g., 700_Wind_Temperature_Frontogenesis.ipynb

  2. Modify the copied notebook to fit your desired visualization.

  3. Don’t forget to perform unit conversions, scaling, and smoothing, as necessary for your product.

Week 13 (Apr. 21, 23)
ERA5 Data and Graphics Preparation for the Science-on-a-Sphere
HW3 Q2 (METAR/NYSM Temperatures) Solution
HW3_Q2_Metpy_METAR_NYSM: Example Solution