This is an initial prototype for the cube-explorer project to integrate Iris, Cartopy and HoloViews to allow easily exploring Iris Cubes of N-dimensional gridded data.
For now we've called the library that interfaces between Iris/Cartopy and HoloViews holocube and we will import it as hc.
In [ ]:
import iris
import numpy as np
import holoviews as hv
import holocube as hc
from cartopy import crs
from cartopy import feature as cf
hv.notebook_extension()
The holocube package provides a library of Element types which make itvery easy to plot data on various geographic projections. Depending on the type of data the plotting code will automatically infer the correct transform and a default projection.
A simple example is the GeoFeature Element, which accepts various cartopy features. We can easily overlay these features by constructing an Overlay of GeoFeature elements.
In [ ]:
%%output size=400
feats = [cf.LAND, cf.OCEAN, cf.RIVERS, cf.LAKES, cf.BORDERS, cf.COASTLINE]
features = hv.Overlay([hc.GeoFeature(feature) for feature in feats])
features
Below is the full list of cartopy projections that can be displayed using matplotlib.
In [ ]:
projections = [crs.RotatedPole, crs.TransverseMercator, crs.Mercator, crs.LambertCylindrical,
crs.Geostationary, crs.AzimuthalEquidistant, crs.OSGB, crs.EuroPP, crs.Gnomonic,
crs.PlateCarree, crs.Mollweide, crs.OSNI, crs.Miller, crs.InterruptedGoodeHomolosine,
crs.LambertConformal, crs.SouthPolarStereo, crs.AlbersEqualArea, crs.Orthographic,
crs.NorthPolarStereo, crs.Robinson, crs.Stereographic]
We can test the different projections by creating a Layout of GeoFeature elements, each with a different projection:
In [ ]:
hv.Layout([hc.GeoFeature(cf.COASTLINE, group=p.__name__)(plot=dict(projection=p()))
for p in projections]).display('all')
To change the projection we can use the call method on HoloViews objects and set it as a plot option, this way we can easily compose plots with different projections. Here we compose a StockImage using the Mollweide projection with an Overlay of a StockImage and Coastlines Element set to a GeoStationary projection.
In [ ]:
%output size=250
In [ ]:
(features(plot=dict(projection=crs.Mollweide())) +
features.relabel(group='Geostationary Overlay')(plot=dict(projection=crs.Geostationary())))
This way we can also use different Element types such as this WMTS Element, which allows wrapping any webmap tilesource:
In [ ]:
%%output backend='matplotlib:nbagg' widgets='live' size=200
url = 'http://map1c.vis.earthdata.nasa.gov/wmts-geo/wmts.cgi'
layer = 'VIIRS_CityLights_2012'
hc.WMTS(url, layer=layer)(plot=dict(projection=crs.PlateCarree()))
In some cases the data does not define the coordinate system the data is in automatically, e.g. when using points. In this case we can supply a coordinate reference system directly. Here we display a GeoTiles object drawn from the MapQuest ordinance survey map of Great Britain, and overlay this tile source with points corresponding to the tube map locations. Since these coordinates are in Ordinance Survery GB coordinates we declare that explicitly via the crs parameter on the GeoPoints.
In [ ]:
%%output size=200
%%opts Overlay [apply_extents=True]
from cartopy.io.img_tiles import MapQuestOSM
from matplotlib.path import Path
def tube_locations():
"""
Returns an (n, 2) array of selected London Tube locations in Ordnance
Survey GB coordinates.
Source: http://www.doogal.co.uk/london_stations.php
"""
return np.array([[531738., 180890.], [532379., 179734.],
[531096., 181642.], [530234., 180492.],
[531688., 181150.], [530242., 180982.],
[531940., 179144.], [530406., 180380.],
[529012., 180283.], [530553., 181488.],
[531165., 179489.], [529987., 180812.],
[532347., 180962.], [529102., 181227.],
[529612., 180625.], [531566., 180025.],
[529629., 179503.], [532105., 181261.],
[530995., 180810.], [529774., 181354.],
[528941., 179131.], [531050., 179933.],
[530240., 179718.]])
theta = np.linspace(0, 2 * np.pi, 100)
circle_verts = np.vstack([np.sin(theta), np.cos(theta)]).T
concentric_circle = Path.make_compound_path(Path(circle_verts[::-1]),
Path(circle_verts * 0.6))
rectangle = Path([[-1.1, -0.2], [1, -0.2], [1, 0.3], [-1.1, 0.3]])
tiles = MapQuestOSM()
hc.GeoTiles(tiles)(plot=dict(projection=tiles.crs, zoom=14)) *\
hc.Points(tube_locations(), crs=crs.OSGB())(style=dict(color='r', s=100, marker=concentric_circle)) *\
hc.Points(tube_locations(), crs=crs.OSGB())(style=dict(color='b', s=100, marker=rectangle))
Now we will load some real data using from the Iris sample datasets:
In [ ]:
import numpy as np
import iris
import iris.coords
def realization_metadata(cube, field, fname):
if not cube.coords('realization'):
realization_number = fname[-6:-3]
realization_coord = iris.coords.AuxCoord(np.int32(realization_number), 'realization')
cube.add_aux_coord(realization_coord)
surface_temp = iris.load_cube(iris.sample_data_path('GloSea4', 'ensemble_???.pp'),
iris.Constraint('surface_temperature', realization=lambda value: True),
callback=realization_metadata)
The Cube Element defined above wraps the Iris Cubes, converting coordinates to HoloViews dimensions and tries to infer the correct order of dimensions:
In [ ]:
cube = hc.HoloCube(surface_temp)
cube
We'll come back to the Cube Element later for now we will slice this cube up manually. By taking slices along latitude and longitude we can slice the data up into 2D chunks and wrap them in GeoImage a subclass of Cube which can be visualized. We place these object into a HoloMap with the remaining dimensions time and realization as key dimensions.
In [ ]:
kdims = ['time', 'realization']
img_hmap = hv.HoloMap(kdims=kdims)
for cb in surface_temp.slices(['longitude', 'latitude']):
key = tuple(cb.coord(kd).points[0] for kd in kdims)
img_hmap[key] = hc.Image(cb)
The HoloMap can summarize the contained data:
In [ ]:
img_hmap.info
A convenient way of accessing a single Element in a HoloMap is the .last attribute. Now that we have a handle on it we can customize it a number of ways using the call method as above or using the options magic:
In [ ]:
%opts Image [colorbar=True projection=crs.Geostationary()] (cmap='viridis')
img_hmap.last * hc.GeoFeature(cf.COASTLINE)
</br>
In [ ]:
%output widgets='live' size=300
Slicing a Cube up in the way we saw before is often very useful but it's also a little bit of effort. To make this easier HoloViews interfaces usually implement a groupby method. Here we show how to achieve the same thing as above but using groupby instead. We may add another clearer interface eventually but groupby will provide the low level API for any such conversion interface.
If we group the cube by realization and time we are left with a bunch of 2D slices of latitude and longitude, by supplying a group_type we wrap each of the sliced 2D cubes in an Image.
In [ ]:
%%opts GeoImage [colorbar=True] (cmap='viridis')
(cube.groupby(['time', 'realization'], group_type=hc.Image) * hc.GeoFeature(cf.COASTLINE))
As you can see it has automatically converted the cube to an widget allowing you to explore this space. We can repeat the same groupby operation this time with a Contours Element as the group_type.
In [ ]:
%%opts Contours [colorbar=True] (cmap='viridis')
(cube.groupby(['time', 'realization'], group_type=hc.Contours) * hc.GeoFeature(cf.COASTLINE))
iris Cubes can also be used as the data source for regular Element types, however unlike their holocube counterparts they won't know anything about the geographic projections of the data.
As the most basic example we can use the conversion interface to display a widget for each time and realization, but this time we declare just the key dimension of the points object letting the interface infer the value and HoloMap dimensions:
In [ ]:
%%opts Points [color_index=2 size_index=None]
cube.to.points(['longitude', 'latitude'])
We can also collapse specific dimensions on the iris Cube first and then view the reduced Cube using regular HoloViews Element types. Here we collapse the longitude and latitude dimensions on the iris Cube by taking the weighted mean, wrap it in a HoloCube and then view the mean surface temperature for each realization and overlay the curves.
In [ ]:
%%opts Curve [aspect=2 xticks=4 ] (linestyle='--') NdOverlay [aspect=2 legend_position='right']
if cube.data.coord('latitude').bounds is None:
cube.data.coord('latitude').guess_bounds()
if cube.data.coord('longitude').bounds is None:
cube.data.coord('longitude').guess_bounds()
grid_weights = iris.analysis.cartography.area_weights(cube.data)
collapsed_cube = cube.data.collapsed(['longitude', 'latitude'], iris.analysis.MEAN, weights=grid_weights)
hc.HoloCube(collapsed_cube).to.curve(['time']).overlay()
Similarly we can collapse the forecast period, leaving just latitude and longitude coordinates and then view slices along each longitude:
In [ ]:
%%opts Curve [aspect=3 xticks=10]
collapsed_cube = cube.data.collapsed('forecast_period', iris.analysis.MEAN, weights=grid_weights)
hc.HoloCube(collapsed_cube).to.curve(['latitude'])
This notebook has outlined a basic API to explore cube datasets using HoloViews. While quite a bit works already there are a large number of issues to still be sorted out: