Plotting contours with ipyleaflet

Question

Trying to get away from Basemap for plotting contours on maps, and in particular doing this from a Jupyter notebook, I came across this discussion related to folium: https://github.com/python-visualization/folium/issues/958. Hey! I thought - this is great, let me try it out with ipyleaflet.

The following code gets me close to a working solution, but I see some obscure contours as demonstrated in the image below.

Note: code for reading data is skipped

import matplotlib.pyplot as plt
from ipyleaflet import Map, basemaps, Polygon

# generate matplotlib contours
cs = plt.contourf(lat, lon, val.T)  # note transposition to obtain correct coordinate order
# set-up map
zoom = 4
center = [50., 0.]
map = Map(basemap=basemaps.CartoDB.Positron, center=center, zoom=zoom)
# add contours as polygons
# get the pathlines of the contours with cs.allsegs
for clev in cs.allsegs:
    polygons = Polygon(
                locations=[p.tolist() for p in clev],
                color="green",
                weight=1,
                opacity=0.4,
                fill_color="green",
                fill_opacity=0.3
    )
    map.add_layer(polygons);
map

For the image below, I loaded data from ECMWF's ERA-5 reanalysis and plotted horizontal wind speeds in February 2020. Only one contour was selected to demonstrate the problem. The left panel in the figure below shows the plt.contourf(lon, lat, val) result, which is correct. The panel on the right shows the contour drawn with leaflet. While much of the contour is displayed correctly (please focus only on the lime green area, i.e. the second-highest contour level), there seems to be some issue with ordering of line segments. Does anyone know how this can be fixed?

Answer 1

With help of the function split_contours below, the code now works as expected. The documentation of Path indicates that the contour attribute allkinds may be None. This is handled by the code below. What is not included, are path segment codes 3 and 4, which correspond to Bezier curves.

cs = plt.contourf(lat, lon, val.T)
plt.close()

# get the pathlines of the contours:
# print(cs.allsegs[0][0][0:12])  # returns list of polygons per contour level; each polygon is a list of [x,y]tuples
# i.e. [level][polygon][x,y]
# print(len(cs.allsegs), len(cs.allsegs[0]))

# other useful information from the contour object:
# print(cs.get_array())  # get contour levels

# plot map
from ipyleaflet import Map, basemaps, Polygon

zoom = 4
center = [50., 0.]
# map = Map(basemap=basemaps.NASAGIBS.ViirsTrueColorCR, center=center, zoom=zoom) # loads current satellite image
# map = Map(basemap=basemaps.NASAGIBS.ViirsEarthAtNight2012, center=center, zoom=zoom)
map = Map(basemap=basemaps.CartoDB.Positron, center=center, zoom=zoom)

def split_contours(segs, kinds=None):
    """takes a list of polygons and vertex kinds and separates disconnected vertices into separate lists.
    The input arrays can be derived from the allsegs and allkinds atributes of the result of a matplotlib
    contour or contourf call. They correspond to the contours of one contour level.
    
    Example:
    cs = plt.contourf(x, y, z)
    allsegs = cs.allsegs
    allkinds = cs.allkinds
    for i, segs in enumerate(allsegs):
        kinds = None if allkinds is None else allkinds[i]
        new_segs = split_contours(segs, kinds)
        # do something with new_segs
        
    More information:
    https://matplotlib.org/3.3.3/_modules/matplotlib/contour.html#ClabelText
    https://matplotlib.org/3.1.0/api/path_api.html#matplotlib.path.Path
    """
    if kinds is None:
        return segs    # nothing to be done
    # search for kind=79 as this marks the end of one polygon segment
    # Notes: 
    # 1. we ignore the different polygon styles of matplotlib Path here and only
    # look for polygon segments.
    # 2. the Path documentation recommends to use iter_segments instead of direct
    # access to vertices and node types. However, since the ipyleaflet Polygon expects
    # a complete polygon and not individual segments, this cannot be used here
    # (it may be helpful to clean polygons before passing them into ipyleaflet's Polygon,
    # but so far I don't see a necessity to do so)
    new_segs = []
    for i, seg in enumerate(segs):
        segkinds = kinds[i]
        boundaries = [0] + list(np.nonzero(segkinds == 79)[0])
        for b in range(len(boundaries)-1):
            new_segs.append(seg[boundaries[b]+(1 if b>0 else 0):boundaries[b+1]])
    return new_segs
        
# add contours as polygons
# hardwired colors for now: these correspons to the 8-level default of matplotlib with an added orange color
colors = ["#48186a", "#424086", "#33638d", "#26828e", "#1fa088", "#3fbc73", "#84d44b", "#d8e219", "#fcae1e"]
allsegs = cs.allsegs
allkinds = cs.allkinds

for clev in range(len(cs.allsegs)):
    kinds = None if allkinds is None else allkinds[clev]
    segs = split_contours(allsegs[clev], kinds)
    polygons = Polygon(
                    locations=[p.tolist() for p in segs],
                    # locations=segs[14].tolist(),
                    color=colors[min(clev, 4)],
                    weight=1,
                    opacity=0.8,
                    fill_color=colors[clev],
                    fill_opacity=0.5
    )
    map.add_layer(polygons);
map

The result looks like this: Now one can of course play with leaflet options and use different background maps etc.