Data Visualization with GeoViews & HvPlot - ScienceCore Climate Risks

The primary tools we’ll use for plotting come from the Holoviz family of Python libraries, principally GeoViews and hvPlot. These are largely built on top of HoloViews and support multiple backends for rendering plots (notably Bokeh for interactive visualization and Matplotlib for static, publication-quality plots).

GeoViews¶

From the GeoViews documentation:

GeoViews is a Python library that makes it easy to explore and visualize geographical, meteorological, and oceanographic datasets, such as those used in weather, climate, and remote sensing research.
GeoViews is built on the HoloViews library for building flexible visualizations of multidimensional data. GeoViews adds a family of geographic plot types based on the Cartopy library, plotted using either the Matplotlib or Bokeh packages. With GeoViews, you can now work easily and naturally with large, multidimensional geographic datasets, instantly visualizing any subset or combination of them, while always being able to access the raw data underlying any plot.

import warnings
warnings.filterwarnings('ignore')
from pathlib import Path
from pprint import pprint

import geoviews as gv
gv.extension('bokeh')
from geoviews import opts

Displaying a basemap¶

A basemap or tile layer is useful when displaying vector or raster data because it allows us to overlay the relevant geospatial data on a familar gepgraphical map as a background. The principal utility is we’ll use is gv.tile_sources. We can use the method opts to specify additional confirguration settings. Below, we use the Open Street Map (OSM) Web Map Tile Service to create the object basemap. When we display the repr for this object in the notebook cell, the Bokeh menu at right enables interactive exploration.

basemap = gv.tile_sources.OSM.opts(width=600, height=400)
basemap # When displayed, this basemap can be zoomed & panned using the menu at the right

Plotting points¶

To get started, let’s define a regular Python tuple for the longitude-latitude coordinates of Tokyo, Japan.

tokyo_lonlat = (139.692222, 35.689722)
print(tokyo_lonlat)

The class geoviews.Points accepts a list of tuples (each of the form (x, y)) & constructs a Points object that can be displayed. We can overlay the point created in the OpenStreetMap tiles from basemap using the * operator in Holoviews. We can also use geoviews.opts to set various display preferences for these points.

tokyo_point   = gv.Points([tokyo_lonlat])
point_opts = opts.Points(
                          size=48,
                          alpha=0.5,
                          color='red'
                        )
print(type(tokyo_point))

# Use Holoviews * operator to overlay plot on basemap
# Note: zoom out to see basemap (starts zoomed "all the way in")
(basemap * tokyo_point).opts(point_opts)

# to avoid starting zoomed all the way in, this zooms "all the way out"
(basemap * tokyo_point).opts(point_opts, opts.Overlay(global_extent=True))

Plotting rectangles¶

Standard way to represent rectangle (also called a bounding box) with corners
$(x_{\mathrm{min}},y_{\mathrm{min}}), (x_{\mathrm{min}},y_{\mathrm{max}}), (x_{\mathrm{max}},y_{\mathrm{min}}), (x_{\mathrm{max}},y_{\mathrm{max}})$
(1)
(assuming $x_{\mathrm{max}}>x_{\mathrm{min}}$ & $y_{\mathrm{max}}>y_{\mathrm{min}}$ ) is as a single 4-tuple
$(x_{\mathrm{min}},y_{\mathrm{min}},x_{\mathrm{max}},y_{\mathrm{max}}),$
(2)
i.e., the lower,left corner coordinates followed by the upper, right corner coordinates.
Let’s create a simple function to generate a rectangle of a given width & height given the centre coordinate.

# simple utility to make a rectangle centered at pt of width dx & height dy
def make_bbox(pt,dx,dy):
    '''Returns bounding box represented as tuple (x_lo, y_lo, x_hi, y_hi)
    given inputs pt=(x, y), width & height dx & dy respectively,
    where x_lo = x-dx/2, x_hi=x+dx/2, y_lo = y-dy/2, y_hi = y+dy/2.
    '''
    return tuple(coord+sgn*delta for sgn in (-1,+1) for coord,delta in zip(pt, (dx/2,dy/2)))

We can test the preceding function using the longitude-latitude coordinates of Marrakesh, Morocco.

# Verify that the function bounds works as intended
marrakesh_lonlat = (-7.93, 31.67)
dlon, dlat = 0.5, 0.25
marrakesh_bbox = make_bbox(marrakesh_lonlat, dlon, dlat)
print(marrakesh_bbox)

The utility geoviews.Rectangles accepts a list of bounding boxes (each described by a tuple of the form (x_min, y_min, x_max, y_max)) to plot. We can again use geoviews.opts to tailor the rectangle as needed.

rectangle = gv.Rectangles([marrakesh_bbox])
rect_opts = opts.Rectangles(
                                line_width=0,
                                alpha=0.1,
                                color='red'
                            )

We can plot a point for Marrakesh as before using geoviews.Points (customized using geoviews.opts).

marrakesh_point   = gv.Points([marrakesh_lonlat])
point_opts = opts.Points(
                          size=48,
                          alpha=0.25,
                          color='blue'
                        )

Finally, we can overlay all these features on the basemap with the options applied.

(basemap * rectangle * marrakesh_point).opts( rect_opts, point_opts )

We’ll use the approach above to visualize Areas of Interest (AOIs) when constructing search queries for NASA EarthData products. In particular, the convention of representing a bounding box by (left, lower, right, upper) ordinates is also used in the PySTAC API.

hvPlot¶

hvPlot is designed to extend the .plot API from Pandas DataFrames.
It works for Pandas DataFrames and Xarray DataArrays/Datasets.

Plotting from a DataFrame with hvplot.pandas¶

The code below loads a Pandas DataFrame of temperature data.

