Interactive maps with Bokeh

Our goal in this lesson is to learn few concepts how we can produce nice looking interactive maps using Geopandas and Bokeh.

Simple interactive point plot

First, we learn the basic logic of plotting in Bokeh by making a simple interactive plot with few points.

Import necessary functionalities from bokeh.

from bokeh.plotting import figure, output_file, show
from bokeh.io import output_notebook

output_notebook()

First we need to initialize our plot by calling the figure object.

# Initialize the plot (p) and give it a title
In [1]: p = figure(title="My first interactive plot!")
# Let's see what it is
p

Next we need to create lists of x and y coordinates that we want to plot.

# Create a list of x-coordinates
In [2]: x_coords = [0,1,2,3,4]

# Create a list of y-coordinates
In [3]: y_coords = [5,4,1,2,0]

Note

In Bokeh drawing points, lines or polygons are always done using list(s) of x and y coordinates.

Now we can plot those as points using a .circle() -object. Let’s give it a red color and size of 10.

# Plot the points
p.circle(x=x_coords, y=y_coords, size=10, color="red")

Finally, we can save our interactive plot into the disk with output_file -function that we imported in the beginning. All interactive plots are typically saved as html files which you can open in a web-browser.

# open the plot html, will open a new tab
show(p)
# Save the plot by passing the plot -object and output path
In [4]: outfp = "points.html"

In [5]: output_file(outfp)

Now open your interactive points.html plot by double-clicking it which should open it in a web browser.

This is how it should look like:

Bokeh Plot

It is interactive. You can drag the plot by clicking with left mouse and dragging. There are also specific buttons on the right side of the plot by default which you can select on and off:

../_images/Pan.png

Pan button allows you to switch the dragging possibility on and off (on by default).

../_images/BoxZoom.png

BoxZoom button allows you to zoom into an area that you define by left dragging with mouse an area of your interest.

../_images/Save.png

Save button allows you to save your interactive plot as a static low resolution .png file.

../_images/WheelZoom.png

WheelZoom button allows you to use mouse wheel to zoom in and out.

../_images/Reset.png

Reset button allows you to use reset the plot as it was in the beginning.

Creating interactive maps using Bokeh and Geopandas

Now we now khow how to make a really simple interactive point plot using Bokeh. What about creating such a map from a Shapefile of points? Of course we can do that, and we can use Geopandas for achieving that goal which is nice!

Creating an interactive Bokeh map from Shapefile(s) contains typically following steps:

  1. Read the Shapefile into GeoDataFrame

  2. Calculate the x and y coordinates of the geometries into separate columns

  3. Convert the GeoDataFrame into a Bokeh DataSource

  4. Plot the x and y coordinates as points, lines or polygons (which are in Bokeh words: circle, multi_line and patches)

Let’s practice these things and see how we can first create an interactive point map, then a map with lines, and finally a map with polygons where we also add those points and lines into our final map.

Point map

Let’s first make a map out schools address points in Tartumaa. That Shapefile is provided for you in the data folder that you downloaded.

Read the data using geopandas which is the first step.

In [6]: import geopandas as gpd

# File path
In [7]: points_fp = "source/_static/data/L6/schools_tartu.shp"

# Read the data
In [8]: points = gpd.read_file(points_fp)

Let’s see what we have.

In [9]: points.head()
Out[9]: 
      id                                     name                                            Aadress            X            y                        geometry
0  13376                   Aakre Lasteaed-Algkool  Valga maakond, Puka vald, Aakre küla, Mõisa te...  629862.0000  6441779.000  POINT (629862.000 6441779.000)
1  13290                       Alatskivi Lasteaed  Tartu maakond, Alatskivi vald, Alatskivi alevi...  682186.9651  6499629.999  POINT (682186.965 6499629.999)
2  13396                 Anna Haava nim Pala Kool    Jõgeva maakond, Pala vald, Pala küla, Koolimaja  675030.0000  6507478.000  POINT (675030.000 6507478.000)
3  13202                          Elva Gümnaasium              Tartu maakond, Elva linn, Tartu mnt 3  641493.0000  6456152.000  POINT (641493.000 6456152.000)
4  13298  Elva Huviala-ja Kultuurikeskus Sinilind               Tartu maakond, Elva linn, Kesk tn 30  641973.9625  6456196.951  POINT (641973.963 6456196.951)

