Interactive plots with Ipyleaflet

Overview:

  1. Interactively display images from remote map servers with ipyleaflet

  2. Overlay gridded data on the interactive map

Older Enivronment Note that we are using an older Python 3 environment in this notebook! The ipyleaflet package is currently incompatible with the version of Jupyterlab that we have been using this semester.

Prerequisites

Concepts

Importance

Notes

Contextily

Xarray

  • Time to learn: 20 minutes


Imports

from datetime import datetime, timedelta
from ipyleaflet import Map, TileLayer, basemaps, basemap_to_tiles
from ipyleaflet.velocity import Velocity
import xarray as xr

Interactively display images from remote map servers with ipyleaflet

While we’ve made some nice looking maps so far in the course, they are all static images … i.e., they can only be viewed … not interacted with. Most of us have browsed interactive map-based websites, such as Google Maps, where we can use our mouse or touch-enabled screen to move around, zoom in or out, or otherwise intereact with the plot. The Javascript programming language is what makes this interactivity possible. Here, we will use a Python package, ipyleaflet, as a “bridge” to Javascript.

Note: The "leaflet" in ipyleaflet refers to the Leaflet open-source Javascript library.

Use a satellite composite as the base image.

Now, specify an image tile server and product, using the same sources as those provided by Contextily. Center the map and specify a zoom level.

Create time strings

We’ll need to select a date and time, and then pass in date/time objects or strings to various functions below; their format will be function-specific.

timeObj = datetime(2024,4,25,0)
timeStr1 = timeObj.strftime("%Y-%m-%d")
timeStr2 = timeObj.strftime("%Y%m%d_%H%M")
timeStr2
'20240425_0000'
m = Map(
    basemap=basemap_to_tiles(basemaps.NASAGIBS.ModisTerraTrueColorCR, timeStr1),
    center=(42.204793, -74.121558),
    zoom=4
)

m

Add a background cartographic layer.

We’ll use ipyleaflet’s add_layer method. We simply add another remote tile server. The default opacity is 1.0, so reduce that a bit so we still see the satellite layer.

basemapLayer = basemap_to_tiles(basemaps.CartoDB.Positron,opacity=0.5)
m.add_layer(basemapLayer)

Interact with the dynamic map

Zoom in and out with the +/- buttons to see both layers dynamically update.

Overlay gridded data on the interactive map

We will display an animated wind velocity flow field, a la https://earth.nullschool.net or https://www.windy.com/.

Note: ipyleaflet can read in Xarray DataArray objects so we will exploit that here.

Read in a dataset based on the 0.5-degree resolution GFS on a given date and time.

ds = xr.open_dataset('https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/Global_0p5deg_ana/GFS_Global_0p5deg_ana_'+ timeStr2 +'.grib2')

Examine the dataset

ds
<xarray.Dataset>
Dimensions:                                                          (lat: 361,
                                                                      lon: 720,
                                                                      time: 1,
                                                                      altitude_above_msl: 3,
                                                                      sigma_layer: 4,
                                                                      pressure_difference_layer: 1,
                                                                      ...
                                                                      sigma: 1,
                                                                      pressure_difference_layer1: 1,
                                                                      hybrid: 1,
                                                                      sigma_layer_bounds_1: 2,
                                                                      pressure_difference_layer_bounds_1: 2,
                                                                      pressure_difference_layer1_bounds_1: 2)
Coordinates: (12/15)
  * lat                                                              (lat) float32 ...
  * lon                                                              (lon) float32 ...
    reftime                                                          datetime64[ns] ...
  * time                                                             (time) datetime64[ns] ...
  * altitude_above_msl                                               (altitude_above_msl) float32 ...
  * sigma_layer                                                      (sigma_layer) float32 ...
    ...                                                               ...
  * height_above_ground                                              (height_above_ground) float32 ...
  * potential_vorticity_surface                                      (potential_vorticity_surface) float32 ...
  * height_above_ground1                                             (height_above_ground1) float32 ...
  * sigma                                                            (sigma) float32 ...
  * pressure_difference_layer1                                       (pressure_difference_layer1) float32 ...
  * hybrid                                                           (hybrid) float32 ...
Dimensions without coordinates: sigma_layer_bounds_1,
                                pressure_difference_layer_bounds_1,
                                pressure_difference_layer1_bounds_1