import pandas as pd, numpy as np
from pathlib import Path
LOCAL_PATH = Path.cwd().parent / 'assets' / 'temperature.csv'

df = pd.read_csv(LOCAL_PATH, index_col=0, parse_dates=[0])
df.head()

Reviewing the Pandas DataFrame.plot API¶

Let’s extract a subset of columns from this DataFrame and generate a plot.

west_coast = df[['Vancouver', 'Portland', 'San Francisco', 'Seattle', 'Los Angeles']]
west_coast.head()

The Pandas DataFrame .plot API provides access to a number of plotting methods. Here, we’ll use .plot.line, but a range of other options are available (e.g., .plot.area, .plot.bar, .plot.hist, .plot.scatter, etc.). This API has been repeated in several libraries due to its convenience.

west_coast.plot.line(); # This produces a static Matplotlib plot

Using the hvPlot DataFrame.hvplot API¶

By importing hvplot.pandas, a similar interactive plot can be generated. The API for .hvplot mimics that for .plot. For instance, we can generate the line plot above using .hvplot.line. In this case, the default plotting backend is Bokeh, so the plot is interactive.

import hvplot.pandas
west_coast.hvplot.line() # This produces an interactive Bokeh plot

It is possible to modify the plot to make it cleaner.

west_coast.hvplot.line(width=600, height=300, grid=True)

The hvplot API also works when chained together with other DataFrame method calls. For instance, we can resample the temperature data and compute a mean to smooth it out.

smoothed = west_coast.resample('2d').mean()
smoothed.hvplot.line(width=600, height=300, grid=True)

Plotting from a DataArray with hvplot.xarray¶

The Pandas .plot API was also extended to Xarray, i.e., for Xarray DataArray.

import xarray as xr
import hvplot.xarray
import rioxarray as rio

To start, load a local GeoTIFF file using rioxarray into an Zarray DataArray structure.

LOCAL_PATH = Path.cwd().parent / 'assets' / 'OPERA_L3_DIST-ALERT-HLS_T10TEM_20220815T185931Z_20220817T153514Z_S2A_30_v0.1_VEG-ANOM-MAX.tif'

data = rio.open_rasterio(LOCAL_PATH)
data

We do some minor reformatting to the DataArray.

data = data.squeeze() # to reduce 3D array with singleton dimension to 2D array
data = data.rename({'x':'easting', 'y':'northing'})
data

Reviewing the Xarray DataArray.plot API¶

The DataArray.plot API by default uses Matplotlib’s pcolormesh to display a 2D array stored within a DataArray. This takes a little time to render for this moderately high-resolution image.

data.plot(); # by default, uses pcolormesh

Using the hvPlot DataArray.hvplot API¶

Again, the DataArray.hvplot API mimics the DataArray.plot API; by default, it uses a subclass derived from holoviews.element.raster.Image.

plot = data.hvplot() # by default uses Image class
print(f'{type(plot)=}')
plot

The result above is an interactive plot, rendered using Bokeh.It’s a bit slow, but we can add some options to speed up rendering. It also requires cleaning up; for instance, the image is not square, the colormap does not highlight useful features, the axes are transposed, and so on.

Building up options incrementally to improve plots¶

Let’s add options to improve the image. To do this, we’ll initialize a Python dictionary image_opts to use within the call to the image method. We’ll use the dict.update method to add key-value pairs to the dictionary incrementally, each time improving the output image.

image_opts = dict(rasterize=True, dynamic=True)
pprint(image_opts)

To start, let’s make the call to hvplot.image explicit & specify the sequence of axes. & apply the options from the dictionary image_opts. We’ll use the dict-unpacking operation **image_opts each time we invoke data.hvplot.image.

plot = data.hvplot.image(x='easting', y='northing', **image_opts)
plot

Next, let’s fix the aspect ratio and image dimensions.

image_opts.update(frame_width=500, frame_height=500, aspect='equal')
pprint(image_opts)
plot = data.hvplot.image(x='easting', y='northing', **image_opts)
plot

The image colormap is a bit unpleasant; let’s change it and use transparency (alpha).

image_opts.update( cmap='hot_r', clim=(0,100), alpha=0.8 )
pprint(image_opts)
plot = data.hvplot.image(x='easting', y='northing', **image_opts)
plot

Before adding a basemap, we need to account for the coordinate system. This is stored in the GeoTIFF file and, when read using rioxarray.open_rasterio, it is available by using the attribute data.rio.crs.

crs = data.rio.crs
crs

We can use the CRS recovered above as an optional argument to hvplot.image. Notice the coordinates have changed on the axes, but the labels are wrong. We can fix that shortly.

image_opts.update(crs=crs)
pprint(image_opts)
plot = data.hvplot.image(x='easting', y='northing', **image_opts)
plot

Let’s now fix the labels. We’ll use the Holoviews/GeoViews opts system to specify these options.

label_opts = dict(title='VEG_ANOM_MAX', xlabel='Longitude (degrees)', ylabel='Latitude (degrees)')
pprint(image_opts)
pprint(label_opts)
plot = data.hvplot.image(x='easting', y='northing', **image_opts).opts(**label_opts)
plot

Let’s overlay the image on a basemap so we can see the terrain underneath.

base = gv.tile_sources.ESRI
base * plot

Finally, the white pixels are distracting; let’s filter them out using the DataArray method where.

plot = data.where(data>0).hvplot.image(x='easting', y='northing', **image_opts).opts(**label_opts)
plot * base

In this notebook, we applied some common strategies to generate plots. We’ll use these extensively in the rest of the tutorial.

2. Geospatial fundamentals

Array Manipulation with Xarray

2. Geospatial fundamentals

Constructing Advanced Visualizations