NOTE: You can find tables listing all of the various parameters in the ERA-5 datasets at
this link"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use your RDA credentials to access these individual NetCDF files. You should have previously registered with RDA, and then saved your user id and password in your home directory as .rdarc , with permissions such that only you have read/write access."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#Retrieve login credential for RDA.\n",
"from pathlib import Path\n",
"HOME = str(Path.home())\n",
"credFile = open(HOME+'/.rdarc','r')\n",
"userId, pw = credFile.read().split()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Connect to the RDA THREDDS server and point to the files of interest. Here, we'll comment out the ones we are not interested in but can always uncomment when we want them. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"session = requests.Session()\n",
"session.auth = (userId, pw)\n",
"#storeU10 = xr.backends.PydapDataStore.open(ds_urlU10, session=session)\n",
"#storeV10 = xr.backends.PydapDataStore.open(ds_urlV10, session=session)\n",
"storeSLP = xr.backends.PydapDataStore.open(ds_urlSLP, session=session)\n",
"#storeQ = xr.backends.PydapDataStore.open(ds_urlQ, session=session)\n",
"#storeT = xr.backends.PydapDataStore.open(ds_urlT, session=session)\n",
"#storeU = xr.backends.PydapDataStore.open(ds_urlU, session=session)\n",
"#storeV = xr.backends.PydapDataStore.open(ds_urlV, session=session)\n",
"storeZ = xr.backends.PydapDataStore.open(ds_urlZ, session=session)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ds_urlSLP"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#dsU10 = xr.open_dataset(storeU10)\n",
"#dsV10 = xr.open_dataset(storeV10)\n",
"dsSLP = xr.open_dataset(storeSLP)\n",
"#dsQ = xr.open_dataset(storeQ)\n",
"#dsT = xr.open_dataset(storeT)\n",
"#dsU = xr.open_dataset(storeU)\n",
"#dsV = xr.open_dataset(storeV)\n",
"dsZ = xr.open_dataset(storeZ)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Examine the datasets"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The surface fields span the entire month. Let's restrict the range so it matches the day of interest."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dsSLP = dsSLP.sel(time=slice(dayBeg,dayEnd))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dsSLP"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create an Xarray DataArray object for SLP and perform unit conversion"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"SLP = dsSLP['MSL']\n",
"SLP = SLP.metpy.convert_units('hPa')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now examine the geopotential data on isobaric surfaces."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dsZ"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To reduce the memory footprint for data on isobaric surfaces, let's request only 6 out of the 37 available levels."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dsZ = dsZ.sel(level = [200,300,500,700,850,1000])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Z = dsZ['Z']"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Z"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Set contour levels for both SLP and Z. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"slpLevels = np.arange(900,1080,4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Choose a particular isobaric surface."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"pLevel = 500"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Redefine Z so it consists only of the data at the level you just chose."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Z = Z.sel(level=pLevel)"
]
},
{
"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": {},
"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": [
"
In this notebook, we did not retrieve the U- and V- wind component grids. However, if you end up modifying this notebook to do so, uncomment the lines in the following cell.
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#print (U.units)\n",
"## Convert wind from m/s to knots; use MetPy to calculate the wind speed.\n",
"#UKts = U.metpy.convert_units('knots')\n",
"#VKts = V.metpy.convert_units('knots')\n",
"## Calculate windspeed.\n",
"#wspdKts = calc.wind_speed(U,V).to('knots')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create objects for the relevant coordinate arrays; in this case, *longitude*, *latitude*, and *time*."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"lons, lats, times= SLP.longitude, SLP.latitude, SLP.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": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"times"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create contour plots of SLP, and overlay Natural Earth coastlines and borders."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
The overall look of our plot depends on many settings. Part of the fun (and challenge) of making plots is playing around with these settings ... such as contour line/fill intervals; line styles and thicknesses, and so on. There really is no \"one size fits all\" set of methods and arguments for a meteorological visualization ... every case is different.
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we will create maps for each hour of our chosen day. We will create *for* loops to do so, iterating through each hour of the day. We will group all of our map plotting code into a couple of functions; one for a global view, and the other regional. Then we will call each function in our loop. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, create a function that will create a global map using the **Plate Carree** projection. This is the projection that is used for visualizations on [NOAA's Science on a Sphere](https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwiy3KHG-tj6AhVIKFkFHTmLAgQQFnoECAwQAQ&url=https%3A%2F%2Fsos.noaa.gov%2F&usg=AOvVaw0Jk0LClAPTR5aNngu9W3Cz), such as the one we have here in ETEC!"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def make_global_map(i, z):\n",
" \n",
" res = '110m'\n",
" fig = plt.figure(figsize=(18,12))\n",
" ax = plt.subplot(1,1,1,projection=ccrs.PlateCarree())\n",
" ax.set_global()\n",
" ax.add_feature(cfeature.COASTLINE.with_scale(res))\n",
" ax.add_feature(cfeature.STATES.with_scale(res))\n",
"\n",
" # Add the title\n",
" valTime = pd.to_datetime(SLP.time[i].values).strftime(\"%m/%d/%Y at %H UTC\")\n",
" tl1 = str(\"ERA-5 SLP and \" + str(pLevel) + \"hPa Z, \" + valTime)\n",
" plt.title(tl1,fontsize=16)\n",
"\n",
" # Height (Z) contour fills\n",
" \n",
" CF = ax.contourf(lons,lats,z,levels=zLevels,cmap=plt.get_cmap('viridis'), extend='both')\n",
" cbar = fig.colorbar(CF,shrink=0.5)\n",
" cbar.set_label(r'Heights (dam)', size='medium')\n",
" \n",
" # SLP contour lines\n",
" CL = ax.contour(lons,lats,SLP.isel(time=i),slpLevels,linewidths=1.25,colors='chocolate')\n",
" ax.clabel(CL, inline_spacing=0.2, fontsize=8, fmt='%.0f'); "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, create a regional map, focusing on our weather event of interest. First, set some values relevant to the regional map projection and extent."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"lonW = -110.0\n",
"lonE = -70.0\n",
"latS = 20.0\n",
"latN = 50.0\n",
"cLon = (lonW + lonE) / 2.\n",
"cLat = (latS + latN) / 2.\n",
"proj_map = ccrs.LambertConformal(central_longitude=cLon, central_latitude=cLat)\n",
"proj_data = ccrs.PlateCarree() # Our data is lat-lon; thus its native projection is Plate Carree.\n",
"res = '50m'\n",
"constrainLon = 3.3"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def make_regional_map(i, z):\n",
" \n",
" fig = plt.figure(figsize=(18, 12))\n",
" ax = plt.subplot(1,1,1,projection=proj_map)\n",
" ax.set_extent ([lonW + constrainLon,lonE - constrainLon,latS,latN])\n",
" ax.add_feature(cfeature.COASTLINE.with_scale(res))\n",
" ax.add_feature(cfeature.STATES.with_scale(res))\n",
"\n",
" # Add the title\n",
" valTime = pd.to_datetime(SLP.time[i].values).strftime(\"%m/%d/%Y at %H UTC\")\n",
" tl1 = str(\"ERA-5 SLP and \" + str(pLevel) + \"hPa Z, \" + valTime)\n",
" plt.title(tl1,fontsize=16)\n",
"\n",
" # Height (Z) contour fills\n",
" \n",
" CF = ax.contourf(lons,lats,z,levels=zLevels,cmap=plt.get_cmap('viridis'), extend='both', transform=proj_data)\n",
" cbar = fig.colorbar(CF,shrink=0.5)\n",
" cbar.set_label(r'Heights (dam)', size='medium')\n",
" \n",
" # SLP contour lines\n",
" CL = ax.contour(lons,lats,SLP.isel(time=i),slpLevels,linewidths=1.25,colors='chocolate', transform=proj_data)\n",
" ax.clabel(CL, inline_spacing=0.2, fontsize=8, fmt='%.0f');\n",
"\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create a `for` loop and loop through the times. In this loop, we will also pass in the geopotential grid, and convert its units from geopotential to decameters."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
In general, we ought to have been able to perform this conversion on the original, multi-hour data array. However, it appears that once again this is not possible when reading from a THREDDS data server such as RDA's; we need to operate on only one single time.
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
Note: In this example, to save time we have set the hour increment to 12 in the cell below ... so we will have two graphics per loop. When you are ready to run this on your own case, change the increment to 1 so you get 24 hours worth of maps!
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"inc = 12 # change to 1 when you are ready\n",
"for time in range(0, 24, inc):\n",
" z = mpcalc.geopotential_to_height(Z.isel(time=time)).metpy.convert_units('dam')\n",
" make_global_map(time, z)\n",
" make_regional_map(time,z) # This will take some time, be patient!\n",
"\n",
" print(time)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"### Things to try\n",
"1. Once you have run through this notebook, make a copy of it and choose another weather event. **The ERA-5 now extends all the way back to 1950!!**\n",
"1. Next, create some different visualizations, using the many variables in the ERA-5. For example:\n",
"- 850 hPa heights and temperature\n",
"- 300 hPa heights and isotachs\n",
"- Snow depth\n",
"- Accumulated precipitation \n",
"- Sea surface temperature\n",
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### \"The sky's the limit!\" (or, \"There is an endless sea of possibilities!\")\n",
"Ultimately, you and your team members will create 4 different 24-hour time loops for the weather event you have chosen."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"### What's next?\n",
"In the next notebook, we will create a similar set of graphics, but save them to disk in a format amenable for loading on our Science on a Sphere!"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"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": 4
}