Geometries

In [1]:
import pandas as pd
import numpy as np
import holoviews as hv
import geoviews as gv
import geoviews.feature as gf
import cartopy
import cartopy.feature as cf
from cartopy import crs as ccrs
hv.notebook_extension()
%output dpi=120

Cartopy and shapely make working with geometries and shapes very simple, and GeoViews provides convenient wrappers for the various geometry types they provide. In addition to basic Path and Polygons types, which draw geometries from lists of arrays, GeoViews also provides the Feature and Shape types, which wrap cartopy Features and shapely geometries respectively.

Feature

The Feature Element provides a very convenient means of overlaying a set of basic geographic features on top of or behind a plot. The cartopy.feature module provides various ways of loading custom features, however geoviews provides a number of default features which we have imported as gf , amongst others this includes coastlines, country borders, and land masses. Here we demonstrate how we can plot these very easily, either in isolation or overlaid:

In [2]:
(gf.ocean + gf.land + gf.ocean * gf.land * gf.coastline * gf.borders).cols(3)
Out[2]:

These deafult features simply wrap around cartopy Features, therefore we can easily load a custom NaturalEarthFeature such as graticules at 30 degree intervals:

In [3]:
%%opts Feature.Lines [projection=ccrs.Robinson()] (facecolor='none' edgecolor='gray')
graticules = cf.NaturalEarthFeature(
    category='physical',
    name='graticules_30',
    scale='110m')
(gf.ocean * gf.land() * gv.Feature(graticules, group='Lines') * gf.borders * gf.coastline)
Out[3]:

Shape

The gv.Shape object wraps around any shapely geometry, allowing finer grained control over each polygon. We can, for example, access the geometries on the LAND feature and display them individually. Here we will get the geometry corresponding to the Australian continent and display it using shapely's inbuilt SVG repr (not yet a HoloViews plot, just a bare SVG displayed by Jupyter directly):

In [4]:
land_geoms = list(gf.land.data.geometries())
land_geoms[21]
Out[4]:

Instead of letting shapely render it as an SVG, we can now wrap it in the gv.Shape object and let matplotlib or bokeh render it, alone or with other GeoViews or HoloViews objects:

In [5]:
%%opts Points (color="black")
australia = gv.Shape(land_geoms[21], crs=ccrs.PlateCarree())

australia * gv.Points([(133.870,-23.700)]) * gv.Text(133.870,-21.5, 'Alice Springs', crs=ccrs.PlateCarree())
Out[5]:

We can also iterate over the geometries and wrap them all in an NdOverlay of gv.Shape Elements:

In [6]:
%%opts NdOverlay [aspect=2]
hv.NdOverlay({i: gv.Shape(s, crs=ccrs.PlateCarree()) for i, s in enumerate(land_geoms)})
Out[6]:

This makes it possible to create choropleth maps, where each part of the geometry is assigned a value that will be used to color it. However, constructing a choropleth by combining a bunch of shapes can be a lot of effort and is error prone. For that reason, the Shape Element provides convenience methods to load geometries from a shapefile. Here we load the boundaries of UK electoral districts directly from an existing shapefile:

In [7]:
shapefile = './assets/boundaries/boundaries.shp'
gv.Shape.from_shapefile(shapefile, crs=ccrs.PlateCarree())
Out[7]:

To combine these shapes with some actual data, we have to be able to merge them with a dataset. To do so we can inspect the records the cartopy shapereader loads:

In [8]:
shapes = cartopy.io.shapereader.Reader(shapefile)
list(shapes.records())[0]
Out[8]:
<Record: <shapely.geometry.multipolygon.MultiPolygon object at 0x124afac10>, {'code': 'E07000007'}, <fields>>

As we can see, the record contains a MultiPolygon together with a standard geographic code , which we can use to match up the geometries with a dataset. To continue we will require a dataset that is also indexed by these codes. For this purpose we load a dataset of the 2016 EU Referendum result in the UK:

In [9]:
referendum = pd.read_csv('./assets/referendum.csv')
referendum = hv.Dataset(referendum)
referendum.data.head()
Out[9]:
leaveVoteshare regionName turnout name code
0 4.100000 Gibraltar 83.500000 Gibraltar BS0005003
1 69.599998 North East 65.500000 Hartlepool E06000001
2 65.500000 North East 64.900002 Middlesbrough E06000002
3 66.199997 North East 70.199997 Redcar and Cleveland E06000003
4 61.700001 North East 71.000000 Stockton-on-Tees E06000004

The from_records function optionally also supports merging the records and dataset directly. To merge them, supply the name of the shared attribute on which the merge is based via the on argument. If the name of attribute in the records and the dimension in the dataset match exactly, you can simply supply it as a string, otherwise supply a dictionary mapping between the attribute and column name. In this case we want to color the choropleth by the 'leaveVoteshare' , which we define via the value argument. By default, the resulting NdOverlay of shapes will be indexed by an integer index. To draw the index from values in the dataset instead, you can request one or more indexes using the index argument. Finally we will declare the coordinate reference system in which this data is stored, which will in most cases be the simple Plate Carree projection. We can then view the choropleth, with each shape colored by the specified value (the percentage who voted to leave the EU):

In [10]:
%%opts Shape (cmap='viridis')
gv.Shape.from_records(shapes.records(), referendum, on='code', value='leaveVoteshare',
                     index=['name', 'regionName'], crs=ccrs.PlateCarree())
Out[10]:

The "Working with Bokeh" GeoViews notebook shows how to enable hover data that displays information about each of these shapes interactively.