bokeh_choropleth_kvb_bike_density.py

import pandas as pd
import fiona
import numpy as np
from bokeh.io import show, output_file
from bokeh.models import ColumnDataSource, HoverTool, LogColorMapper
from bokeh.plotting import figure, save
from shapely.geometry import Polygon, Point, MultiPoint, MultiPolygon
from shapely.prepared import prep

""" from: https://www.offenedaten-koeln.de/dataset/59a8a033-5ac8-4240-ab06-608a7f542472/resource/a677cd63-d887-4f0f-95cf-25781641e576
 converted with http://www.gdal.org/ogr2ogr.html:
 ogr2ogr -lco ENCODING=UTF-8 -t_srs EPSG:4326 Stadtteil_WGS84.shp Stadtteil.shp """
SHAPEFILE="~/maps/offene/Stadtteil_WGS84.shp" 

def read_data(filename):
    colnames =["scrape_date","scrape_time","scrape_weekday","u_id","bike_id","lat","lon","bike_name"]
    with open(filename,"r") as f:
        data = pd.read_csv(filename, names=colnames,nrows=1000 )
    data = data[data.bike_name.str.contains("BIKE")].reset_index() 
    data.scrape_date = pd.to_datetime(data.scrape_date) # Convert to datetime object
    data.scrape_time = pd.to_datetime(data.scrape_time, format="%H-%M-%S")
    return data

def calc_points_per_poly(poly, points): # Returns number of points contained
    poly = prep(poly)
    return int(len(list(filter(poly.contains, points))))

# bike data scraped: https://data-dive.com/cologne-bike-rentals-getting-data
data = read_data("~/kvb/data/2017-03-01.csv") # Read bike location data
map_points = [Point(x,y) for x,y in zip(data.lon, data.lat)] # Convert Points to Shapely Points
all_points = MultiPoint(map_points) # all bike points

shp = fiona.open(SHAPEFILE)
# Extract features from shapefile
district_name = [ feat["properties"]["STT_NAME"] for feat in shp]
district_area = [ feat["properties"]["SHAPE_AREA"] for feat in shp]
district_x = [ [x[0] for x in feat["geometry"]["coordinates"][0]] for feat in shp]
district_y = [ [y[1] for y in feat["geometry"]["coordinates"][0]] for feat in shp]
district_xy = [ [ xy for xy in feat["geometry"]["coordinates"][0]] for feat in shp] 
district_poly = [ Polygon(xy) for xy in district_xy] # coords to Polygon

num_bikes = [ calc_points_per_poly(x, all_points) for x in district_poly]
bikes_per_area = [  x/y*10000 for x,y in zip(num_bikes, district_area) ]

# prepare plotting with bokeh
custom_colors = ['#f2f2f2', '#fee5d9', '#fcbba1', '#fc9272', '#fb6a4a', '#de2d26']
color_mapper = LogColorMapper(palette=custom_colors)    
source = ColumnDataSource(data=dict(
    x=district_x, y=district_y,
    name=district_name, rate=bikes_per_area,
))

TOOLS = "pan,wheel_zoom,reset,hover,save"
p = figure(
    title="KVB bike density per district, Mar. 2017", tools=TOOLS,
    x_axis_location=None, y_axis_location=None
)
p.grid.grid_line_color = None
p.patches('x', 'y', source=source,
          fill_color={'field': 'rate', 'transform': color_mapper},
          fill_alpha=0.8, line_color="black", line_width=0.3)

hover = p.select_one(HoverTool)
hover.point_policy = "follow_mouse"
hover.tooltips = [("District", "@name"),("Bikes per km²", "@rate"),("(Long, Lat)", "($x, $y)")]

output_file("kvb_interactive.html")
show(p)