Contents

Interactive plots with Ipyleaflet

Interactive plots with Ipyleaflet

Overview:

  1. Interactively display images from remote map servers with ipyleaflet

  2. Overlay gridded data on the interactive map

Prerequisites

Concepts

Importance

Notes

Contextily

Xarray

  • Time to learn: 20 minutes

Imports

from datetime import datetime
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.

First, create an ipyleaflet Map object with default options. Zoom out with the “-” tool button several times until you see some geographic features. Then, interactively pan around with your mouse/touchpad . At what lat-lon coordinate do you think the map is initially centered?

Note: A bug in the current release of ipyleaflet requires that we specify a zoomlevel. It will likely be fixed in an upcoming release (see the relevant pull request). 
m = Map(zoom=10)
m
Did you figure it out? The default center point is at the equator, on the Greenwich meridian. In geospatial lingo, this is termed Null Island.

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(2022,9,28,18)
timeStr1 = timeObj.strftime("%Y-%m-%d")
timeStr2 = timeObj.strftime("%Y%m%d_%H%M")

timeStr2

'20220928_1800'

m = Map(
    basemap=basemap_to_tiles(basemaps.NASAGIBS.ModisTerraTrueColorCR, timeStr1),
    center=(28.5, -85),
    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)
# After executing this cell, note that the map has been updated with this layer.

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('/spare11/atm533/data/2022092818_gfs.nc')

Examine the dataset

ds

<xarray.Dataset>
Dimensions:                                                                   (
                                                                               lat: 361,
                                                                               lon: 720,
                                                                               time: 1,
                                                                               height_above_ground: 1,
                                                                               pressure_difference_layer: 1,
                                                                               ...
                                                                               height_above_ground4: 1,
                                                                               height_above_ground5: 2,
                                                                               potential_vorticity_surface: 2,
                                                                               sigma: 1,
                                                                               hybrid: 1,
                                                                               hybrid1: 2)
Coordinates: (12/25)
  * lat                                                                       (lat) float32 ...
  * lon                                                                       (lon) float32 ...
    reftime                                                                   datetime64[ns] ...
  * time                                                                      (time) datetime64[ns] ...
  * height_above_ground                                                       (height_above_ground) float32 ...
  * pressure_difference_layer                                                 (pressure_difference_layer) float32 ...
    ...                                                                        ...
  * height_above_ground4                                                      (height_above_ground4) float32 ...
  * height_above_ground5                                                      (height_above_ground5) float32 ...
  * potential_vorticity_surface                                               (potential_vorticity_surface) float32 ...
  * sigma                                                                     (sigma) float32 ...
  * hybrid                                                                    (hybrid) float32 ...
  * hybrid1                                                                   (hybrid1) float32 ...
Dimensions without coordinates: anon0, anon1, anon2, anon3, anon4, anon5, anon6
Data variables: (12/133)
    LatLon_Projection                                                         int32 ...
    pressure_difference_layer_bounds                                          (pressure_difference_layer, anon0) float32 ...
    height_above_ground_layer_bounds                                          (height_above_ground_layer, anon1) float32 ...
    height_above_ground_layer1_bounds                                         (height_above_ground_layer1, anon2) float32 ...
    pressure_difference_layer1_bounds                                         (pressure_difference_layer1, anon3) float32 ...
    pressure_difference_layer2_bounds                                         (pressure_difference_layer2, anon4) float32 ...
    ...                                                                        ...
    v-component_of_wind_isobaric                                              (time, isobaric, 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_ground2, 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:                                                             ...

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_ground2=20)
v = ds['v-component_of_wind_height_above_ground'].sel(time=timeObj,height_above_ground2=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] 2022-09-28T18:00:00
    time                  datetime64[ns] 2022-09-28T18:00:00
    height_above_ground2  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 = (27.0, -85.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.

previous

ERA-5 End-to-End Data and Graphics Preparation for the Science-on-a-Sphere

next

Week 12: 11/7