Okey, so we have several columns plus the geometry column as attributes.

Now, as a second step, we need to calculate the x and y coordinates of those points. Unfortunately there is not a ready made function in geopandas to do that.

Thus, let’s create our own function called getPointCoords() which will return the x or y coordinate of a given geometry. It shall have two parameters: geom and coord_type where the first one should be a Shapely geometry object and coord_type should be either 'x' or 'y'.

In [10]: def getPointCoords(row, geom, coord_type):
   ....:     """Calculates coordinates ('x' or 'y') of a Point geometry"""
   ....:     if coord_type == 'x':
   ....:         return row[geom].x
   ....:     elif coord_type == 'y':
   ....:         return row[geom].y
   ....: 

Let’s then use our function using .apply().

# Calculate x coordinates
In [11]: points['x'] = points.apply(getPointCoords, geom='geometry', coord_type='x', axis=1)

# Calculate y coordinates
In [12]: points['y'] = points.apply(getPointCoords, geom='geometry', coord_type='y', axis=1)

# Let's see what we have now
In [13]: points.head()
Out[13]: 
      id                                     name                                            Aadress            X            y                        geometry            x
0  13376                   Aakre Lasteaed-Algkool  Valga maakond, Puka vald, Aakre küla, Mõisa te...  629862.0000  6441779.000  POINT (629862.000 6441779.000)  629862.0000
1  13290                       Alatskivi Lasteaed  Tartu maakond, Alatskivi vald, Alatskivi alevi...  682186.9651  6499629.999  POINT (682186.965 6499629.999)  682186.9651
2  13396                 Anna Haava nim Pala Kool    Jõgeva maakond, Pala vald, Pala küla, Koolimaja  675030.0000  6507478.000  POINT (675030.000 6507478.000)  675030.0000
3  13202                          Elva Gümnaasium              Tartu maakond, Elva linn, Tartu mnt 3  641493.0000  6456152.000  POINT (641493.000 6456152.000)  641493.0000
4  13298  Elva Huviala-ja Kultuurikeskus Sinilind               Tartu maakond, Elva linn, Kesk tn 30  641973.9625  6456196.951  POINT (641973.963 6456196.951)  641973.9625

Okey great! Now we have the x and y columns in our GeoDataFrame.

The third step, is to convert our DataFrame into a format that Bokeh can understand. Thus, we will convert our DataFrame into ColumnDataSource which is a Bokeh-specific way of storing the data.

Note

Bokeh ColumnDataSource does not understand Shapely geometry -objects. Thus, we need to remove the geometry -column before convert our DataFrame into a ColumnDataSouce.

Let’s make a copy of our points GeoDataFrame where we drop the geometry column.

# Make a copy and drop the geometry column
In [14]: p_df = points.drop('geometry', axis=1).copy()

# See head
In [15]: p_df.head(2)
Out[15]: 
      id                    name                                            Aadress            X            y            x
0  13376  Aakre Lasteaed-Algkool  Valga maakond, Puka vald, Aakre küla, Mõisa te...  629862.0000  6441779.000  629862.0000
1  13290      Alatskivi Lasteaed  Tartu maakond, Alatskivi vald, Alatskivi alevi...  682186.9651  6499629.999  682186.9651

Now we can convert that pandas DataFrame into a ColumnDataSource.

In [16]: from bokeh.models import ColumnDataSource

# Point DataSource
In [17]: psource = ColumnDataSource(p_df)
# What is it?
psource

Okey, so now we have a ColumnDataSource object that has our data stored in a way that Bokeh wants it.

Finally, we can make a Point map of those points in a fairly similar manner as in the first example. Now instead of passing the coordinate lists, we can pass the data as a source for the plot with column names containing those coordinates.

# Initialize our plot figure
In [18]: p = figure(title="A map of school location points from a Shapefile")

# Add the points to the map from our 'psource' ColumnDataSource -object
In [19]: p.circle('x', 'y', source=psource, color='red', size=10)
Out[19]: GlyphRenderer(id='1078', ...)
# open the plot html, will open a new tab
show(p)

