We can easily display velocities as field of arrows in a leaflet map using mplleaflet, but what if we need to access more leaflet functionalities? Like adding a WMS layer and the layer control panel? We cannot use WMS with mplleaflet. Folium, on the other hand, can do that WMS layers and etc, but falls short when the job is to represent matplotlib elements...

So why not take the bet of two worlds? mplleaflet is built on top of mplexporter and works by exporting the figure elements as JSONs. If the element can be represent as a GeoJSON simplestyle-spec things are simple: just save the JSON and add it as a layer. When the element is more complex, and cannot be represent by the simplestyle-spec, mplexporter stores an SVG image that can be displayed as leaflet's divIcon.

All we need to do is to export the JSON with mplleaflet and implement leaflet divIcons on folium. In the example below we will try a quiver layer together with a WMS layer and a control panel.

First let's create the matplotlib figure!

In [3]:

import numpy as np
from netCDF4 import Dataset
from oceans import wrap_lon180


nc = Dataset('./data/mdt_cnes_cls2009_global_v1.1.nc')

u = nc['Grid_0002'][:]
v = nc['Grid_0003'][:]
lon = nc['NbLongitudes'][:]
lat = nc['NbLatitudes'][:]
err = nc['Grid_0004'][:]

lon = wrap_lon180(lon)
lon, lat = np.meshgrid(lon, lat)

In [4]:

mask_x = np.logical_and(lon > -40, lon < -34)
mask_y = np.logical_and(lat > -16, lat < -12)
mask = np.logical_and(mask_x, mask_y)

In [5]:

x = lon[mask]
y = lat[mask]
U = u.T[mask].filled(fill_value=0)
V = v.T[mask].filled(fill_value=0)

In [6]:

import matplotlib.pyplot as plt

fig, ax = plt.subplots()

kw = dict(color='black', alpha=0.8, scale=1)
q = ax.quiver(x, y, U, V, **kw)

We are plotting Mean Geostrophic Velocity derived from AVISO altimetry.

Now we need to export the figure to a GeoJSON folium can consume. The magic happens in the line gj = mplleaflet.fig_to_geojson(fig=fig). The gj is a "GeoJSON-like" Python dictionary where each Point has a html property containing its own SVG icon. In this case a velocity vector.

In [7]:

import mplleaflet

gj = mplleaflet.fig_to_geojson(fig=fig)

Now we can put the map together using folium.features.DivIcon.

In [8]:

import folium

feature_group0 = folium.FeatureGroup(name='quiver')

mapa = folium.Map(location=[y.mean(), x.mean()], tiles="Cartodb Positron",
                  zoom_start=7)

for feature in gj['features']:
    if feature['geometry']['type'] == 'Point':
        lon, lat = feature['geometry']['coordinates']
        div = feature['properties']['html']

        icon_anchor = (feature['properties']['anchor_x'],
                       feature['properties']['anchor_y'])

        icon = folium.features.DivIcon(html=div,
                                       icon_anchor=icon_anchor)
        marker = folium.Marker([lat, lon], icon=icon)
        feature_group0.add_children(marker)
    else:
        msg = "Unexpected geometry {}".format
        raise ValueError(msg(feature['geometry']))

url = 'http://gmrt.marine-geo.org/cgi-bin/mapserv?map=/public/mgg/web/gmrt.marine-geo.org/htdocs/services/map/wms_merc.map'
wms = folium.features.WmsTileLayer(url,
                                   name='GMRT',
                                   format='image/png',
                                   layers='topo',
                                   transparent=True)
feature_group1 = folium.FeatureGroup(name='Topo')
feature_group1.add_children(wms)

In [9]:

mapa.add_children(feature_group0)
mapa.add_children(feature_group1)
mapa.add_children(folium.map.LayerControl())
mapa

Out[9]:

Hover the mouse over the layer control (top right) to enable/disable layers. The control layer is useful to explore the data in the light of various useful layers, like bathymetry, SST, SSH, etc.

Note that the layers in this case do not fit perfectly. That could be a projection problem! We are naively assuming that everything fits nicely in a Web Mercator projection.

Be sure the check the crowding funding campaign on the side bar of the site!

Update: If you cannot see the side bar it is because your browser cannot display mixed HTTP/HTTPS. In that case click here.

In [10]:

HTML(html)

Out[10]:

This post was written as an IPython notebook. It is available for download or as a static html.

python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.

python4oceanographers

Turning ripples into waves

Displaying a field of arrows with folium

Comments