{ "cells": [ { "cell_type": "markdown", "id": "adc8fe25-feea-4121-8743-e786349e0ab3", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "id": "clinical-console", "metadata": {}, "source": [ "# Plotting archived HRRR data: 1km reflectivity" ] }, { "cell_type": "markdown", "id": "a32df8af-ecc4-4db7-a211-277eb46e9133", "metadata": {}, "source": [ "## Overview\n", "1. Access archived HRRR data hosted on AWS in Zarr format\n", "2. Visualize one of the variables (1km reflectivity) at an analysis time\n", "3. Use one of **MetPy**'s customized color tables\n", "4. Create and visualize a reflectivity time loop of all forecast hours in an archived HRRR run" ] }, { "cell_type": "markdown", "id": "1ec4a6c8-6042-4e91-920f-1e709d14a8b8", "metadata": {}, "source": [ "## Prerequisites\n", "\n", "| Concepts | Importance | Notes |\n", "| --- | --- | --- |\n", "| Zarr, Dask, S3 storage | 2mT notebook| Necessary | |\n", "\n", "* **Time to learn**: 30 minutes\n", "***" ] }, { "cell_type": "markdown", "id": "860a2761-73dc-4fa3-9f09-8336802f2c55", "metadata": {}, "source": [ "## Imports" ] }, { "cell_type": "code", "execution_count": null, "id": "fluid-transfer", "metadata": {}, "outputs": [], "source": [ "import xarray as xr\n", "import s3fs\n", "import metpy\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import cartopy.crs as ccrs\n", "import cartopy.feature as cfeature\n", "from metpy.plots import colortables\n", "import pandas as pd" ] }, { "cell_type": "markdown", "id": "parallel-strike", "metadata": {}, "source": [ "## Access the Zarr-formatted data on AWS" ] }, { "cell_type": "markdown", "id": "004b16ec-c133-4662-aa88-a43bd8f69ab4", "metadata": {}, "source": [ "As in the 2m temperature notebook, create objects pointing to the HRRR S3 bucket and object of interest." ] }, { "cell_type": "markdown", "id": "behavioral-apparatus", "metadata": {}, "source": [ "
Reminder: To interactively browse the contents of this archive, go to the HRRR Zarr File browser
" ] }, { "cell_type": "markdown", "id": "927ac8c8-36d7-4b70-9773-b5b5c628b842", "metadata": {}, "source": [ "To access Zarr-formatted data stored in an S3 bucket, we follow a 3-step process:\n", "1. Create URL(s) pointing to the bucket and object(s) that contain the data we want\n", "1. Create *map(s)* to the object(s) with the **s3fs** library's `S3Map` method\n", "1. Pass the *map(s)* to Xarray's `open_dataset` or `open_mfdataset` methods, and specify `zarr` as the format, via the `engine` argument." ] }, { "cell_type": "markdown", "id": "52dded5b-f813-4177-9400-f800905f4000", "metadata": {}, "source": [ "Create the URLs" ] }, { "cell_type": "code", "execution_count": null, "id": "cc62deb0-cfd2-4720-ab84-2e0bbece8c3a", "metadata": {}, "outputs": [], "source": [ "date = '20221014'\n", "hour = '00'\n", "var = 'REFD'\n", "level = '1000m_above_ground'\n", "url1 = 's3://hrrrzarr/sfc/' + date + '/' + date + '_' + hour + 'z_anl.zarr/' + level + '/' + var + '/' + level\n", "url2 = 's3://hrrrzarr/sfc/' + date + '/' + date + '_' + hour + 'z_anl.zarr/' + level + '/' + var\n", "print(url1)\n", "print(url2)" ] }, { "cell_type": "markdown", "id": "8a286eda-418b-4b95-a0de-fceb6feb393e", "metadata": {}, "source": [ "
\n", " In this case, hrrrzarr is the S3 bucket name. 1000m_above_ground and REFD are both objects within the bucket. The former object has the 1 km reflectivity array, while the latter contains the coordinate arrays of the spatial dimensions of 1 km reflectivity (i.e., x and y).\n", "
" ] }, { "cell_type": "markdown", "id": "d661c992-5db9-40af-8ad4-fac7b28dad2f", "metadata": {}, "source": [ "Create the S3 maps from the S3 object store." ] }, { "cell_type": "code", "execution_count": null, "id": "530efd16-4779-4ace-9e71-0bbf6119056e", "metadata": {}, "outputs": [], "source": [ "fs = s3fs.S3FileSystem(anon=True)\n", "file1 = s3fs.S3Map(url1, s3=fs)\n", "file2 = s3fs.S3Map(url2, s3=fs)" ] }, { "cell_type": "markdown", "id": "d76f2f8b-059c-43bc-9bc6-0f7070c254c6", "metadata": {}, "source": [ "Use Xarray's `open_mfdataset` to create a `Dataset` from these two S3 objects." ] }, { "cell_type": "code", "execution_count": null, "id": "23db3344-cf7e-48c2-b907-039a195dc73d", "metadata": {}, "outputs": [], "source": [ "ds = xr.open_mfdataset([file1,file2], engine='zarr')" ] }, { "cell_type": "markdown", "id": "processed-importance", "metadata": {}, "source": [ "Examine the dataset." ] }, { "cell_type": "code", "execution_count": null, "id": "ceramic-letter", "metadata": {}, "outputs": [], "source": [ "ds" ] }, { "cell_type": "markdown", "id": "e45ee292-a90b-4716-894e-96932027f5ca", "metadata": {}, "source": [ "Remind ourselves of the HRRR's projection info:" ] }, { "cell_type": "markdown", "id": "worldwide-thumbnail", "metadata": {}, "source": [ "### HRRR Grid Navigation: \n", " PROJECTION: LCC \n", " ANGLES: 38.5 -97.5 38.5\n", " GRID SIZE: 1799 1059\n", " LL CORNER: 21.1381 -122.7195\n", " UR CORNER: 47.8422 -60.9168" ] }, { "cell_type": "code", "execution_count": null, "id": "embedded-lafayette", "metadata": {}, "outputs": [], "source": [ "lon1 = -97.5\n", "lat1 = 38.5\n", "slat = 38.5\n", "projData= ccrs.LambertConformal(central_longitude=lon1,\n", " central_latitude=lat1,\n", " standard_parallels=[slat,slat],globe=ccrs.Globe(semimajor_axis=6371229,\n", " semiminor_axis=6371229))" ] }, { "cell_type": "markdown", "id": "16efe72d-53ea-4f24-a5cf-ca8b46a7a2d5", "metadata": {}, "source": [ "
Note: \n", " The HRRR's projection assumes a spherical earth, whose semi-major/minor axes are both equal to 6371.229 km. We therefore need to explicitly define a Globe in Cartopy with these values.\n", "
" ] }, { "cell_type": "markdown", "id": "typical-canvas", "metadata": {}, "source": [ "Examine the dataset's coordinate variables. Each x- and y- value represents distance in meters from the central latitude and longitude." ] }, { "cell_type": "code", "execution_count": null, "id": "rotary-logic", "metadata": {}, "outputs": [], "source": [ "ds.coords" ] }, { "cell_type": "markdown", "id": "entertaining-telescope", "metadata": {}, "source": [ "Create an object pointing to the dataset's data variable." ] }, { "cell_type": "code", "execution_count": null, "id": "interior-springer", "metadata": {}, "outputs": [], "source": [ "scalar = ds.REFD" ] }, { "cell_type": "markdown", "id": "sustainable-divorce", "metadata": {}, "source": [ "When we examine the object, we see that it is a special type of `DataArray` ... a *Dask* array." ] }, { "cell_type": "code", "execution_count": null, "id": "universal-syndrome", "metadata": {}, "outputs": [], "source": [ "x = scalar.projection_x_coordinate\n", "y = scalar.projection_y_coordinate" ] }, { "cell_type": "markdown", "id": "split-trail", "metadata": {}, "source": [ "## Visualize 1km Reflectivity at an analysis time" ] }, { "cell_type": "markdown", "id": "cdd75f32-faad-47c2-bb8b-a32430b5daae", "metadata": {}, "source": [ "
\n", "We'll select one of MetPy's color maps, following one of MetPy's gallery examples.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "81275e5a-01c9-4dfb-9a1b-3317eab02046", "metadata": {}, "outputs": [], "source": [ "norm, cmap = colortables.get_with_range('NWSReflectivity', 5,80)" ] }, { "cell_type": "markdown", "id": "25b9fd8a-530b-48c1-a9e0-a0719f9c14f2", "metadata": {}, "source": [ "Examine the `cmap` object" ] }, { "cell_type": "code", "execution_count": null, "id": "0330e34d-3e90-40c4-9b9a-11c3e26c691f", "metadata": {}, "outputs": [], "source": [ "cmap" ] }, { "cell_type": "markdown", "id": "858f95a3-3be3-4483-b9d1-68eff0317df2", "metadata": {}, "source": [ "Set the color map so values under the specified range are shaded in white" ] }, { "cell_type": "code", "execution_count": null, "id": "213f06c1-2adc-497c-a277-c15ab2ef3305", "metadata": {}, "outputs": [], "source": [ "cmap.set_under('w')\n", "cmap" ] }, { "cell_type": "markdown", "id": "utility-peripheral", "metadata": {}, "source": [ "## Plot the map\n", "We'll define the plot extent to nicely encompass the HRRR's spatial domain." ] }, { "cell_type": "code", "execution_count": null, "id": "43df0c08-1538-41e0-95bd-d351f8c6a918", "metadata": { "tags": [ "remove-stderr" ] }, "outputs": [], "source": [ "latN = 50.4\n", "latS = 24.25\n", "lonW = -123.8\n", "lonE = -71.2\n", "\n", "res = '50m'\n", "\n", "fig = plt.figure(figsize=(18,12))\n", "ax = plt.subplot(1,1,1,projection=projData)\n", "ax.set_extent ([lonW,lonE,latS,latN],crs=ccrs.PlateCarree())\n", "ax.add_feature(cfeature.COASTLINE.with_scale(res))\n", "ax.add_feature(cfeature.STATES.with_scale(res))\n", "\n", "# Add the title\n", "tl1 = str('HRRR 1km Reflectivity (dBz)')\n", "tl2 = str('Analysis valid at: '+ hour + '00 UTC ' + date )\n", "plt.title(tl1+'\\n'+tl2,fontsize=16)\n", "\n", "CM = ax.pcolormesh(x, y, scalar, norm=norm, cmap=cmap)\n", "# Make a colorbar for the color mesh.\n", "cbar = fig.colorbar(CM,shrink=0.5)\n", "cbar.set_label(r'1km Reflectivity (dBz)', size='large')" ] }, { "cell_type": "markdown", "id": "4cad87c2-0ecf-4186-9fb6-94139ba4211d", "metadata": {}, "source": [ "## Now, create a loop of forecast radar reflectivity for a specific date and time in the archive." ] }, { "cell_type": "markdown", "id": "4938114f-8506-4c14-b4bf-fda07995dad9", "metadata": {}, "source": [ "The bucket/object are very similar; it's `fcst` (forecast) instead of `anl` (analysis), though!" ] }, { "cell_type": "code", "execution_count": null, "id": "64d63d14-1c4f-4fcf-bd52-244ba8cdfc27", "metadata": {}, "outputs": [], "source": [ "date = '20221013'\n", "hour = '18'\n", "var = 'REFD'\n", "level = '1000m_above_ground'\n", "url1 = 's3://hrrrzarr/sfc/' + date + '/' + date + '_' + hour + 'z_fcst.zarr/' + level + '/' + var + '/' + level\n", "url2 = 's3://hrrrzarr/sfc/' + date + '/' + date + '_' + hour + 'z_fcst.zarr/' + level + '/' + var\n", "print(url1)\n", "print(url2)" ] }, { "cell_type": "code", "execution_count": null, "id": "f7e80e54-8e22-464b-90e7-a4f3d5837c51", "metadata": {}, "outputs": [], "source": [ "fs = s3fs.S3FileSystem(anon=True)\n", "file1 = s3fs.S3Map(url1, s3=fs)\n", "file2 = s3fs.S3Map(url2, s3=fs)" ] }, { "cell_type": "code", "execution_count": null, "id": "87731208-b51a-490c-b5df-12c38a68c4b9", "metadata": {}, "outputs": [], "source": [ "ds = xr.open_mfdataset([file1,file2], engine='zarr')" ] }, { "cell_type": "code", "execution_count": null, "id": "44fa691a-bd11-4070-944e-08302d97b3fb", "metadata": {}, "outputs": [], "source": [ "ds" ] }, { "cell_type": "markdown", "id": "f75636bd-5b19-4b5b-ab84-0d4b525979c3", "metadata": {}, "source": [ "The `Dataset` is very similar to the analysis; it just has a time dimension too! What times do we have available?" ] }, { "cell_type": "code", "execution_count": null, "id": "6f74fc82-b692-4a23-a118-e076b4f94a7f", "metadata": {}, "outputs": [], "source": [ "times = ds.time\n", "times" ] }, { "cell_type": "code", "execution_count": null, "id": "f652cc3a-ebda-4e77-ad3a-ad38d769b159", "metadata": {}, "outputs": [], "source": [ "def make_regional_map(i, z):\n", " \n", " latN = 50.4\n", " latS = 24.25\n", " lonW = -123.8\n", " lonE = -71.2\n", "\n", " res = '50m'\n", "\n", " fig = plt.figure(figsize=(18,12))\n", " ax = plt.subplot(1,1,1,projection=projData)\n", " ax.set_extent ([lonW,lonE,latS,latN],crs=ccrs.PlateCarree())\n", " ax.add_feature(cfeature.COASTLINE.with_scale(res))\n", " ax.add_feature(cfeature.STATES.with_scale(res))\n", "\n", "# Add the title\n", " validTime = pd.to_datetime(z.time[i].values).strftime(\"%m/%d/%Y at %H UTC\")\n", "\n", " tl1 = str('HRRR 1km Reflectivity (dBz)')\n", " tl2 = str('Forecast valid at: '+ validTime )\n", " ax.set_title(tl1+'\\n'+tl2,fontsize=16)\n", "\n", " CM = ax.pcolormesh(z.projection_x_coordinate, z.projection_y_coordinate, z.isel(time=i), norm=norm, cmap=cmap)\n", " # Make a colorbar for the color mesh.\n", " cbar = fig.colorbar(CM,shrink=0.5)\n", " cbar.set_label(r'1km Reflectivity (dBz)', size='large')\n", " plt.show()\n", " " ] }, { "cell_type": "markdown", "id": "64aba569-83dd-459c-8a99-c2f912546047", "metadata": {}, "source": [ "In the interest of *time* (pun intended), let's just make a six-hour forecast loop. Revise as you desire!" ] }, { "cell_type": "code", "execution_count": null, "id": "868a0777-1ef3-497e-8cd0-0bff330cf3d1", "metadata": {}, "outputs": [], "source": [ "parameter = 'REFD'\n", "hrStart = 0\n", "hrStop = 6\n", "hrInc = 1\n", "for i in range (hrStart, hrStop, hrInc):\n", " make_regional_map(i, ds[parameter])\n", " " ] }, { "cell_type": "markdown", "id": "20d36311-c750-4368-b266-70dbd6e3fec7", "metadata": {}, "source": [ "---\n", "## Summary\n", "* With a minimum of changes, we can adapt our initial 2-meter temperature notebook to different HRRR variables, such as 1 km reflectivity.\n", "* **Matplotlib**'s `ax.pcolormesh` function is ideal for plotting raster fields, such as radar reflectivity.\n", "* **Metpy**'s `colortables` provide several meteorologically-relevant color tables, which can be easily incorporated and modified in **matplotlib**.\n", "\n", "## Things to try\n", "* Save your maps to disk, instead of (or in addition to) within your notebook.\n", "* Create a map whose domain covers a smaller region ... e.g. the Albany region, or another area of interest to you. If the region is small enough, consider adding [MetPy's county line](https://unidata.github.io/MetPy/latest/examples/plots/US_Counties.html) outlines.\n", "* Create a near-realtime HRRR forecast loop. The most recent HRRR tends to appear on the AWS *hrrrzarr* bucket with an approximate 90-120 minute lag. Think about how you might dynamically determine the most recently available time, no matter when you run the notebook.\n", "* Besides the single-level parameters, located in the *sfc* folder of the *hrrrzarr* bucket, explore the ones on isobaric surfaces, in the *prs* folder (these are analysis times only, no forecasts). \n" ] }, { "cell_type": "code", "execution_count": null, "id": "31b47596-87d1-488b-acb6-0566fd0f1d7d", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "execution": { "allow_errors": false }, "kernelspec": { "display_name": "Python 3 August 2022 Environment", "language": "python", "name": "aug22" }, "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.10.5" } }, "nbformat": 4, "nbformat_minor": 5 }