And now the last thing is to save our map as html file into our computer.

# Output filepath
In [20]: outfp = "point_map.html"

# Save the map
In [21]: output_file(outfp)

Now you can open your point map in the browser. Your map should look like following:

Bokeh Plot

Adding interactivity to the map

In Bokeh there are specific set of plot tools that you can add to the plot. Actually all the buttons that you see on the right side of the plot are exactly such tools. It is e.g. possible to interactively show information about the plot objects to the user when placing mouse over an object as you can see from the example on top of this page. The tool that shows information from the plot objects is an inspector called HoverTool that annotate or otherwise report information about the plot, based on the current cursor position.

Let’s see now how this can be done.

First we need to import the HoverTool from bokeh.models that includes .

In [22]: from bokeh.models import HoverTool

Next, we need to initialize our tool.

In [23]: my_hover = HoverTool()

Then, we need to tell to the HoverTool that what information it should show to us. These are defined with tooltips like this:

In [24]: my_hover.tooltips = [('Name of the School', '@name')]

From the above we can see that tooltip should be defined with a list of tuple(s) where the first item is the name or label for the information that will be shown, and the second item is the column-name where that information should be read in your data.

The @ character in front of the column-name is important because it tells that the information should be taken from a column named as the text that comes after the character.

Lastly we need to add this new tool into our current plot.

In [25]: p.add_tools(my_hover)
# open the plot html, will open a new tab
show(p)

Great! Let’s save this enhanced version of our map as point_map_hover.html and see the result.

In [26]: outfp = "point_map_hover.html"

In [27]: output_file(outfp)
Bokeh Plot

As you can see now the plot shows information about the points and the content is the information derived from column address.

Hint

Of course, you can show information from multiple columns at the same time. This is achieved simply by adding more tooltip variables when defining the tooltips, such as:

my_hover2.tooltips = [('Label1', '@col1'), ('Label2', '@col2'), ('Label3', '@col3')]

Line map

Okey, now we have made a nice point map out of a Shapefile. Let’s see how we can make an interactive map out of a Shapefile that represents roads lines in Tartumaa. We follow the same steps than before, i.e. 1) read the data, 2) calculate x and y coordinates, 3) convert the DataFrame into a ColumnDataSource and 4) make the map and save it as html.

Read the data using geopandas which is the first step.

In [28]: import geopandas as gpd

# File path
In [29]: roads_fp = "source/_static/data/L6/roads.shp"

# Read the data
In [30]: roads = gpd.read_file(roads_fp)

Let’s see what we have.

In [31]: roads.head()
Out[31]: 
            TYYP                                           geometry
