{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## ERA5 Data and Graphics Preparation for the Science-on-a-Sphere" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook will create all the necessary files and directories for a visualization on NOAA's Science-on-a-Sphere. For this example, we will produce a map of 850 hPa geopotential heights, temperature, and horizontal wind and sea-level pressure, using the ERA5, for the March 13-14 1993 winter storm, aka **Superstorm '93**." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Overview:\n", "\n", "1. Set up output directories for SoS\n", "1. Specify the date/time range for your case\n", "1. Use Cartopy's `add_cyclic_point` method to avoid a blank seam at the dateline\n", "1. Create small and large thumbnail graphics.\n", "1. Create a standalone colorbar\n", "1. Create a set of labels for each plot\n", "1. Create 5 days worth of Science-on-a-Sphere-ready plots\n", "1. Create an SoS playlist file\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Prerequisites\n", "\n", "| Concepts | Importance | Notes |\n", "| --- | --- | --- |\n", "| Matplotlib | Necessary | |\n", "| Cartopy | Necessary | |\n", "| Xarray | Necessary | |\n", "| Metpy | Necessary | |\n", "| Linux command line / directory structure | Helpful | |\n", "\n", "* **Time to learn**: 30 minutes\n", "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Imports" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import xarray as xr\n", "import numpy as np\n", "import pandas as pd\n", "import cartopy.feature as cfeature\n", "import cartopy.crs as ccrs\n", "import cartopy.util as cutil\n", "from datetime import datetime as dt\n", "import warnings\n", "import metpy.calc as mpcalc\n", "from metpy.units import units\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Do not output warning messages" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "warnings.simplefilter(\"ignore\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Set up output directories for SoS" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The software that serves output for the Science-on-a-Sphere expects a directory structure as follows:\n", "- Top level directory: choose a name that is consistent with your case, e.g.: **SS93**\n", "- 2nd level directory: choose a name that is consistent with the graphic, e.g.: **SLP_500Z**\n", "- 3rd level directories: \n", " - **2048**: Contains the graphics (resolution: 2048x1024) that this notebook generates\n", " - **labels**: Contains one or two files:\n", " 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. \n", " 1. (optional) A PNG file, `colorbar.png`, which colorbar which will get overlaid on your map graphic.\n", " - **media**: Contains large and small thumbnails (`thumbnail_small.jpg`, `thumbnail_large.jpg`) that serve as icons on the SoS iPad and SoS computer apps\n", " - **playlist**: A text file, `playlist.sos`, which tells the SoS how to display your product\n", " \n", "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.\n", "\n", "```\n", "sos@sos1:/shared/sos/media/site-custom/ERA5/cle_superbomb/SLP_500Z$ ls -R\n", ".:\n", "2048 labels media playlist\n", "\n", "./2048:\n", "ERA5_1978012600-fs8.png ERA5_1978012606-fs8.png ERA5_1978012612-fs8.png ERA5_1978012618-fs8.png\n", "ERA5_1978012601-fs8.png ERA5_1978012607-fs8.png ERA5_1978012613-fs8.png ERA5_1978012619-fs8.png\n", "ERA5_1978012602-fs8.png ERA5_1978012608-fs8.png ERA5_1978012614-fs8.png ERA5_1978012620-fs8.png\n", "ERA5_1978012603-fs8.png ERA5_1978012609-fs8.png ERA5_1978012615-fs8.png ERA5_1978012621-fs8.png\n", "ERA5_1978012604-fs8.png ERA5_1978012610-fs8.png ERA5_1978012616-fs8.png ERA5_1978012622-fs8.png\n", "ERA5_1978012605-fs8.png ERA5_1978012611-fs8.png ERA5_1978012617-fs8.png ERA5_1978012623-fs8.png\n", "\n", "./labels:\n", "colorbar.png labels.txt\n", "\n", "./media:\n", "thumbnail_big.jpg thumbnail_small.jpg\n", "\n", "./playlist:\n", "playlist.sos\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the 1st and 2nd-level directories." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# You define these\n", "caseDir = 'SS93'\n", "prodDir = '850_T_Z_Wind'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The 3rd-level directories follow from the 1st and 2nd." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# These remain as is\n", "graphicsDir = caseDir + '/' + prodDir + '/2048/'\n", "labelsDir = caseDir + '/' + prodDir + '/labels/'\n", "mediaDir = caseDir + '/' + prodDir + '/media/'\n", "playlistDir = caseDir + '/' + prodDir + '/playlist/'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create these directories via a Linux command" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "! mkdir -p {graphicsDir} {labelsDir} {mediaDir} {playlistDir}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Note: In a Jupyter notebook, the ! magic indicates that what follows is a Linux command.
\n", "The -p option for mkdir will create all subdirectories, and also will do nothing if the directories already exist.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Remove all PNGs in the graphics directory - COMMENT OUT if you don't want to do this!\n", "! rm -f {graphicsDir}*.png" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1) Specify a starting and ending date/time, and access ERA5 data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Note: The sphere animation will be more effective with at least a minimum number of frames. Since the ERA5 repositories have a timestep of 6 hours (4x/day), let's set our time range to 5 days ... i.e. 20 frames.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "# Date/Time specification\n", "startYear = 1993\n", "startMonth = 3\n", "startDay = 11\n", "startHour = 0\n", "startMinute = 0\n", "startDateTime = dt(startYear,startMonth,startDay, startHour, startMinute)\n", "startDateTimeStr = dt.strftime(startDateTime, format=\"%Y%m%d\")\n", "\n", "endYear = 1993\n", "endMonth = 3\n", "endDay = 15\n", "endHour = 18\n", "endMinute = 0\n", "endDateTime = dt(endYear,endMonth,endDay, endHour, endMinute)\n", "\n", "# Vertical level specification\n", "pLevel = 850\n", "levStr = f'{pLevel}'\n", "\n", "delta_time = endDateTime - startDateTime\n", "time_range_max = 5*86400\n", "\n", "if (delta_time.total_seconds() > time_range_max):\n", " raise RuntimeError(\"Your time range must not exceed 5 days. Go back and try again.\")\n", "\n", "if (delta_time.total_seconds() < 0):\n", " raise RuntimeError(\"Your end time must not be earlier than your start time! Go back and try again.\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Wb2EndDate = dt(2023,1,10)\n", "if (endDateTime <= Wb2EndDate):\n", " ds = xr.open_dataset(\n", " 'gs://weatherbench2/datasets/era5/1959-2023_01_10-wb13-6h-1440x721.zarr', \n", " chunks={'time': 48},\n", " consolidated=True,\n", " engine='zarr'\n", ")\n", "else: \n", " import glob, os\n", " input_directory = '/free/ktyle/era5'\n", " files = glob.glob(os.path.join(input_directory,'*.nc'))\n", " varDict = {'valid_time': 'time', \n", " 'pressure_level': 'level',\n", " 'msl': 'mean_sea_level_pressure',\n", " 'q': 'specific_humidity',\n", " 't': 'temperature',\n", " 'u': 'u_component_of_wind',\n", " 'v': 'v_component_of_wind',\n", " 'w': 'vertical_velocity',\n", " 'z': 'geopotential'}\n", " dimDict = {'valid_time': 'time',\n", " 'pressure_level': 'level'}\n", " ds = xr.open_mfdataset(files).rename_dims(dimDict).rename_vars(varDict)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Examine the `Dataset`" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ds" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 2) Specify a date/time range, and subset the desired `Dataset`s along their dimensions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create a list of date and times based on what we specified for the initial and final times, using Pandas' `date_range` function" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dateList = pd.date_range(startDateTime, endDateTime,freq=\"6h\")\n", "dateList" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Now create objects for our desired DataArrays based on the coordinates we have subsetted." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Lat-lon subsetting: We will be plotting over the entire globe, so we will not restrict the range of latitude/longitude coordinates. However, instead of loading the full 1/4 degree resolution, we will sample every other point (i.e., 1/2 degree).
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lonRange = np.arange(0,360,0.5)\n", "latRange = np.arange(-90,90.1,0.5) # Why 90.1?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
To avoid contouring problems across the dateline, create a function so longitudes run from -180 to 180 instead of 0 to 360.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def reorder_longitudes(da):\n", " \"\"\"\n", " This function reorders the longitudes in a global grid from 0 --> 360 to -180 --> 180\n", " \"\"\"\n", " da = da.assign_coords(longitude=(((da.longitude + 180 ) % 360 ) - 180))\n", " return da.sortby(da['longitude'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Read in the desired variables; subset dimensions; and re-order longitudes." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Modify the next cell so you are choosing only those parameters you need!
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Data variable selection-Modify depending on what variables you need\n", "\n", "#%time Q = reorder_longitudes(ds['specific_humidity'].sel(time=dateList,level=pLevel,latitude=latRange,longitude=lonRange).compute())\n", "%time U = reorder_longitudes(ds['u_component_of_wind'].sel(time=dateList,level=pLevel,latitude=latRange,longitude=lonRange).compute())\n", "%time V = reorder_longitudes(ds['v_component_of_wind'].sel(time=dateList,level=pLevel, latitude=latRange,longitude=lonRange).compute())\n", "#%time W = reorder_longitudes(ds['vertical_velocity'].sel(time=dateList,level=pLevel, latitude=latRange,longitude=lonRange).compute())\n", "%time T = reorder_longitudes(ds['temperature'].sel(time=dateList,level=pLevel, latitude=latRange,longitude=lonRange).compute())\n", "%time Z = reorder_longitudes(ds['geopotential'].sel(time=dateList,level=pLevel,latitude=latRange,longitude=lonRange).compute())\n", "#%time SLP = reorder_longitudes(ds['mean_sea_level_pressure'].sel(time=dateList,latitude=latRange,longitude=lonRange).compute())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
%time magic directive: This Jupyter lab directive will display how long each line beginning with the directive takes to execute.
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set a variable corresponding to the number of times in the desired time range" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "nTimes = len(dateList)\n", "nTimes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Examine one of the data arrays that you have read in" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Z" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Perform unit conversions: modify depending on what variables you have read in!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Z = mpcalc.geopotential_to_height(Z).metpy.convert_units('dam')\n", "U = U.metpy.convert_units('kts')\n", "V = V.metpy.convert_units('kts')\n", "T = T.metpy.convert_units('degC')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Depending on what variables you have read in, set appropriate contour value ranges. The cells below are for temperature and geopotential height. Modify to fit the context of what you are visualizing!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set contour levels for T, after first inspecting its range of values." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "T.min(),T.max()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "tLevels = np.arange(-45,39,3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the purposes of plotting geopotential heights in decameters, choose an appropriate contour interval and range of values ... for geopotential heights, a common convention is: from surface up through 700 hPa: 3 dam; above that, 6 dam to 400 and then 9 or 12 dam from 400 and above." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "if (pLevel == 1000):\n", " zLevels= np.arange(-90,63, 3)\n", "elif (pLevel == 850):\n", " zLevels = np.arange(60, 183, 3)\n", "elif (pLevel == 700):\n", " zLevels = np.arange(201, 339, 3)\n", "elif (pLevel == 500):\n", " zLevels = np.arange(468, 606, 6)\n", "elif (pLevel == 300):\n", " zLevels = np.arange(765, 1008, 9)\n", "elif (pLevel == 200): \n", " zLevels = np.arange(999, 1305, 9)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create objects for the relevant coordinate arrays; in this case, *longitude*, *latitude*, and *time*.\n", "#### Modify so you are using a variable you have read in! The cell below assumes we have read in temperature as T." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lons, lats, times= T.longitude, T.latitude, T.time" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Take a peek at a couple of these coordinate arrays." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lons" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Notice that the longitudinal array extends to 179.5, not 180. We will use Cartopy's add_cyclic_point method to eliminate the resulting seam.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "times" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create your plot.\n", "#### In this example, we will create contour lines of geopotential height, contour fills of temperature, and wind barbs.\n", "Let's create a plot for a single time, just to ensure all looks good." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Note: The Science on a Sphere treats titles and colorbars as separate layers. Thus, in these next cells we will not generate nor include them in the figure.
\n", "Additionally, the sphere expects its graphics to have a resolution of 2048x1024.
\n", "Finally, by default, Matplotlib includes a border frame around each Axes. We don't want that included on the sphere-projected graphic either.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "timeIdx = 0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Subset the DataArrays for that single time (thus, eliminating the time dimension)\n", "Once again, modify the next cell depending on what variables you have read in." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "Z0 = Z.isel(time=timeIdx)\n", "T0 = T.isel(time=timeIdx)\n", "U0 = U.isel(time=timeIdx)\n", "V0 = V.isel(time=timeIdx)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "proj_data = ccrs.PlateCarree() # The dataset's x- and y- coordinates are lon-lat" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "res = '110m'\n", "dpi = 100\n", "fig = plt.figure(figsize=(2048/dpi, 1024/dpi))\n", "ax = plt.subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=180), frameon=False)\n", "ax.set_global()\n", "ax.add_feature(cfeature.COASTLINE.with_scale(res))\n", "ax.add_feature(cfeature.BORDERS.with_scale(res))\n", "ax.add_feature(cfeature.STATES.with_scale(res))\n", "\n", "# Temperature (T) contour fills\n", "\n", "# Note we don't need the transform argument since the map/data projections are the same, but we'll leave it in\n", "CF = ax.contourf(lons,lats,T0,levels=tLevels,cmap=plt.get_cmap('coolwarm'), extend='both', transform=proj_data) \n", "\n", "# Height (Z) contour lines\n", "CL = ax.contour(lons,lats,Z0,zLevels,linewidths=1.25,colors='yellow', transform=proj_data)\n", "ax.clabel(CL, inline_spacing=0.2, fontsize=8, fmt='%.0f')\n", "fig.tight_layout(pad=.01)\n", "\n", "# 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.\n", "# Also need to sample (skip) a selected # of points to keep the plot readable.\n", "skip = 8\n", "ax.barbs(lons[::skip],lats[::skip],U0[::skip,::skip].values, V0[::skip,::skip].values, transform=proj_data);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hmmm, notice something? Let's now try plotting on an orthographic projection, somewhat akin to how it would look on the sphere." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "res = '110m'\n", "dpi = 100\n", "fig = plt.figure(figsize=(2048/dpi, 1024/dpi))\n", "ax = plt.subplot(1,1,1,projection=ccrs.Orthographic(central_longitude=180), frameon=False)\n", "ax.set_global()\n", "ax.add_feature(cfeature.COASTLINE.with_scale(res))\n", "ax.add_feature(cfeature.BORDERS.with_scale(res))\n", "ax.add_feature(cfeature.STATES.with_scale(res))\n", "\n", "# Temperature (T) contour fills\n", "\n", "# Note we don't need the transform argument since the map/data projections are the same, but we'll leave it in\n", "CF = ax.contourf(lons,lats,T0,levels=tLevels,cmap=plt.get_cmap('coolwarm'), extend='both', transform=proj_data) \n", "\n", "# Height (Z) contour lines\n", "CL = ax.contour(lons,lats,Z0,zLevels,linewidths=1.25,colors='yellow', transform=proj_data)\n", "ax.clabel(CL, inline_spacing=0.2, fontsize=8, fmt='%.0f')\n", "fig.tight_layout(pad=.01)\n", "\n", "# 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.\n", "# Also need to sample (skip) a selected # of points to keep the plot readable.\n", "skip = 8\n", "ax.barbs(lons[::skip],lats[::skip],U0[::skip,::skip].values, V0[::skip,::skip].values, transform=proj_data);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There \"seams\" to be a problem close to the dateline. This is because the longitude ends at 179.5 ° E ... leaving a gap between 179.5 and 180 °. Fortunately, Cartopy includes the `add_cyclic_point` method in its `utilities` class to deal with this rather common phenomenon in global gridded datasets." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "res = '110m'\n", "dpi = 100\n", "fig = plt.figure(figsize=(2048/dpi, 1024/dpi))\n", "ax = plt.subplot(1,1,1,projection=ccrs.Orthographic(central_longitude=180), frameon=False)\n", "ax.set_global()\n", "ax.add_feature(cfeature.COASTLINE.with_scale(res))\n", "ax.add_feature(cfeature.BORDERS.with_scale(res))\n", "ax.add_feature(cfeature.STATES.with_scale(res))\n", "\n", "# add cyclic points to data and longitudes\n", "# latitudes are unchanged in 1-dimension\n", "Zcyc, clons = cutil.add_cyclic_point(Z0, lons)\n", "Tcyc, clons = cutil.add_cyclic_point(T0, lons)\n", "Ucyc, clons = cutil.add_cyclic_point(U0, lons)\n", "Vcyc, clons = cutil.add_cyclic_point(V0, lons)\n", "\n", "# Temperature (T) contour fills\n", "\n", "# Note we don't need the transform argument since the map/data projections are the same, but we'll leave it in\n", "CF = ax.contourf(clons,lats,Tcyc,levels=tLevels,cmap=plt.get_cmap('coolwarm'), extend='both', transform=proj_data) \n", "\n", "# Height (Z) contour lines\n", "CL = ax.contour(clons,lats,Zcyc,zLevels,linewidths=1.25,colors='yellow', transform=proj_data)\n", "ax.clabel(CL, inline_spacing=0.2, fontsize=8, fmt='%.0f')\n", "fig.tight_layout(pad=.01)\n", "\n", "# 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.\n", "# Also need to sample (skip) a selected # of points to keep the plot readable.\n", "skip = 8\n", "\n", "# The U and V arrays produced by add_cyclic_point do not have units attached, so we do not \n", "# extract the values attribute as in the previous cells.\n", "\n", "ax.barbs(clons[::skip],lats[::skip],Ucyc[::skip,::skip], Vcyc[::skip,::skip], transform=proj_data);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
The seam has disappeared! Now we can go back to a PlateCarree projection, centered at the Greenwich Meridian.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "res = '110m'\n", "dpi = 100\n", "fig = plt.figure(figsize=(2048/dpi, 1024/dpi))\n", "ax = plt.subplot(1,1,1,projection=ccrs.PlateCarree(), frameon=False)\n", "ax.set_global()\n", "ax.add_feature(cfeature.COASTLINE.with_scale(res), edgecolor='brown', linewidth=2.5)\n", "ax.add_feature(cfeature.BORDERS.with_scale(res), edgecolor='brown', linewidth=2.5)\n", "ax.add_feature(cfeature.STATES.with_scale(res), edgecolor='brown')\n", "\n", "# add cyclic points to data and longitudes\n", "# latitudes are unchanged in 1-dimension\n", "Zcyc, clons = cutil.add_cyclic_point(Z0, lons)\n", "Tcyc, clons = cutil.add_cyclic_point(T0, lons)\n", "Ucyc, clons = cutil.add_cyclic_point(U0, lons)\n", "Vcyc, clons = cutil.add_cyclic_point(V0, lons)\n", "\n", "# Temperature (T) contour fills\n", "\n", "# Note we don't need the transform argument since the map/data projections are the same, but we'll leave it in\n", "CF = ax.contourf(clons,lats,Tcyc,levels=tLevels,cmap=plt.get_cmap('coolwarm'), extend='both', transform=proj_data) \n", "\n", "# Height (Z) contour lines\n", "CL = ax.contour(clons,lats,Zcyc,zLevels,linewidths=1.25,colors='yellow', transform=proj_data)\n", "ax.clabel(CL, inline_spacing=0.2, fontsize=8, fmt='%.0f')\n", "fig.tight_layout(pad=.01)\n", "\n", "# 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.\n", "# Also need to sample (skip) a selected # of points to keep the plot readable.\n", "skip = 8\n", "\n", "# The U and V arrays produced by add_cyclic_point do not have units attached, so we do not \n", "# extract the values attribute as in the previous cells.\n", "\n", "# Default length of wind barbs is 7. Shorten them a bit.\n", "ax.barbs(clons[::skip],lats[::skip],Ucyc[::skip,::skip], Vcyc[::skip,::skip], length = 5, transform=proj_data)\n", "# Save this figure to your current directory\n", "fig.savefig('test_ERA5_SOS.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## If you'd like, go to https://globe-3d.vercel.app/ and view how your graphic might look on the sphere." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create small (128x128) and large (800x800) thumbnails. These will serve as icons in the SoS iPad and computer apps that go along with your particular product.\n", "We'll use the orthographic projection and omit the contour lines and some of the cartographic features, and add a title string." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Note: This cell may take a while to complete. Be patient!
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "res = '110m'\n", "dpi = 100\n", "for size in (128, 800):\n", " if (size == 128):\n", " sizeStr = 'small'\n", " else:\n", " sizeStr ='big'\n", " \n", " fig = plt.figure(figsize=(size/dpi, size/dpi))\n", " ax = plt.subplot(1,1,1,projection=ccrs.Orthographic(central_longitude=-90), frameon=False)\n", " ax.set_global()\n", " ax.add_feature(cfeature.COASTLINE.with_scale(res))\n", " tl1 = caseDir\n", " tl2 = prodDir\n", " ax.set_title(f\"{tl1}\\n{tl2}\", color='purple', fontsize=8)\n", "\n", " # Contour fills\n", "\n", " CF = ax.contourf(clons,lats,Tcyc,levels=tLevels,cmap=plt.get_cmap('coolwarm'), extend='both', transform=proj_data) \n", "\n", " fig.tight_layout(pad=.01)\n", " fig.savefig(f'{mediaDir}thumbnail_{sizeStr}.jpg')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create a standalone colorbar\n", "Visualizations on the Science-on-a-Sphere consist of a series of image files, layered on top of each other. In this example, instead of having the colorbar associated with 500 hPa heights adjacent to the map, let's use a technique by which we remove the contour plot, leaving only the colorbar to be exported as an image. We change the colorbar's orientation to horizontal, and also change its tick label colors so they will be more visible on the sphere's display." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# draw a new figure and replot the colorbar there\n", "fig,ax = plt.subplots(figsize=(14,2), dpi=125)\n", "# set tick and ticklabel color\n", "tick_color='black'\n", "label_color='orange'\n", "cbar = fig.colorbar(CF, ax=ax, orientation='horizontal')\n", "cbar.ax.xaxis.set_tick_params(color=tick_color, labelcolor=label_color, labelsize=8)\n", "# Remove the Axes object ... essentially the contours and cartographic features ... from the figure.\n", "ax.remove()\n", "# All that remains is the colorbar ... save it to disk. Make the background transparent.\n", "fig.savefig(f'{labelsDir}colorbar.png',transparent=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create a set of labels for each plot" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "labelFname = f'{labelsDir}labels.txt'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define a function to create the plot for each hour. \n", "The function accepts a time array element as its argument." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def make_sos_map (timeIdx):\n", " \n", " res = '110m'\n", " dpi = 100\n", " fig = plt.figure(figsize=(2048/dpi, 1024/dpi))\n", " ax = plt.subplot(1,1,1,projection=ccrs.PlateCarree(), frameon=False)\n", " ax.set_global()\n", " ax.add_feature(cfeature.COASTLINE.with_scale(res), edgecolor='brown', linewidth=1.5)\n", " ax.add_feature(cfeature.BORDERS.with_scale(res), edgecolor='brown', linewidth=1.5)\n", " ax.add_feature(cfeature.STATES.with_scale(res), edgecolor='brown')\n", "\n", "# add cyclic points to data and longitudes\n", "# latitudes are unchanged in 1-dimension\n", " Zcyc, clons= cutil.add_cyclic_point(Z.isel(time=timeIdx), lons)\n", " Tcyc, clons = cutil.add_cyclic_point(T.isel(time=timeIdx), lons)\n", " Ucyc, clons = cutil.add_cyclic_point(U.isel(time=timeIdx), lons)\n", " Vcyc, clons = cutil.add_cyclic_point(V.isel(time=timeIdx), lons)\n", "\n", "# Temperature (T) contour fills\n", "\n", "# Note we don't need the transform argument since the map/data projections are the same, but we'll leave it in\n", " CF = ax.contourf(clons,lats,Tcyc,levels=tLevels,cmap=plt.get_cmap('coolwarm'), extend='both', transform=proj_data) \n", "\n", "# Height (Z) contour lines\n", " CL = ax.contour(clons,lats,Zcyc,zLevels,linewidths=1.25,colors='yellow', transform=proj_data)\n", " ax.clabel(CL, inline_spacing=0.2, fontsize=8, fmt='%.0f')\n", " fig.tight_layout(pad=.01)\n", "\n", "# 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.\n", "# Also need to sample (skip) a selected # of points to keep the plot readable.\n", " skip = 8\n", "# Default length of wind barbs is 7. Shorten them a bit.\n", " ax.barbs(clons[::skip],lats[::skip],Ucyc[::skip,::skip], Vcyc[::skip,::skip], length = 5, transform=proj_data)\n", " frameNum = f'{timeIdx}'.zfill(2)\n", " figName = f'{graphicsDir}ERA5_{caseDir}_{prodDir}_{frameNum}.png'\n", " fig.savefig(figName)\n", " \n", " # Reduce the size of the PNG image via the Linux pngquant utility. The -f option overwites the resulting file if it already exists.\n", " # The output file will end in \"-fs8.png\"\n", " ! pngquant -f {figName}\n", " \n", " # Remove the original PNG\n", " ! rm -f {figName}\n", " \n", " # Do not show the graphic in the notebook\n", " plt.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create the graphics and titles\n", "- Open a handle to the labels file\n", "- Define the time dimension index's start, end, and increment values\n", "- Loop over the period of interest\n", " - Perform any necessary unit conversions\n", " - Create each timestep's graphic\n", " - Write the title line to the text file.\n", "- Close the handle to the labels file" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
For demonstration purposes, we will use only produce graphics for the first two timesteps ... thus only 2 graphics will be produced. When you are ready, change the nFrames to nTimes so all times are processed.
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Create a title string for the SoS labels file" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## Modify the title to fit your case!\n", "SOS_label_title = f'ERA5 {pLevel} hPa Z/T/Wind'" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "labelFileObject = open(labelFname, 'w')\n", "nFrames = 2\n", "#nFrames = nTimes\n", "\n", "for timeIdx in range(0, nFrames, 1):\n", " \n", " make_sos_map(timeIdx)\n", " \n", " # Construct the title string and write it to the file\n", " valTime = pd.to_datetime(times[timeIdx].values).strftime(\"%m/%d/%Y %H UTC\")\n", " tl1 = f'{SOS_label_title} {valTime} \\n' # \\n is the newline character\n", " labelFileObject.write(tl1)\n", " print(tl1)\n", " \n", "# Close the text file\n", "labelFileObject.close()\n", "print(\"Finished with SoS graphic loop\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create an SoS playlist file\n", "We follow the guidelines in https://sos.noaa.gov/support/sos/manuals/datasets/playlists/#dataset-playlist.sos-files" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "playlistFname = f'{playlistDir}playlist.sos'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Open a file handle to the playlist file\n", "- Add each line of the playlist. *Modify as needed for your case and product; in general you will only need to modify the **creaName** value!*\n", "- Close the file handle " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "plFileObject = open(playlistFname, 'w')\n", "subCat = \"ERA5\"\n", "cRet = \"\\n\" # New line character code\n", "\n", "creaName = \"Rovin Lazyle\" # Put your name here!\n", "\n", "plFileObject.write(f\"name = {SOS_label_title}{cRet}\")\n", "plFileObject.write(f\"description = {SOS_label_title}{cRet}\")\n", "\n", "plFileObject.write(f\"pip = ../labels/colorbar.png{cRet}\")\n", "plFileObject.write(f\"pipheight = 10{cRet}\")\n", "plFileObject.write(f\"pipvertical = -35{cRet}\")\n", "\n", "plFileObject.write(f\"label = ../labels/labels.txt{cRet}\")\n", "plFileObject.write(f\"layer = Grids{cRet}\")\n", "plFileObject.write(f\"layerdata = ../2048{cRet}\")\n", "\n", "plFileObject.write(f\"firstdwell = 2000{cRet}\")\n", "plFileObject.write(f\"lastdwell = 3000{cRet}\")\n", "plFileObject.write(f\"fps = 8{cRet}\")\n", "\n", "plFileObject.write(f\"zrotationenable = 1{cRet}\")\n", "plFileObject.write(f\"zfps = 30{cRet}\")\n", "\n", "plFileObject.write(f\"subcategory = {subCat}{cRet}\")\n", "plFileObject.write(f\"Creator = {creaName}{cRet}\")\n", " \n", "plFileObject.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Display the contents of the playlist file" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "! cat {playlistFname}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
The playlist essentially describes the following:\n", "
  1. The name and description of the product, which are used in the SoS iPad and computer apps
    2. \n", "
  2. The three components of what is displayed on the screen:
    3. \n", "
    1. The colorbar (a picture-in-picture, aka pip), with its height and vertical position
      2. \n", "
    2. The label, or title, which describes the graphic
      3. \n", "
    3. The graphic layer itself. Multiple graphic layers, each with its own layer name and directory path, could be included; in this case, there is just one.
    \n", "
  3. The dwell times, in ms, for the first frame, last frame, and all frames in-between.
\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### We're done! The directory tree you have created can then be copied/synced to the correct directory on the SoS computer (Ross and Kevin will take care of this)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## Summary\n", "\n", "We now have an end-to-end workflow that will create all that is necessary for a custom SoS visualization, using ERA5 reanalysis data.\n", "\n", "### What's Next?\n", "\n", "Use this notebook as a template for your own case and its accompanying visualizations." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 Jan. 2025 Environment", "language": "python", "name": "jan25" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.8" } }, "nbformat": 4, "nbformat_minor": 4 }