{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Rasters 2: Overlays" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Overview\n", "\n", "\n", "1. Read in a raster file in GeoTIFF format of surface elevation using the Rioxarray package\n", "1. Explore the attributes of the GeoTIFF file\n", "1. Plot the GeoTIFF using Matplotlib and Cartopy\n", "1. Overlay the GeoTIFF on a plot of realtime forecast reflectivity data\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Prerequisites\n", "\n", "| Concepts | Importance | Notes |\n", "| --- | --- | --- |\n", "| Rasters 1| Necessary | |\n", "\n", "* **Time to learn**: 30 minutes\n", "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
WARNING: due to the high-resolutions of both the raster image and the HRRR forecast reflectivity, this notebook will use approx. 8GB of system memory as you run it! You may wish to simply use this already-populated notebook as a reference, and only modify/run through it when you need to. \n", "
Be sure to close and shutdown the notebook when finished!\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Imports" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from datetime import datetime, timedelta\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import cartopy.crs as ccrs\n", "import cartopy.feature as cfeat\n", "import numpy as np\n", "import xarray as xr\n", "\n", "from metpy.plots import USCOUNTIES\n", "from metpy.plots import ctables\n", "import rioxarray as rio\n", "from matplotlib.colors import BoundaryNorm\n", "from matplotlib.ticker import MaxNLocator" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Note: Besides the rioxarray library, we've also imported USCOUNTIES from MetPy. This allows for the plotting of US county line boundaries, similar to other Cartopy cartographic features such as states.
\n", "Additonally, we've imported a set of color tables, also from MetPy.
\n", "We've also imported a couple additional methods from Matplotlib, which we'll use to customize the look of our plot when using the imshow method.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read in a raster file in GeoTIFF format" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Process the topography GeoTIFF file, downloaded from Natural Earth, using the [rioxarray](https://corteva.github.io/rioxarray/stable/#) package. This package provides an xarray interface while also presenting functionality similar to the [rasterio](https://rasterio.readthedocs.io/) library." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "msr = rio.open_rasterio('/spare11/atm533/data/raster/US_MSR_10M/US_MSR.tif')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Explore the GeoTIFF file. Remember, it's a TIFF file that is georeferenced." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "msr" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Note: \n", " rioxarray returns a DataArray that has a unique coordinate variable, spatial_ref. This coordinate variable's attributes includes CRS info in what's called well-known text (WKT) format, as well as the Affine transformation parameters. See references at the end of this notebook.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Get the bounding box coordinates of the raster file. We will need to pass them to Matplotlib's `imshow` function." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "left,bottom,right,top = msr.rio.bounds()\n", "left,bottom,right,top" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These values represent the coordinates, in units specific to the GeoTIFF's geographic coordinate reference system, of the lower-left and upper-right corners of the raster file." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Get the projection info, which we can later use in Cartopy. A rioxarray dataset object has an attribute called `crs` for this purpose." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "msr.rio.crs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Assign this CRS to a Cartopy projection object." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "projMSR = ccrs.epsg('3857')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Note: What does this coordinate reference system refer to? It is called Web Mercator, aka Pseudo Mercator, which is what is used in most web-based interactive mapping applications, such as Google Earth. See https://epsg.io/3857 " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Examine the `spatial_ref` attribute of this Dataset." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "msr.spatial_ref" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This information makes it possible for GIS applications, such as Cartopy, to transform data that's represented in one CRS to another." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Select the first (and only) band and redefine the DataArray object, so it only has dimensions of x and y." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "msr = msr.sel(band=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Determine the min and max values. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "msr.min(),msr.max()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Note: For this raster file, the values (which were manually determined by its author) are inversely proportional to elevation. Thus, lower values correspond to higher elevation.
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The next cell will determine how the topography field will be displayed on the figure. These settings work well for this dataset, but obviously must be customized for use with other raster images. This is akin to setting contour fill intervals, but instead of using Matplotlib's `contourf` method, we will use `imshow`, and take advantage of the two Matplotlib methods we imported at the start of the notebook." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Pick a binning (which is analagous to the # of fill contours, and set a min/max range)\n", "levels = MaxNLocator(nbins=30).tick_values(msr.min(), msr.max()-20)\n", "\n", "# Pick the desired colormap, sensible levels, and define a normalization\n", "# instance which takes data values and translates those into levels.\n", "# In this case, we take a color map of 256 discrete colors and normalize them down to 30.\n", "cmap = plt.get_cmap('gist_gray')\n", "norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Examine the levels." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "levels" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot the GeoTIFF using Matplotlib and Cartopy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create a region over which we will plot ... in this case, centered over NY State." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "latN_sub = 45.2\n", "latS_sub = 40.2\n", "lonW_sub = -80.0\n", "lonE_sub = -71.5\n", "cLat_sub = (latN_sub + latS_sub)/2.\n", "cLon_sub = (lonW_sub + lonE_sub )/2.\n", "\n", "proj_sub = ccrs.LambertConformal(central_longitude=cLon_sub,\n", " central_latitude=cLat_sub,\n", " standard_parallels=[cLat_sub])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the GeoTIFF using `imshow`. Just as we've previously transformed data from PlateCarree (lat/lon), here we transform from the Web Mercator projection that the elevation raster file uses." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "res = '50m' # Use the medium resolution Natural Earth shapefiles so plotting goes a bit faster than with the highest (10m)\n", "fig = plt.figure(figsize=(18,12))\n", "ax = plt.subplot(1,1,1,projection=proj_sub)\n", "ax.set_extent((lonW_sub,lonE_sub,latS_sub,latN_sub),crs=ccrs.PlateCarree())\n", "ax.add_feature (cfeat.LAND.with_scale(res))\n", "ax.add_feature (cfeat.OCEAN.with_scale(res))\n", "ax.add_feature(cfeat.COASTLINE.with_scale(res))\n", "ax.add_feature (cfeat.LAKES.with_scale(res), alpha = 0.5)\n", "ax.add_feature (cfeat.STATES.with_scale(res))\n", "ax.add_feature(USCOUNTIES,edgecolor='grey', linewidth=1 );\n", "im = ax.imshow(msr,cmap=plt.get_cmap('gist_gray'),\n", " alpha=0.4,norm=norm,transform=projMSR, zorder=10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
WTF! Well ... that does not seem to accurately represent the topography in the NE US!\n", " One additional parameter we must pass into imshow is extent. The default values are not what we want; instead, we use the Bounding Box values that we read in earlier. See the documentation for extent in imshow's documentation.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "res = '50m'\n", "fig = plt.figure(figsize=(18,12))\n", "ax = plt.subplot(1,1,1,projection=proj_sub)\n", "ax.set_extent((lonW_sub,lonE_sub,latS_sub,latN_sub),crs=ccrs.PlateCarree())\n", "ax.add_feature (cfeat.LAND.with_scale(res))\n", "ax.add_feature (cfeat.OCEAN.with_scale(res))\n", "ax.add_feature(cfeat.COASTLINE.with_scale(res))\n", "ax.add_feature (cfeat.LAKES.with_scale(res), alpha = 0.5)\n", "ax.add_feature (cfeat.STATES.with_scale(res))\n", "ax.add_feature(USCOUNTIES,edgecolor='grey', linewidth=1)\n", "# Pass in the bounding box coordinates as the value for extent\n", "im = ax.imshow(msr, extent = [left,right,bottom,top],cmap=plt.get_cmap('gist_gray'),\n", " alpha=0.4,norm=norm,transform=projMSR);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's more like it! Now, we can plot maps that include one (or more) GeoTIFF raster images, while also overlaying meteorological data, such as point observations or gridded model output." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Overlay the GeoTIFF on a plot of realtime forecast reflectivity data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below, we go through the same steps that we did in the first Raster notebook to gather the latest HRRR model output, but plot filled contours of forecasted 1km above ground level radar reflectivity over the background of topography. We'll plot over the CONUS." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "now = datetime.now()\n", "year = now.year\n", "month = now.month\n", "day = now.day\n", "hour = now.hour\n", "minute = now.minute\n", "print (year, month, day, hour,minute)\n", "\n", "if (minute < 45): \n", " hrDelta = 3\n", "else:\n", " hrDelta = 2\n", "runTime = now - timedelta(hours=hrDelta)\n", "runTimeStr = runTime.strftime('%Y%m%d %H00 UTC')\n", "modelDate = runTime.strftime('%Y%m%d')\n", "modelHour = runTime.strftime('%H')\n", "modelDay = runTime.strftime('%D')\n", "print (modelDay)\n", "print (modelHour)\n", "print (runTimeStr)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "data_url = 'https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/HRRR/CONUS_2p5km/HRRR_CONUS_2p5km_'+ modelDate + '_' + modelHour + '00.grib2'" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ds = xr.open_dataset(data_url)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "refl = ds.Reflectivity_height_above_ground" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "refl" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Select the time index (corresponding to the forecast hour; index 0 corresponds to the analysis time) and the (single) height index. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "idxTime = 0\n", "idxHght = 0\n", "validTime = refl.time1.isel(time1 = idxTime).values # Returns validTime as a NumPy Datetime 64 object\n", "validTime = pd.Timestamp(validTime) # Converts validTime to a Pandas Timestamp object\n", "timeStr = validTime.strftime(\"%Y-%m-%d %H%M UTC\") # We can now use Datetime's strftime method to create a time string.\n", "tl1 = \"HRRR 1km AGL Reflectivity (dBZ)\"\n", "tl2 = str('Valid at: '+ timeStr)\n", "title_line = (tl1 + '\\n' + tl2 + '\\n')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Note: The names of some of the coordinate variables, particularly the vertical and time ones, can vary from run-to-run. Although MetPy's Xarray accessor has a clever means of dealing with that, for now, pay attention to the names of the coordinate variables and substitute when necessary. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "refl = ds.Reflectivity_height_above_ground.isel(time1=idxTime, height_above_ground1=idxHght)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Get the projection info from the HRRR Dataset" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "proj = ds.LambertConformal_Projection\n", "\n", "lat1 = proj.latitude_of_projection_origin\n", "lon1 = proj.longitude_of_central_meridian\n", "slat = proj.standard_parallel\n", "projHRRR = ccrs.LambertConformal(central_longitude=lon1,\n", " central_latitude=lat1,\n", " standard_parallels=[slat])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use a set of color tables that mimics what the National Weather Services uses for its reflectivity plots." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ref_norm, ref_cmap = ctables.registry.get_with_steps('NWSReflectivity', 5, 5)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "refl_range = np.arange(5,80,5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Be units-aware!\n", "Since we will be using Matplotlib's `contourf` method, we need to pass in the x- and y- coordinates. In the HRRR dataset we obtained from the THREDDS server, these coordinates are in units of kilometers. For a projection such as Lambert Conformal Conic, the units need to be in meters. We could simply mulptiply the x- and y- coordinate arrays by 1000, but in this case, let's have MetPy do the unit translation for us.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = refl.x.metpy.convert_units('m')\n", "y = refl.y.metpy.convert_units('m')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make the map. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Note:\n", "
  1. For speed, choose the coarsest resolution for the cartographic features. It's fine for a larger view such as the CONUS.
    2. \n", "
  2. For readability, do not plot county outlines for a regional extent as large as the CONUS
    3. \n", "
  3. This step is particularly memory-intensive, since we are reading in two high-resolution rasters (elevation and reflectivity)!\n", "
\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "res = '110m'\n", "fig = plt.figure(figsize=(18,12))\n", "ax = plt.subplot(1,1,1,projection=projHRRR)\n", "ax.set_extent((-125,-68,20,52),crs=ccrs.PlateCarree())\n", "ax.add_feature (cfeat.LAND.with_scale(res))\n", "ax.add_feature (cfeat.OCEAN.with_scale(res))\n", "ax.add_feature(cfeat.COASTLINE.with_scale(res))\n", "ax.add_feature (cfeat.LAKES.with_scale(res), alpha = 0.5)\n", "ax.add_feature (cfeat.STATES.with_scale(res))\n", "# Usually avoid plotting counties once the regional extent is large\n", "#ax.add_feature(USCOUNTIES,edgecolor='grey', linewidth=1 );\n", "im = ax.imshow(msr,extent=[left, right, bottom, top],cmap=plt.get_cmap('gist_gray'),\n", " alpha=0.4,norm=norm,transform=projMSR)\n", "CF = ax.contourf(x,y,refl,levels=refl_range,cmap=ref_cmap,alpha=0.5,transform=projHRRR)\n", "cbar = plt.colorbar(CF,fraction=0.046, pad=0.03,shrink=0.5)\n", "cbar.ax.tick_params(labelsize=10)\n", "cbar.ax.set_ylabel(\"Reflectivity (dBZ)\",fontsize=10);\n", "title = ax.set_title(title_line,fontsize=16)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## Summary\n", "* The Rioxarray package allows for raster-specific methods on geo-referenced Xarray objects. \n", "* When using `imshow` with Cartopy's geo-referenced axes in Matplotlib, one must specify the `extent`, derived from the GeoTIFF's bounding box.\n", "* GeoTIFFs can be overlayed with other geo-referenced datasets, such as gridded fields or point observations.\n", "* The higher the resolution, the higher the memory usage.\n", "\n", "### What's Next?\n", "In the next notebook, we'll take a look at a multi-band GeoTIFF that is optimized for cloud storage." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n", "1. [US Counties in MetPy](https://unidata.github.io/MetPy/latest/examples/plots/US_Counties.html?highlight=uscounties)\n", "1. [CRS definitions](https://www.earthdatascience.org/courses/use-data-open-source-python/intro-vector-data-python/spatial-data-vector-shapefiles/epsg-proj4-coordinate-reference-system-formats-python)\n", "1. WKT definitions:
\n", " a. [From opengeospatial.org](https://docs.opengeospatial.org/is/18-010r7/18-010r7.html)
\n", " b. [From Wikipedia](https://en.wikipedia.org/wiki/Well-known_text_representation_of_coordinate_reference_systems)\n", "1. Affine transformation defs:
\n", " a. [From Wikipedia](https://en.wikipedia.org/wiki/Affine_transformation)
\n", " b. [From Univ of Minnesota GIS Course](https://giscourses.cfans.umn.edu/sites/giscourses.cfans.umn.edu/files/5480chapter4_tranf_excerpt.pdf)\n", "1. [imshow's documentation](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.imshow.html)" ] } ], "metadata": { "anaconda-cloud": {}, "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 }