0  Kõrvalmaantee  LINESTRING (628395.967 6437374.941, 628230.700...
1  Kõrvalmaantee  LINESTRING (627016.000 6438859.770, 626942.160...
2  Kõrvalmaantee  LINESTRING (630462.518 6439009.709, 630496.794...
3    Põhimaantee  LINESTRING (630059.642 6438844.378, 630092.690...
4  Kõrvalmaantee  LINESTRING (630092.690 6439015.520, 630462.518...

Okey, so we have a road type plus the geometry column as attributes.

Second step is where calculate the x and y coordinates of the nodes of our lines.

Let’s create our own function called getLineCoords() in a similar manner as previously but now we need to modify it a bit so that we can get coordinates out of the Shapely LineString object.

In [32]: from shapely.geometry import LineString, MultiLineString

In [33]: def getLineCoords(row, geom, coord_type):
   ....:     if isinstance(row[geom], MultiLineString):
   ....:         empty_l = []
   ....:         return empty_l
   ....:     else:
   ....:         if coord_type == 'x':
   ....:             return list( row[geom].coords.xy[0] )
   ....:         elif coord_type == 'y':
   ....:             return list( row[geom].coords.xy[1] )
   ....: 

Note

Wondering about what happens here? Take a tour to our earlier materials about LineString attributes. By default Shapely returns the coordinates as a numpy array of the coordinates. Bokeh does not understand arrays, hence we need to convert the array into a list which is why we apply the Python built/in list() -function.

Let’s now apply our function in a similar manner as previously.

# Calculate x coordinates of the line
In [34]: roads['x'] = roads.apply(getLineCoords, geom='geometry', coord_type='x', axis=1)

# Calculate y coordinates of the line
In [35]: roads['y'] = roads.apply(getLineCoords, geom='geometry', coord_type='y', axis=1)

# Let's see what we have now
In [36]: roads.head()
Out[36]: 
            TYYP                                           geometry                                                  x                                                  y
0  Kõrvalmaantee  LINESTRING (628395.967 6437374.941, 628230.700...  [628395.967023925, 628230.7000000704, 628170.1...  [6437374.940890101, 6437599.020451155, 6437663...
1  Kõrvalmaantee  LINESTRING (627016.000 6438859.770, 626942.160...              [627016.0000000595, 626942.160000056]              [6438859.770451147, 6438884.15045115]
2  Kõrvalmaantee  LINESTRING (630462.518 6439009.709, 630496.794...             [630462.5180000596, 630496.7935589477]             [6439009.709451153, 6438938.423092961]
3    Põhimaantee  LINESTRING (630059.642 6438844.378, 630092.690...              [630059.6421779364, 630092.690000058]             [6438844.378158718, 6439015.520451145]
4  Kõrvalmaantee  LINESTRING (630092.690 6439015.520, 630462.518...              [630092.690000058, 630462.5180000596]             [6439015.520451145, 6439009.709451153]

Yep, now we have the x and y columns in our GeoDataFrame.

The third step. Convert the DataFrame (without geometry column) into a ColumnDataSource which, as you remember, is a Bokeh-specific way of storing the data.

# Make a copy and drop the geometry column
In [37]: m_df = roads.drop('geometry', axis=1).copy()

# Point DataSource
In [38]: msource = ColumnDataSource(m_df)

Finally, we can make a map of the roads line and save it in a similar manner as earlier but now instead of plotting circle we need to use a .multiline() -object. Let’s define the line_width to be 3.

# Initialize our plot figure
In [39]: p = figure(title="A map of Tartumaa Roads")

# Add the lines to the map from our 'msource' ColumnDataSource -object
In [40]: p.multi_line('x', 'y', source=msource, color='red', line_width=3)
Out[40]: GlyphRenderer(id='1121', ...)
# open the plot html, will open a new tab
show(p)

Great! Let’s save this and check the result.

In [41]: outfp = "roads_map.html"

# Save the map
In [42]: output_file(outfp)

Now you can open your point map in the browser and it should look like following:

Bokeh Plot

Polygon map with Points and Lines

It is of course possible to add different layers on top of each other. Let’s visualize a map showing Population per km2 in Tartumaa compared to Road network in Tartu, and locations of schools.

1st step: Import necessary modules and read the Shapefiles.

In [43]: from bokeh.plotting import figure, show, output_file, output_notebook

In [44]: from bokeh.models import ColumnDataSource, HoverTool, LogColorMapper

In [45]: import geopandas as gpd

In [46]: import geopandas as gpd

# Filepaths
In [47]: grid_fp = "source/_static/data/L6/population_square_km.shp"

In [48]: roads_fp = "source/_static/data/L6/roads.shp"

In [49]: schools_fp = "source/_static/data/L6/schools_tartu.shp"

# Read files
In [50]: grid = gpd.read_file(grid_fp)

In [51]: roads = gpd.read_file(roads_fp)

In [52]: schools = gpd.read_file(schools_fp)

As usual, we need to make sure that the coordinate reference system is the same in every one of the layers. Let’s use the CRS of the grid layer and apply it to our schools and roads line.

# Get the CRS of our grid
In [53]: CRS = grid.crs

In [54]: print(CRS)
epsg:3301

# Convert the geometries of roads line and schools into that one
In [55]: schools['geometry'] = schools['geometry'].to_crs(crs=CRS)

In [56]: roads['geometry'] = roads['geometry'].to_crs(crs=CRS)

Okey now, the geometries should have similar values:

In [57]: print(schools['geometry'].head(1))
0    POINT (629862.000 6441779.000)
Name: geometry, dtype: geometry

In [58]: print(roads['geometry'].head(1))
0    LINESTRING (628395.967 6437374.941, 628230.700...
Name: geometry, dtype: geometry

In [59]: print(grid['geometry'].head(1))
0    POLYGON ((625000.000 6433000.000, 624000.000 6...
Name: geometry, dtype: geometry

Indeed, they do. Let’s proceed and parse the x and y values of our grid. Let’s create own function for that as well.

In [60]: def getPolyCoords(row, geom, coord_type):
   ....:     """Returns the coordinates ('x' or 'y') of edges of a Polygon exterior"""
   ....:     exterior = row[geom].exterior
   ....:     if coord_type == 'x':
   ....:         return list( exterior.coords.xy[0] )
   ....:     elif coord_type == 'y':
   ....:         return list( exterior.coords.xy[1] )
   ....: 

2nd step: Let’s now apply the functions that we have created and parse the x and y coordinates for all of our datasets.

# Get the Polygon x and y coordinates
In [61]: grid['x'] = grid.apply(getPolyCoords, geom='geometry', coord_type='x', axis=1)

In [62]: grid['y'] = grid.apply(getPolyCoords, geom='geometry', coord_type='y', axis=1)

# Calculate x and y coordinates of the line
In [63]: roads['x'] = roads.apply(getLineCoords, geom='geometry', coord_type='x', axis=1)

In [64]: roads['y'] = roads.apply(getLineCoords, geom='geometry', coord_type='y', axis=1)

# Calculate x and y coordinates of the schools
In [65]: schools['x'] = schools.apply(getPointCoords, geom='geometry', coord_type='x', axis=1)

In [66]: schools['y'] = schools.apply(getPointCoords, geom='geometry', coord_type='y', axis=1)

Great, now we have x and y coordinates for all of our datasets. Let’s see how our grid coordinates look like.

# Show only head of x and y columns
In [67]: grid[['x', 'y']].head(2)
Out[67]: 
                                                   x                                                  y
0  [625000.0000000006, 624000.000000004, 624000.0...  [6432999.999999995, 6432999.999999995, 6434000...
1  [622999.9999999993, 622000.0000000027, 622000....  [6434000.0, 6434000.0, 6434999.999999996, 6434...

Let’s now classify the population data of our grid into 7 classes using a pysal classifier called Quantiles.

In [68]: import pysal.viz.mapclassify as mc

# Initialize the classifier and apply it
In [69]: classifier = mc.Quantiles.make(k=5)

In [70]: grid['pop_km2'] = grid[['Population']].apply(classifier)

# What do we have now?
In [71]: grid.head(2)
Out[71]: 
    STAMP_CRE  Population                                           geometry                                                  x                                                  y  pop_km2
0  1899-12-30           0  POLYGON ((625000.000 6433000.000, 624000.000 6...  [625000.0000000006, 624000.000000004, 624000.0...  [6432999.999999995, 6432999.999999995, 6434000...        0
1  1899-12-30           0  POLYGON ((623000.000 6434000.000, 622000.000 6...  [622999.9999999993, 622000.0000000027, 622000....  [6434000.0, 6434000.0, 6434999.999999996, 6434...        0

Okey, so we have many columns but the new one that we just got is the last one, i.e. pop_km2 that contains the classes that we reclassified based on the population per grid square.

3rd step: Let’s now convert our GeoDataFrames into Bokeh ColumnDataSources (without geometry columns)

# Make a copy, drop the geometry column and create ColumnDataSource
In [72]: m_df = roads.drop('geometry', axis=1).copy()

In [73]: msource = ColumnDataSource(m_df)

# Make a copy, drop the geometry column and create ColumnDataSource
In [74]: p_df = schools.drop('geometry', axis=1).copy()

In [75]: psource = ColumnDataSource(p_df)

# Make a copy, drop the geometry column and create ColumnDataSource
In [76]: g_df = grid.drop('geometry', axis=1).copy()

In [77]: gsource = ColumnDataSource(g_df)

Okey, now we are ready to roll and visualize our layers.

4th step: For visualizing the Polygons we need to define the color palette that we are going to use. There are many different ones available but we are now going to use a palette called RdYlBu and use eleven color-classes for the values (defined as RdYlBu11).

Let’s prepare our color_mapper.

# Let's first do some coloring magic that converts the color palet into map numbers (it's okey not to understand)
In [78]: from bokeh.palettes import RdYlBu11 as palette

In [79]: from bokeh.models import LogColorMapper

# Create the color mapper
In [80]: color_mapper = LogColorMapper(palette=palette)

Now we are ready to visualize our polygons and add the roads line and the schools on top of that. Polygons are visualized using patches objects in Bokeh.

# Initialize our figure
In [81]: p = figure(title="Population per km2 compared to Road network in Tartu")

# Plot grid
In [82]: p.patches('x', 'y', source=gsource,
   ....:          fill_color={'field': 'pop_km2', 'transform': color_mapper},
   ....:          fill_alpha=1.0, line_color="black", line_width=0.05)
   ....: 
Out[82]: GlyphRenderer(id='1165', ...)

# Add roads on top of the same figure
In [83]: p.multi_line('x', 'y', source=msource, color="red", line_width=2)
Out[83]: GlyphRenderer(id='1172', ...)

# Add schools on top (as yellow points)
In [84]: p.circle('x', 'y', size=3, source=psource, color="yellow")
Out[84]: GlyphRenderer(id='1179', ...)

# let's also add the hover over info tool
In [85]: tooltip = HoverTool()

In [86]: tooltip.tooltips = [('Name  of the school', '@name'),
   ....:                     ('Type of road', '@TYYP'),
   ....:                     ('Population density', '@Population')]
   ....: 

In [87]: p.add_tools(tooltip)
# open the plot html, will open a new tab
show(p)

Let’s save this and check the result.

In [88]: outfp = "roads_pop_km2_map.html"

# Save the map
In [89]: output_file(outfp)
Bokeh Plot

Cool, now we have an interactive map with three layers!

We can do one a few little more things. Let’s get a bit more experimental. As we are developing basically an interactive web map, let’s add a background base map.

We can use a so called tile source. In our example we use the CartoDB Positron free tile layer as our basemap background.

from bokeh.tile_providers import get_provider, Vendors

# instatiate the tile source provider
tile_provider = get_provider(Vendors.CARTODBPOSITRON)

# add the back ground basemap
p.add_tile(tile_provider)

One important aspect here is, that these web tile layers are often only available in WGS84 (EPSG:4326) or Web Mercator (EPSG:3857). We need to reproject our geometries if we want them to align with the background web map.

In addition, we create a legend scale bar that we can add to the plot. We use a different ColorMapper to visualise our population density grid.

from bokeh.palettes import RdPu9
from bokeh.models import LinearColorMapper, BasicTicker, ColorBar, Legend, LegendItem

# instantiate a new colormapper, provide min and max values for the scale
color_mapper_v = LinearColorMapper(palette=RdPu9,
                                low=grid_proj['pop_km2'].min(),
                                high=grid_proj['pop_km2'].max())

# build a legend scale for the grid population classes, using the LinearColorMapper from above
color_bar = ColorBar(color_mapper=color_mapper_v, ticker=BasicTicker(),
                    label_standoff=12, border_line_color=None, location=(0,0))

# and add this additional legend scale bar to the plot
p.add_layout(color_bar, 'right')

And in the end we build a custom legend. For that we need create separate LegendItems. The tricky bit is to link the legend item to the actual rendered (plotted) layer in the plot. We can access these separate renderers via the p.renderers[<pos>]-list. The order in this list is based on the order in which we place our layers on the map. The tile source background map is also rendered by a separate renderer. But it is our basemap, so we skip the zero’th element.

# now let's develop a custom legend for our 3 active layers
# the first rendered layer (p.renderers[0]) is the TileRenderer of our background map!
# here we link the rendered layer to a legend item
li1 = LegendItem(label='Pop Km2 Grid', renderers=[p.renderers[1]])
li2 = LegendItem(label='Roads', renderers=[p.renderers[2]])
li3 = LegendItem(label='Schools', renderers=[p.renderers[3]])

# then we put them together into an actual legend and place it on the plot
legend1 = Legend(items=[li1, li2, li3], location='top_right')
p.add_layout(legend1)

# we can even make this interactive and de-/activate layers
p.legend.click_policy="hide"

Altogether, that workflow would look like that:

from bokeh.tile_providers import get_provider, Vendors
from bokeh.palettes import RdPu9
from bokeh.models import LinearColorMapper, BasicTicker, ColorBar, Legend, LegendItem

# Convert the geometries
grid_proj = grid.to_crs(epsg=3857)
schools_proj = schools.to_crs(epsg=3857)
roads_proj = roads.to_crs(epsg=3857)

# Get the Polygon x and y coordinates
grid_proj['x'] = grid_proj.apply(getPolyCoords, geom='geometry', coord_type='x', axis=1)
grid_proj['y'] = grid_proj.apply(getPolyCoords, geom='geometry', coord_type='y', axis=1)

# Calculate x and y coordinates of the line
roads_proj['x'] = roads_proj.apply(getLineCoords, geom='geometry', coord_type='x', axis=1)
roads_proj['y'] = roads_proj.apply(getLineCoords, geom='geometry', coord_type='y', axis=1)

# Calculate x and y coordinates of the schools
schools_proj['x'] = schools_proj.apply(getPointCoords, geom='geometry', coord_type='x', axis=1)
schools_proj['y'] = schools_proj.apply(getPointCoords, geom='geometry', coord_type='y', axis=1)

# Make a copy, drop the geometry column and create ColumnDataSource
grid_proj = grid_proj.drop('geometry', axis=1).copy()
grid_proj_source = ColumnDataSource(grid_proj)

# Make a copy, drop the geometry column and create ColumnDataSource
roads_proj = roads_proj.drop('geometry', axis=1).copy()
roads_proj_source = ColumnDataSource(roads_proj)

# Make a copy, drop the geometry column and create ColumnDataSource
schools_proj = schools_proj.drop('geometry', axis=1).copy()
schools_proj_source = ColumnDataSource(schools_proj)

# Initialize our figure
p = figure(title="Population per km2 compared to Road network in Tartu",
        x_range=(2899000,3065000), y_range=(7980000,8140000),
        x_axis_type="mercator", y_axis_type="mercator")

# instatiate the tile source provider
tile_provider = get_provider(Vendors.CARTODBPOSITRON)

# add the back ground basemap
p.add_tile(tile_provider)

# instantiate a new colormapper
color_mapper_v = LinearColorMapper(palette=RdPu9,
                                low=grid_proj['pop_km2'].min(),
                                high=grid_proj['pop_km2'].max())

# Plot grid
p.patches('x', 'y', source=grid_proj_source,
        fill_color={'field': 'pop_km2', 'transform': color_mapper_v},
        fill_alpha=0.7, line_color="black", line_width=0.05)

# build a legend scale for the grid population classes
color_bar = ColorBar(color_mapper=color_mapper_v, ticker=BasicTicker(),
                    label_standoff=12, border_line_color=None, location=(0,0))

# and add this additional legend scale bar to the plot
p.add_layout(color_bar, 'right')

# Add roads on top of the same figure
p.multi_line('x', 'y', source=roads_proj_source, color="red", line_width=2)

# Add schools on top (as yellow points)
p.circle('x', 'y', size=3, source=schools_proj_source, color="yellow")

# now let's develop a custom legend for our 3 active layers
# the first rendered layer (p.renderers[0]) is the TileRenderer of our background map!
# here we link the rendered layer to a legend item
li1 = LegendItem(label='Pop Km2 Grid', renderers=[p.renderers[1]])
li2 = LegendItem(label='Roads', renderers=[p.renderers[2]])
li3 = LegendItem(label='Schools', renderers=[p.renderers[3]])

# then we put them together into an actual legend and place it on the plot
legend1 = Legend(items=[li1, li2, li3], location='top_right')
p.add_layout(legend1)

# we can even make this interactive and de-/activate layers
p.legend.click_policy="hide"

# let's also add the hover over info tool
tooltip = HoverTool()

tooltip.tooltips = [('Name  of the school', '@name'),
                    ('Type of road', '@TYYP')]
p.add_tools(tooltip)

show(p)
Bokeh Plot

Launch in the web/MyBinder:

https://mybinder.org/badge_logo.svg