Rasters 2: Overlays¶

Overview¶

Read in a raster file in GeoTIFF format of surface elevation using the Rioxarray package
Explore the attributes of the GeoTIFF file
Plot the GeoTIFF using Matplotlib and Cartopy
Overlay the GeoTIFF on a plot of realtime forecast reflectivity data

Prerequisites¶

Concepts	Importance	Notes
Rasters 1	Necessary

Time to learn: 30 minutes

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.
Be sure to close and shutdown the notebook when finished!

Imports¶

from datetime import datetime, timedelta
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeat
import numpy as np
import xarray as xr

from metpy.plots import USCOUNTIES
from metpy.plots import ctables
import rioxarray as rio
from matplotlib.colors import BoundaryNorm
from matplotlib.ticker import MaxNLocator

ERROR 1: PROJ: proj_create_from_database: Open of /knight/anaconda_aug22/envs/aug22_env/share/proj failed

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.
Additonally, we've imported a set of color tables, also from MetPy.
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.

Read in a raster file in GeoTIFF format¶

Process the topography GeoTIFF file, downloaded from Natural Earth, using the rioxarray package. This package provides an xarray interface while also presenting functionality similar to the rasterio library.

msr = rio.open_rasterio('/spare11/atm533/data/raster/US_MSR_10M/US_MSR.tif')

Explore the GeoTIFF file. Remember, it’s a TIFF file that is georeferenced.

Note: 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.

Get the bounding box coordinates of the raster file. We will need to pass them to Matplotlib’s imshow function.

left,bottom,right,top = msr.rio.bounds()
left,bottom,right,top

(-14920200.0, 666726.142827049, -5787726.0672549475, 7560000.0)

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.

Get the projection info, which we can later use in Cartopy. A rioxarray dataset object has an attribute called crs for this purpose.

msr.rio.crs

CRS.from_epsg(3857)

Assign this CRS to a Cartopy projection object.

projMSR = ccrs.epsg('3857')

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

Examine the spatial_ref attribute of this Dataset.

msr.spatial_ref

<xarray.DataArray 'spatial_ref' ()>
array(0)
Coordinates:
    spatial_ref  int64 0
Attributes:
    crs_wkt:       PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["W...
    spatial_ref:   PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["W...
    GeoTransform:  -14920200.0 726.0672549487241 0.0 7560000.0 0.0 -726.14282...

This information makes it possible for GIS applications, such as Cartopy, to transform data that’s represented in one CRS to another.

Select the first (and only) band and redefine the DataArray object, so it only has dimensions of x and y.

msr = msr.sel(band=1)

Determine the min and max values.

msr.min(),msr.max()

(<xarray.DataArray ()>
 array(113, dtype=uint8)
 Coordinates:
     band         int64 1
     spatial_ref  int64 0,
 <xarray.DataArray ()>
 array(249, dtype=uint8)
 Coordinates:
     band         int64 1
     spatial_ref  int64 0)

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.

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.

# Pick a binning (which is analagous to the # of fill contours, and set a min/max range)
levels = MaxNLocator(nbins=30).tick_values(msr.min(), msr.max()-20)

# Pick the desired colormap, sensible levels, and define a normalization
# instance which takes data values and translates those into levels.
# In this case, we take a color map of 256 discrete colors and normalize them down to 30.
cmap = plt.get_cmap('gist_gray')
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)

Examine the levels.

levels

array([112., 116., 120., 124., 128., 132., 136., 140., 144., 148., 152.,
       156., 160., 164., 168., 172., 176., 180., 184., 188., 192., 196.,
       200., 204., 208., 212., 216., 220., 224., 228., 232.])

Plot the GeoTIFF using Matplotlib and Cartopy¶

Create a region over which we will plot … in this case, centered over NY State.

latN_sub = 45.2
latS_sub = 40.2
lonW_sub = -80.0
lonE_sub = -71.5
cLat_sub = (latN_sub + latS_sub)/2.
cLon_sub = (lonW_sub + lonE_sub )/2.

proj_sub = ccrs.LambertConformal(central_longitude=cLon_sub,
                             central_latitude=cLat_sub,
                             standard_parallels=[cLat_sub])

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.

res = '50m' # Use the medium resolution Natural Earth shapefiles so plotting goes a bit faster than with the highest (10m)
fig = plt.figure(figsize=(18,12))
ax = plt.subplot(1,1,1,projection=proj_sub)
ax.set_extent((lonW_sub,lonE_sub,latS_sub,latN_sub),crs=ccrs.PlateCarree())
ax.add_feature (cfeat.LAND.with_scale(res))
ax.add_feature (cfeat.OCEAN.with_scale(res))
ax.add_feature(cfeat.COASTLINE.with_scale(res))
ax.add_feature (cfeat.LAKES.with_scale(res), alpha = 0.5)
ax.add_feature (cfeat.STATES.with_scale(res))
ax.add_feature(USCOUNTIES,edgecolor='grey', linewidth=1 );
im = ax.imshow(msr,cmap=plt.get_cmap('gist_gray'),
                 alpha=0.4,norm=norm,transform=projMSR, zorder=10)

WTF! Well ... that does not seem to accurately represent the topography in the NE US! 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.

res = '50m'
fig = plt.figure(figsize=(18,12))
ax = plt.subplot(1,1,1,projection=proj_sub)
ax.set_extent((lonW_sub,lonE_sub,latS_sub,latN_sub),crs=ccrs.PlateCarree())
ax.add_feature (cfeat.LAND.with_scale(res))
ax.add_feature (cfeat.OCEAN.with_scale(res))
ax.add_feature(cfeat.COASTLINE.with_scale(res))
ax.add_feature (cfeat.LAKES.with_scale(res), alpha = 0.5)
ax.add_feature (cfeat.STATES.with_scale(res))
ax.add_feature(USCOUNTIES,edgecolor='grey', linewidth=1)
# Pass in the bounding box coordinates as the value for extent
im = ax.imshow(msr, extent = [left,right,bottom,top],cmap=plt.get_cmap('gist_gray'),
                 alpha=0.4,norm=norm,transform=projMSR);

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.

ATM433/533 Fall 2022

Rasters 2: Overlays

Contents

Rasters 2: Overlays¶

Overview¶

Prerequisites¶

Imports¶

Read in a raster file in GeoTIFF format¶

Plot the GeoTIFF using Matplotlib and Cartopy¶

Overlay the GeoTIFF on a plot of realtime forecast reflectivity data¶

Summary¶

What’s Next?¶

References¶