In [1]:
from bokeh.io import output_notebook, show
output_notebook()
Loading BokehJS ...

Plotting WRF output with Bokeh

I recently started a Bokeh tutorial, and I'm starting to think it's a great package. With its focus on the browser and its integration with javascript to make fast and interactive visualizations, it might really become a useful alternative to matplotlib. However, it is only just starting to implement tools for visualization of graphic data, and the documentation about this aspect is still a bit poor. So I decided to try and plot some WRF output data using Bokeh, perhaps it can be useful for others as well.

Note that there is also Holoview and Geoview (which is also very new and poorly documented), but I wanted to keep it simple and stayed close to the Bokeh basics.

Setting up a map

The first step is to set up a map. Basically, that boils down to opening an empty figure and putting on top of it a set of tiles that constitute the map. These tiles are taken from the web, for example from google earth or openstreetmap. These are referred to as Web Map Tile Service (WMTS). All these web services use the web mercator coordinate system, which roughly ranges from -20M to 20M meters (in fact I think they're not exactly meters, but let's ignore that for now) in both the west-east and the north-south direction.

A number of WMTS links are readily available in Bokeh, other will need to be loaded using a specific URL. To produce an empty map, we can use:

In [2]:
from bokeh.plotting import figure
from bokeh.tile_providers import STAMEN_TERRAIN, STAMEN_TONER, STAMEN_TONER_LABELS, STAMEN_TONER_BACKGROUND

# Create a figure with a really large extent
bound = 20000000 # meters
fig = figure(tools='pan, wheel_zoom', x_range=(-bound, bound), y_range=(-bound, bound))

# Hide the axes and add the map tile to the plot
fig.axis.visible = False
fig.add_tile(STAMEN_TERRAIN)

show(fig)

or, alternatively, supply a URL to WMTSTileSource

In [3]:
from bokeh.tile_providers import WMTSTileSource

wikimap = WMTSTileSource(url='https://maps.wikimedia.org/osm-intl/{Z}/{X}/{Y}@2x.png')

# Create a figure with a really large extent
bound = 20000000 # meters
fig = figure(tools='pan, wheel_zoom', x_range=(-bound, bound), y_range=(-bound, bound))

# Hide the axis and add the map tile
fig.axis.visible = False
fig.add_tile(wikimap)

show(fig)

Adding data to the map

To add data to the map is as simple as plotting, for example, some circles to the right grid coordinates. Because the axis extent is really large, we also need really large circles and put the far enough apart to actually be able to see them. Let's put some circles on the equator:

In [4]:
# Create a figure with a really large extent
bound = 20000000 # meters
fig = figure(tools='pan, wheel_zoom', x_range=(-bound, bound), y_range=(-bound, bound))

# Hide the axis and add the map tile
fig.axis.visible = False
fig.add_tile(wikimap)

# Add some circles to the map
fig.circle([0,1000000,2000000,3000000],[0,0,0,0],radius=100000)

show(fig)

Transforming your data to the web mercator coordinate system

In order to visualize georeferenced data, it is necessary to know the exact coordinates in the appropriate projection, i.e. the web mercator projection. This can be done manually, or by using Cartopy. In both cases I'll asume the original data is in lat/lon coordinates.

For the manual approach, I took the transformation function from this link. I then use this function to set the extent of the map and plot a center point on it.

In [5]:
import numpy as np

# http://philippjfr.com/work/work-in-progress/bokeh-maptiles/
def toWebMercator(xLon, yLat):
    # Check if coordinate out of range for Latitude/Longitude
    if (np.abs(xLon) > 180) and (np.abs(yLat) > 90):
        return
 
    semimajorAxis = 6378137.0  # WGS84 spheriod semimajor axis
    east = xLon * 0.017453292519943295
    north = yLat * 0.017453292519943295
 
    northing = 3189068.5 * np.log((1.0 + np.sin(north)) / (1.0 - np.sin(north)))
    easting = semimajorAxis * east
 
    return easting, northing

l, b = toWebMercator(-5, 45)
r, t = toWebMercator(5, 60)
cx,cy = toWebMercator(0,52.5)

extents = (l, b, r, t)

fig = figure(tools='pan, wheel_zoom', x_range=(l,r), y_range=(b,t))
fig.axis.visible = False
fig.add_tile(wikimap)

# Add center point
fig.circle(cx,cy,radius=100000,color='firebrick')

show(fig)

With cartopy, it is very easy to transform between coordinate systems. However, it took me some time to find out that there is actually a projection available for the web mercator projection under the name GOOGLE_MERCATOR. In Cartopy, the lat/lon coordinate system is called PlateCarree, so to transform between the two we use the following:

In [6]:
from cartopy import crs as ccrs
from bokeh.models import Range1d # an alternative way to define the axis extent, not sure what the difference is.

x1,y1 = -5,45
x2,y2 = 5,60

l,b = ccrs.GOOGLE_MERCATOR.transform_point(x1,y1,ccrs.PlateCarree())
r,t = ccrs.GOOGLE_MERCATOR.transform_point(x2,y2,ccrs.PlateCarree())
cx,cy = ccrs.GOOGLE_MERCATOR.transform_point(0,52.5,ccrs.PlateCarree())

fig = figure(tools='pan, wheel_zoom', x_range=Range1d(l,r), y_range=Range1d(b,t))
fig.axis.visible = False
fig.add_tile(wikimap)

# Add center point
fig.circle(cx,cy,radius=100000,color='firebrick')

show(fig)

I figured that it should also be possible to load the web-mercator projection into cartopy by using ccrs.epsg(3857), but it seems this redirects to the default mercator projection, so this does not work as I expected.

Loading some data from my WRF output

Here I use the xarray package to quickly load some WRF data from my output. I will need latitude, longitude and some toy variable, for which I simply use T2. Note that I don't load all the points (600x600), but use every 10th point instead to save some memory.

In [11]:
import xarray as xr

wrfout = xr.open_dataset('wrfout_d01')
lon = wrfout.XLONG.isel(Time=1).values[::10,::10]
lat = wrfout.XLAT.isel(Time=1).values[::10,::10]
T2 = wrfout.T2.isel(Time=1).values[::10,::10]

Plotting the data on a map

Now, in matplotlib we have the possibility to plot 2-dimensional geographic data as contour lines, pcolormeshes, etc. This functionality is not (yet?) available for Bokeh, but a suitable alternative is to plot circles and color the according to their value. To this end, we convert the T2 array to an array with color codes (using matplotlib here for simplicity). Then, we use cartopy's convert_points function to get the right coordinates in web-mercator projection. Then, all we need to do is plot all the points on the map and we're done!

In [15]:
import matplotlib as mpl
from matplotlib.cm import viridis
# Transform T2 to color codes
colors = [
    "#%02x%02x%02x" % (int(r), int(g), int(b)) for r, g, b, _ in 255*viridis(mpl.colors.Normalize()(T2.flatten()))
]

# Transform lat and lon to coordinates in web-mercator projection
points = ccrs.GOOGLE_MERCATOR.transform_points(ccrs.PlateCarree(),lon,lat)
x = points[:,:,0].flatten()
y = points[:,:,1].flatten()

# Make a figure, add the map
fig = figure(tools='pan, wheel_zoom')

# Plot the points (this will automatically set a range for the figure)
fig.circle(x,y,color=colors)

# Add the map
fig.axis.visible = False
fig.add_tile(wikimap)

show(fig)