Data variables: (12/83)
    LatLon_Projection                                                int32 ...
    sigma_layer_bounds                                               (sigma_layer, sigma_layer_bounds_1) float32 ...
    pressure_difference_layer_bounds                                 (pressure_difference_layer, pressure_difference_layer_bounds_1) float32 ...
    pressure_difference_layer1_bounds                                (pressure_difference_layer1, pressure_difference_layer1_bounds_1) float32 ...
    Absolute_vorticity_isobaric                                      (time, isobaric, lat, lon) float32 ...
    Cloud_mixing_ratio_hybrid                                        (time, hybrid, lat, lon) float32 ...
    ...                                                               ...
    v-component_of_wind_potential_vorticity_surface                  (time, potential_vorticity_surface, lat, lon) float32 ...
    v-component_of_wind_altitude_above_msl                           (time, altitude_above_msl, lat, lon) float32 ...
    v-component_of_wind_maximum_wind                                 (time, lat, lon) float32 ...
    v-component_of_wind_height_above_ground                          (time, height_above_ground1, lat, lon) float32 ...
    v-component_of_wind_tropopause                                   (time, lat, lon) float32 ...
    v-component_of_wind_sigma                                        (time, sigma, lat, lon) float32 ...
Attributes:
    Originating_or_generating_Center:                                        ...
    Originating_or_generating_Subcenter:                                     ...
    GRIB_table_version:                                                      ...
    Type_of_generating_process:                                              ...
    Analysis_or_forecast_generating_process_identifier_defined_by_originating...
    file_format:                                                             ...
    Conventions:                                                             ...
    history:                                                                 ...
    featureType:                                                             ...
    _CoordSysBuilder:                                                        ...

Specify an analysis time, and then create DataArray objects that are subset for the specified time and the 20m AGL level (this is the lowest level in the dataset)

u = ds['u-component_of_wind_height_above_ground'].sel(time=timeObj,height_above_ground1=20)
v = ds['v-component_of_wind_height_above_ground'].sel(time=timeObj,height_above_ground1=20)

Create a Dataset that contains just the wind components

ipyleaflet’s Velocity method produces an animated streamline effect. It accepts an Xarray dataset that contains just two arrays … specifically, the u- and v- components of the wind vector we have chosen for visualization.

(Perhaps there is a more efficient way of doing this, but the following cell will create a new Xarray Dataset that consists of just two arrays … u and v.)

# First create the Dataset out of the 'u' Dataarray
dsWind = u.to_dataset()

# Append v to the Dataset
dsWind['v'] = v

# Rename the u DataArray to 'u'
name_dict = {'u-component_of_wind_height_above_ground':'u'}
dsWind = dsWind.rename_vars(name_dict)

# Verify that the Dataset looks ok
dsWind
<xarray.Dataset>
Dimensions:               (lat: 361, lon: 720)
Coordinates:
  * lat                   (lat) float32 90.0 89.5 89.0 ... -89.0 -89.5 -90.0
  * lon                   (lon) float32 0.0 0.5 1.0 1.5 ... 358.5 359.0 359.5
    reftime               datetime64[ns] ...
    time                  datetime64[ns] 2024-04-25
    height_above_ground1  float32 20.0
Data variables:
    u                     (lat, lon) float32 ...
    v                     (lat, lon) float32 ...

Create the animation

Center it, specify a zoom level, and select a background map; create and add a velocity layer from the Dataset that contains u and v.

center = (40.0, -70.0)
zoom = 6
m = Map(center=center, zoom=zoom, interpolation='nearest', basemap=basemaps.CartoDB.DarkMatter)


display_options = {
    'velocityType': 'Global Wind',
    'displayPosition': 'bottomleft',
    'displayEmptyString': 'No wind data'
}
wind = Velocity(data=dsWind,
                zonal_speed='u',
                meridional_speed='v',
                latitude_dimension='lat',
                longitude_dimension='lon',
                velocity_scale=0.01,
                max_velocity=120,
                display_options=display_options)
m.add_layer(wind)

Display the map.

m

If you move your mouse/touchpad, you will see direction and speed at the nearest gridpoint displayed at the bottom left of the frame.

Things to try:
  1. Create a velocity map using gridded data from another source ... such as RDA or another datasource from the Unidata THREDDS Server.
  2. Make this notebook dynamic in terms of date so it will always provide a recent (e.g. one day ago) visualization.

Summary

  • Ipyleaflet provides Javascript-enabled interactivity in a Jupyter notebook.

  • Ipyleaflet’s Velocity function creates a flow-vector field that can be layered onto a basemap.

  • The Velocity function accepts an Xarray dataset that contains the u- and v- wind components.

What’s Next?

We’ll continue to explore interactive visualizations, next with the Holoviz suite of packages.