ECMWF logo Google logo

Interactive Visualization using `GeoViews`

Overview

A team at Google Research & Cloud are making parts of the ECMWF Reanalysis version 5 (aka ERA-5) accessible in a Analysis Ready, Cloud Optimized (aka ARCO) format.

In this notebook, we will do the following:

Access the ERA-5 ARCO catalog
Select a particular dataset and variable from the catalog
Convert the data from Gaussian to Cartesian coordinates
Plot a map of sea-level pressure contours and 2-meter temperature mesh using Geoviews.

Prerequisites

Concepts	Importance	Notes
Cartopy	Necessary
Xarray	Necessary
[Geoviews]	Necessary

Time to learn: 30 minutes

Imports

import fsspec
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import scipy.spatial
import numpy as np
import geoviews as gv
from geoviews import opts
from pyproj import Transformer
from holoviews.operation.datashader import rasterize

Access the ARCO ERA-5 catalog on Google Cloud

Let’s open the single-level-reanalysis Zarr file and save two variables from it

reanalysis = xr.open_zarr(
    'gs://gcp-public-data-arco-era5/co/single-level-reanalysis.zarr', 
    chunks={'time': 48},
    consolidated=True,
)

msl = reanalysis.msl
t2m = reanalysis.t2m

Regrid to Cartesian coordinates

These reanalyses are in their native, Guassian coordinates. We will define and use several functions to convert them to a lat-lon grid, as documented in the ARCO ERA-5 GCP example notebooks

def mirror_point_at_360(ds):
  extra_point = (
      ds.where(ds.longitude == 0, drop=True)
      .assign_coords(longitude=lambda x: x.longitude + 360)
  )
  return xr.concat([ds, extra_point], dim='values')

def build_triangulation(x, y):
  grid = np.stack([x, y], axis=1)
  return scipy.spatial.Delaunay(grid)

def interpolate(data, tri, mesh):
  indices = tri.find_simplex(mesh)
  ndim = tri.transform.shape[-1]
  T_inv = tri.transform[indices, :ndim, :]
  r = tri.transform[indices, ndim, :]
  c = np.einsum('...ij,...j', T_inv, mesh - r)
  c = np.concatenate([c, 1 - c.sum(axis=-1, keepdims=True)], axis=-1)
  result = np.einsum('...i,...i', data[:, tri.simplices[indices]], c)
  return np.where(indices == -1, np.nan, result)

Select a particular time range from the dataset.

msl93 = msl.sel(time=slice('1993-03-13T18:00:00','1993-03-14T00:00:00')).compute().pipe(mirror_point_at_360)
t2m93 = t2m.sel(time=slice('1993-03-13T18:00:00','1993-03-14T00:00:00')).compute().pipe(mirror_point_at_360)

Regrid to a lat-lon grid.

tri = build_triangulation(msl93.longitude, msl93.latitude)

longitude = np.linspace(0, 360, num=360*4+1)
latitude = np.linspace(-90, 90, num=180*4+1)

mesh = np.stack(np.meshgrid(longitude, latitude, indexing='ij'), axis=-1)

grid_mesh_msl = interpolate(msl93.values, tri, mesh)
grid_mesh_t2m = interpolate(t2m93.values, tri, mesh)

Construct an Xarray DataArray from the regridded array.

da_msl = xr.DataArray(data=grid_mesh_msl,
                dims=["time", "longitude", "latitude"],
                coords=[('time', msl93.time.data), ('longitude', longitude), ('latitude', latitude)],
                name='msl')

da_t2m = xr.DataArray(data=grid_mesh_t2m,
                dims=["time", "longitude", "latitude"],
                coords=[('time', t2m93.time.data), ('longitude', longitude), ('latitude', latitude)],
                name='t2m')

Convert to hPa and deg C.

slp = da_msl/100
t2m = da_t2m-273.15

Plot the data using geoviews

Add an interactive element by using the holoviz/geoviews libraries.

We need to choose the rendering backend that we want to use in Geoviews, in this case we are using Bokeh.

gv.extension('bokeh')

If you want the plot restricted to a specific part of the world, you can specify where, however, we need to transform the coordinates from PlateCarree to WebMercator (aka EPSG-3857), which is the default projection for Geoviews with the Bokeh backend.

lonW, lonE, latS, latN = -130, -60, 20, 55 

transformer = Transformer.from_crs(ccrs.PlateCarree(), "EPSG:3857")
lon1, lat1 = transformer.transform(lonW,latS)
lon2, lat2 = transformer.transform(lonE, latN)

Convert our Xarray datasets/data arrays into Geoviews dataset objects, specifying the different dimensions of the dataset (lat, lon, time) as well as the variable we want to plot. For this first example, we’ll just visualize the first time in the selected time range.

Note:

By default, Geoviews expects a dataset whose coordinates are lat-lon (i.e., using a PlateCarree projection) to have the longitude dimension vary from -180 to 180. In this case, the ERA-5’s longitudes range from 0 to 360.
In this case, ensure that the crs of the element declares the actual central_longitude, e.g. 180 degrees rather than the default 0 degrees.

dataset_era = gv.Dataset(slp.isel(time=0), kdims=['longitude','latitude'],vdims='msl') #create a geoviews dataset
contour = gv.project(dataset_era.to(gv.LineContours, ['longitude', 'latitude'],crs=ccrs.PlateCarree(central_longitude=180))) #create a Geoviews contour object

Plotting MSLP

cint = np.arange(900,1080,4)
(gv.tile_sources.OSM * contour).opts(
    opts.LineContours(tools=['hover'], levels=cint, frame_width=700, frame_height=400,show_legend=False, line_width=3,xlim=(lon1,lon2),ylim=(lat1,lat2)))

Interactive Visualization using GeoViews