bokehで地図を表示する際のあれこれ

bokeh学習の記録。

2022/6/18: Googleマップの表示、メッシュデータを重ねて表示を追加。
2023/1/15: Bokehバージョン3でのタイルサービスの使い方を追加。

注):Bokehバージョン3ではplot_width, plot_heightが廃止されwidth, heightのみになった。

最初に

使うもの

• bokeh
• geopandas

参考サイト

• bokeh
• GeoPandas
• Shapely

• ipywidgetsとBokeh使ってインタラクティブな可視化をする

ページ内のコード

• JupyterLabで作成

バージョン情報

• python 3.9.6
• bokeh 2.4.3
• geopandas 0.9.0
• jupyterlab 3.1.4

---

bokehについて

ブラウザやjupyter notebook上でデータのインタラクティブな可視化ができるライブラリ。
pythonだけでグリグリ動かせるグラフをつくり、データを差し替えたりできる。

GeoPandasについて

pandasを拡張して地理空間データを扱えるようにしたライブラリ

---

やること

• geopandasに入っている世界地図のデータセットを使い、bokehで地図を表示する。
• bokehの機能を使ってマップタイルを表示する。

import geopandas as gpd

from bokeh.plotting import figure, output_notebook, show 
from bokeh.models import GeoJSONDataSource
from bokeh.layouts import column, row

output_notebook()  # bokehのプロットをjupyter上で表示するようにする。

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.head()

	pop_est	continent	name	iso_a3	gdp_md_est	geometry
0	920938	Oceania	Fiji	FJI	8374.0	MULTIPOLYGON (((180.00000 -16.06713, 180.00000...
1	53950935	Africa	Tanzania	TZA	150600.0	POLYGON ((33.90371 -0.95000, 34.07262 -1.05982...
2	603253	Africa	W. Sahara	ESH	906.5	POLYGON ((-8.66559 27.65643, -8.66512 27.58948...
3	35623680	North America	Canada	CAN	1674000.0	MULTIPOLYGON (((-122.84000 49.00000, -122.9742...
4	326625791	North America	United States of America	USA	18560000.0	MULTIPOLYGON (((-122.84000 49.00000, -120.0000...

geometry列が地図データ。中身はshapelyのPolygon, MultiPolygon。その他は名前、人口、GDPなど。

地図の表示や設定

とにかく地図を表示する

patchesメソッドとGeoJSONDataSourceを使う。地図の拡大縮小はツールのWheelZoomをアクティブにして行う。p.patches('xs', 'ys')はGeoJSONDataSourceを使う時の決まり文句。

source = GeoJSONDataSource(geojson=world.to_json())

p = figure()
p.patches('xs', 'ys', source=source)

show(p)

見栄えを整える

形が崩れてしまっているので、ぱっと見違和感のないようにする。match_aspcet=Trueでx軸とy軸の1が同じ長さになる。正方形がちゃんと正方形に表示されるようになる。aspect_ratioは縦と横の比率。縦に対しての横の長さ。

source = GeoJSONDataSource(geojson=world.to_json())

p = figure(
        match_aspect=True,
        aspect_ratio=2,
        plot_width=800  # 横幅指定
)

p.patches('xs', 'ys', source=source)

show(p)

さらに調整

形は整ったので、色を変えたり、サイズを変えたりしてみる。sizing_modeでプロットのサイズをウィンドウの縦横に合わせて可変にできる。おおまかに言うと、'fixed'が固定、'stretch_'は形を無視して可変、'scale_'は形を保って可変。詳しくは公式参照。x_range.range_paddingで描画範囲内のパディングを設定できる。デフォルトは0.1。x_rangeとy_rangeのパディングが違うと、形が崩れる。.axis, .gridでxとyの両方の軸やグリッドを設定できる。.xaxis,.ygridのように個別の設定もできる。

source = GeoJSONDataSource(geojson=world.to_json())

p = figure(
        match_aspect=True,
        aspect_ratio=2,
        sizing_mode='scale_both',
        min_height=300,  # 縦幅の最低値指定
        background_fill_color='#222222',  # 背景色
        border_fill_color='#fafafa',  # 縁の色
        active_scroll='wheel_zoom'  # wheel_zoomを最初からアクティブに
)
p.x_range.range_padding = 0.02 # paddingの変更。0でパディングなし
p.y_range.range_padding = 0.02

p.grid.grid_line_color = 'grey'
p.xaxis.minor_tick_line_color = None
p.yaxis.minor_tick_line_color = None  # 色設定 Noneで消える

p.patches('xs', 'ys',
          source=source,
          fill_color='orange',  # 色
          fill_alpha=0.7,  # 透明度
          line_color='black',  # 輪郭線の色
          line_alpha=0.8,
          line_width=0.4  # 輪郭線の幅
          )

show(p)

点を重ねて表示する

gdf.centroidで各ポリゴンの重心を取得できる。マルチポリゴンだと点がポリゴンの外になることがある。representative_point()メソッドだと点はポリゴン内になる。centroidとは位置が異なる。centroidを使うと警告が出るが、点の位置に不満がなければ無視しても大丈夫。警告の中身に関してはstackoverflowの質問とstackexchangeの質問を自動翻訳するとわかったような気になれた。
circle('x', 'y')で円を重ねて表示する。GeoJSONDataSourceに渡したgeojson形式のデータはbokehで使える形に変換される。Pointは'x', 'y'の列、Polygonは'xs', 'ys'の列になるため同じ名前を使った列があると上書きされる。

points = world.representative_point()

print(points.head())

source = GeoJSONDataSource(geojson=world.to_json())
points_source = GeoJSONDataSource(geojson=points.to_json())

p = figure(
        match_aspect=True,
        aspect_ratio=2.2,
        plot_width=800
)

p.patches('xs', 'ys',
          source=source,
          fill_color='lightblue',
          line_width=0.4
         )

circle = p.circle('x', 'y',
                  source=points_source,
                  size=6,
                  fill_color='orange',
                  alpha=0.8
                  )
show(p)

0    POINT (177.97595 -17.93762)
1      POINT (34.14207 -6.20783)
2     POINT (-12.57202 24.23056)
3    POINT (-110.24381 56.70192)
4     POINT (-99.31483 37.23675)
dtype: geometry

表示範囲を指定する

x_range,y_rangeで描画範囲を設定できる。x_range=(100, 160)のようにタプルで渡せる。
描画範囲が固定されるのでmatch_aspectは効かなくなる。

形を崩さないようにするには自分で縦横の長さを設定する必要がある。
中央のフレーム部分の縦横を指定するframe_heightとframe_widthを使うのが簡単だがsizing_modeとの相性が悪い。sizing_modeを'fixed'以外にすると思った通りの形にならない可能性がある。plot_width, plot_heightやaspect_ratioだとボーダー部分も考慮する必要がある。

2022/6/18 追記:わざわざDataRange1dを使う意味がないと気付いたのでコードを修正。

source = GeoJSONDataSource(geojson=world.to_json())

p = figure(
        frame_width=300,
        frame_height=250,
        x_range=(100, 160),  # 横の範囲を指定
        y_range=(20, 70),  # 縦の範囲を指定
)
p.patches('xs', 'ys', source=source)

p2 = figure(
        plot_width=400,
        aspect_ratio= 1.2852,
        x_range=(100, 160),
        y_range=(20, 70),
)
p2.patches('xs', 'ys', source=source)

show(row(p, p2))

ソースのデータをもとに色をつける

国の名前ごとに塗り分ける

ソースに色の値を持つ列を加えるかfactor_cmap()を使う。第一引数から適用するソースの列名、色のリスト、色を適用する値のリストとなる。第三引数が['Japan', 'China']だけだった場合、日本と中国だけに指定した色がつき、それ以外はnan_colorの値になる。ここでは全ての国を塗り分けるので、そのままworld['name'](値はすべて固有値)を渡している。

from bokeh.transform import factor_cmap
from bokeh.palettes import magma, turbo

# ソースに色の列を追加
world['color'] = magma(177)

source = GeoJSONDataSource(geojson=world.to_json())

p = figure(
        match_aspect=True,
        aspect_ratio=2.2,
        plot_width=800,
        title='ソースの列名を指定して',
        x_axis_type=None,  # 軸とグリッドの表示なし
        y_axis_type=None,
)

p.patches('xs', 'ys',
          source=source,
          fill_color='color',  # sourceの列名を指定
          alpha=0.6)

p2 = figure(
        match_aspect=True,
        aspect_ratio=2.2,
        plot_width=800,
        title='factor_cmapを使って',
        x_axis_type=None,
        y_axis_type=None,
)    

mapper = factor_cmap('name', turbo(177), world['name'])

p2.patches('xs', 'ys',
           source=source,
           fill_color=mapper,
           alpha=0.6)

show(column(p, p2, background='white'))

値に応じて塗り分ける

linear_cmap()またはlog_cmap()を使う。第一引数から適用するソースの列名、色のリスト。lowとhighで範囲を指定する。linear_cmap()は指定した範囲を色数で等分。log_cmap()は指定した範囲を対数スケールで分割。数学?の自分はよく理解していないが、例えば、low=1, high=2**10とした場合、色の数を指数と同じ10にすれば、ちょうど2の何乗で色がわかれるようになる。下の例ではlow=10**3,high=10**10としているので、色の数を7にしている。

ColorBarでカラーバーを作成する。color_mapper引数にはColorMapperオブジェクトを渡す。log_cmap()の戻り値は辞書で、'field'と'transform'のkeyを持っている。'field'の値は適用する列名、'transform'はLogColorMapperオブジェクトになる。LogColorMapperを渡した時のカラーバーのデフォルトの目盛は10の累乗。例えば、2の累乗の目盛にしたい場合はticker=bokeh.models.LogTicker(base=2)のようにTickerを渡す。

p.add_layout(color_bar, 'right')でカラーバーを右側に追加する。他に'above', 'below', 'left','center'を指定できる。'center'は描画範囲内、それ以外は縁の位置。

from bokeh.transform import log_cmap
from bokeh.palettes import Cividis7
from bokeh.models import ColorBar

source = GeoJSONDataSource(geojson=world.to_json())

p = figure(
        title='人口による塗分け',
        match_aspect=True,
        aspect_ratio=2.2,
        plot_width=800,
)

mapper = log_cmap('pop_est', Cividis7[::-1], low=10**3, high=10**10)

patches = p.patches('xs', 'ys',
                    source=source,
                    fill_color=mapper,
                    line_color='#cdcdcd',
                    line_width=0.2
                    )

color_bar = ColorBar(color_mapper=mapper['transform'])

p.add_layout(color_bar, 'right')

show(p)

ホバーツールを使う

HoverToolクラスを使う。tooltipsがホバーツールの中身になる。渡すのはlistかstring。listの値は(str, str)のタプル。左が見出し、右が適用する列名か特殊な変数。列名は頭に'@'をつける。特殊な変数は'$'をつける。変数の詳細は公式のHoverToolを参照。列名にスペースが含まれたり、日本語だったりする場合は'{}'で囲う。列名のあとの'{}'はformatの設定。{0,0}で数字に'，'を挿入。詳細は公式を参照。
stringはhtml形式で書ける。HoverTool()のインスタンス作成時にtooltips以外の引数を渡さないなら、figure()の引数で直接tooltips=と渡せる。

ツールについて

右側の縦に並んでいるアイコンがツール。デフォルトはPanTool, BoxZoomTool, WheelZoomTool, SaveTool, ResetTool, HelpTool。このツールを通じて、地図をグリグリ動かしたり、拡大縮小したりしている。例えば、TapToolを追加してアクティブにすると、描画範囲内をクリック、タップしてポリゴンを選択できるようになる。figure()の引数としてtools='pan, wheel_zoom, tap'のように文字列で渡せば、任意のツールを設定できる。p.tools属性でツール類の入ったリストを参照できる。p.add_tools()でツールを追加。

from bokeh.models import HoverTool

world = world.rename(columns={'pop_est': '人口'})
source = GeoJSONDataSource(geojson=world.to_json())

p = figure(
        match_aspect=True,
        aspect_ratio=2.2,
        plot_width=800,
        tools='pan, wheel_zoom, tap'
)

hover_tool = HoverTool(
    tooltips=[('', '@name'),
              ('人口', '@{人口}{0,0}人'),
              ('緯度', '$y'),
              ('経度', '$x')],
    point_policy='follow_mouse'  # ホバーツールの表示位置
)

p.add_tools(hover_tool)
print(p.tools)

mapper = factor_cmap('name', magma(177), world['name'])

p.patches('xs', 'ys', source=source, fill_color=mapper, alpha=.6)

show(p)

[PanTool(id='2111', ...), WheelZoomTool(id='2112', ...), TapTool(id='2113', ...), HoverTool(id='2118', ...)]

ウィジェットを使う

jupyter上でウィジェットを使う場合、ipywidgetsとbokeh.models.widgetsの2つが選択肢になる。基本的にはどちらか一方を選んで使う。ipywidgets_bokehやjupyter_bokehを使うと、jupyterとbokehアプリケーションの両方で2つとも使えるらしいが、試していないのでここでは書かない。
ipywidgetsとbokeh.models.widgetsを使って人口とGDPの塗り分けを切り替えるドロップダウンを作る。書き方は異なるが、ドロップダウンのオブジェクトを作り、コールバック関数を作り、登録するという流れは同じ。

# ipywidgets

from bokeh.io import push_notebook
from bokeh.plotting import figure, show
from bokeh.models import GeoJSONDataSource, ColorBar
from bokeh.palettes import linear_palette, Cividis256, magma
from bokeh.transform import log_cmap

import ipywidgets as widgets
from IPython.display import display

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

p = figure(
        title='人口',
        match_aspect=True,
        aspect_ratio=2.2,
        plot_width=800,
)

# 人口データのカラーマッパー
colors = linear_palette(Cividis256, 100)[::-1]  # 256色のリストを100色に減らす
mapper_pop = log_cmap('pop_est', colors, low=1, high=10**10)

# GDPデータのカラーマッパー
colors = magma(80)
mapper_gdp = log_cmap('gdp_md_est', colors, low=1, high=10**8)

# カラーバーは共通
color_bar = ColorBar(color_mapper=mapper_pop['transform'])

p.add_layout(color_bar, 'right')

source = GeoJSONDataSource(geojson=world.to_json())

patches = p.patches('xs', 'ys',
                    source=source,
                    fill_color=mapper_pop,
                    line_color='#cdcdcd',
                    line_width=0.3
                    )
# ドロップダウン
dropdown = widgets.Dropdown(options=['人口', 'GDP'],
                            value='人口',
                            layout={'width': '100px'}
                            )

def callback(event):
    if event['new'] == '人口':
        p.title.text = '人口'
        color_bar.color_mapper = mapper_pop['transform']
        patches.glyph.fill_color = mapper_pop
    else:
        p.title.text = 'GDP'
        color_bar.color_mapper = mapper_gdp['transform']
        patches.glyph.fill_color = mapper_gdp
    
    push_notebook(handle=t)
    
dropdown.observe(callback, names='value')

display(dropdown)
t = show(p, notebook_handle=True)

bokehアプリケーション形式

jupyter上で動かす場合は、引数にDocumentオブジェクト受け取る関数を定義する。この関数の中身はほぼそのままbokehサーバーを使うbokehアプリケーションに転用できる。doc.add_root()を使ってDocumentにプロットやウィジェットなど、一連の部品を全て放り込むことで、ブラウザに反映される。 add_root()したものはdoc.rootsで確認できる。
dropdown.on_click(callback)はドロップダウンのメニューがクリックされた時にcallbackを呼びだす。on_change(attr, callback)はモデルの属性に変化があった時にcallbackを呼びだす。今回は未使用。
column()で縦に並べるレイアウトを作れる。他にrow()、gridplot()、layout()がある。詳細は公式参照。作ったオブジェクトは.childrenで確認できる。
show()に定義した関数を引数で渡して実行すると、プロットが表示される。jupyterのポートがデフォルトの8888でない時は、notebook_url='http://localhost:8889'のように現在のポート番号を引数で渡す。

# bokehアプリケーション形式
from bokeh.plotting import figure, show
from bokeh.models import GeoJSONDataSource, ColorBar
from bokeh.palettes import linear_palette, Cividis256, magma
from bokeh.transform import log_cmap
from bokeh.layouts import column
import bokeh.models.widgets as bkw

import geopandas as gpd


def bkapp(doc):
    world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
    
    p = figure(
            title='人口',
            match_aspect=True,
            aspect_ratio=2.2,
            plot_width=800,
    )
    
    # 人口データのカラーマッパー
    colors = linear_palette(Cividis256, 100)[::-1]
    mapper_pop = log_cmap('pop_est', colors, low=1, high=10**10)

    # GDPデータのカラーマッパー
    colors = magma(80)
    mapper_gdp = log_cmap('gdp_md_est', colors, low=1, high=10**8)

    # カラーバーは共通
    color_bar = ColorBar(color_mapper=mapper_pop['transform'])

    p.add_layout(color_bar, 'right')

    source = GeoJSONDataSource(geojson=world.to_json())

    patches = p.patches('xs', 'ys',
                        source=source,
                        fill_color=mapper_pop,
                        line_color='#cdcdcd',
                        line_width=0.3
                        )
    # ドロップダウン
    dropdown = bkw.Dropdown(menu=['人口', 'GDP'], width=100)

    def callback(event):
        if event.item == '人口':
            p.title.text = '人口'
            color_bar.color_mapper = mapper_pop['transform']
            patches.glyph.fill_color = mapper_pop

        else:
            p.title.text = 'GDP'
            color_bar.color_mapper = mapper_gdp['transform']
            patches.glyph.fill_color = mapper_gdp

    dropdown.on_click(callback)  # コールバックの登録
    
    col = column(dropdown, p)
    doc.add_root(col)

show(bkapp)

マップタイルの地図を表示

プリセットのタイルサービスを使う

bokehではあらかじめいくつかのタイルサービスを利用できるようになっている。
get_provider関数を使ってタイルサービスを選ぶ。選べる名前は公式を参照。
rangeの範囲を指定する必要がある。x_range, y_rangeともにだいたい-2000万(-2e+7)から2000万(2e+7)の間。
figure()にx_axis_type='mercator'、y_axis_type='mercator'の引数を渡すと、軸の表示をいい感じにしてくれる。
p.add_tile()で追加。

from bokeh.tile_providers import get_provider

p = figure(
        x_range=(-2e+7, 2e+7),
        y_range=(-2e+7, 2e+7),
        x_axis_type='mercator',
        y_axis_type='mercator',
        tools='pan,wheel_zoom,reset',
        active_scroll='wheel_zoom'
)

tile_provider = get_provider('CARTODBPOSITRON')
p.add_tile(tile_provider)

show(p)

Bokehバージョン3でのタイルサービスの使い方

add_tile()を使う。

引数は文字列かxyzservicesライブラリのTileProviderインスタンス。文字列はxyzservicesが認識できる名前。TileProviderはxyzservices.providersで定義されているものか自分で作成して使う。

retinaキーワードで解像度を選べる。高解像度のタイル地図が提供されていれば有効。'CartoDB Positron retina'のように名前に含めることも可能。

p = figure(
        width=500, height=500,
        x_range=(1.35e+7, 1.65e+7), y_range=(2.8e+6, 6.1e+6),
        x_axis_type='mercator', y_axis_type='mercator',
)
p.add_tile('CartoDB Positron', retina=True)

show(p)

import xyzservices.providers as xyz

p = figure(
        width=500, height=500,
        x_range=(1.35e+7, 1.65e+7), y_range=(2.8e+6, 6.1e+6),
        x_axis_type='mercator', y_axis_type='mercator',
)
p.add_tile(xyz.OpenStreetMap.Mapnik)

show(p)

# 国土地理院の地図タイルを使う
from xyzservices import TileProvider

gsimap_standard = dict(
    url='https://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png',
    attribution='出典:国土地理院',
    html_attribution='出典:<a href="https://maps.gsi.go.jp/development/ichiran.html">国土地理院</a>',
    name='GSIMap.Standard',
)
provider = TileProvider(gsimap_standard)

p = figure(
    width=500, height=500,
    x_range=(1.52e+7, 1.57e+7), y_range=(4.1e+6, 4.6e+6),
    x_axis_type='mercator', y_axis_type='mercator'
)
p.grid.grid_line_color = None

p.add_tile(provider)

show(p)

グリフを重ねて表示

マップタイルの座標は緯度経度とは違うので、これまで使ってきた世界地図のデータを重ねて表示するには座標を変換する必要がある。そのままでは表示できないということではなく、横幅4000万の軸に横幅360の地図を表示することになるので世界が豆粒になる。
GeoDataFrameのto_crs()メソッドで変換できる。

2021/11/03 追記: to_crs()だけだと上手くいかないことがわかったので処理を追加。jupyter上では南極がおかしなことになり、ブラウザでは描画自体ができない。変換後の座標では緯度約85度以上は地図の範囲外となるので、変換の前にこの範囲外の部分を除いておく必要があるらしい。pyproj.CRSを使って返還後の範囲を取得し、geopandas.overlay()で範囲外を取り除く。

参考

• stackoverflowの記事
• メルカトル図法｜特徴とWebメルカトルについて解説

import pyproj
from shapely.geometry import box

crs = pyproj.CRS.from_epsg(3857)
bounds = crs.area_of_use.bounds
print(bounds)

epsg_3857_bounds = gpd.GeoDataFrame([box(*bounds)], columns=['geometry'], crs=4326)
clipped = gpd.overlay(world, epsg_3857_bounds, how='intersection')

epsg_3857 = clipped.to_crs('epsg:3857')
epsg_3857.crs

(-180.0, -85.06, 180.0, 85.06)

<Projected CRS: EPSG:3857>
Name: WGS 84 / Pseudo-Mercator
Axis Info [cartesian]:
- X[east]: Easting (metre)
- Y[north]: Northing (metre)
Area of Use:
- name: World - 85°S to 85°N
- bounds: (-180.0, -85.06, 180.0, 85.06)
Coordinate Operation:
- name: Popular Visualisation Pseudo-Mercator
- method: Popular Visualisation Pseudo Mercator
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

source = GeoJSONDataSource(geojson=epsg_3857.to_json())

p = figure(
        match_aspect=True,
        aspect_ratio=1.2,    
        x_axis_type='mercator',
        y_axis_type='mercator',
)

p.patches('xs', 'ys', source=source)

show(p)

japan = epsg_3857.loc[epsg_3857['name'] == 'Japan']
source = GeoJSONDataSource(geojson=japan.to_json())

p = figure(
        x_range=(1.35e+7, 1.65e+7),
        y_range=(2.8e+6, 6.1e+6),
        x_axis_type='mercator',
        y_axis_type='mercator'
)
p.grid.grid_line_color = None

tile_provider = get_provider('STAMEN_TERRAIN')
p.add_tile(tile_provider)

p.patches('xs', 'ys', source=source, alpha=0.5)
p.star(x=15558880, y=4256800, size=26, color='gold', alpha=0.9)

show(p)

Googleマップの表示

GMapOptionsで設定、figureの代わりにgmap関数を使ってプロットを作成する。Google APIキーが必要。

from bokeh.models import GMapOptions, Slider, Select, Spinner
from bokeh.plotting import gmap

map_options = GMapOptions(
    lat=35.6809591,  # 緯度経度で指定
    lng=139.7673068,
    zoom=10,  # ズームレベル
    map_type='roadmap',
)

p = gmap('GoogleAPIキー',
         map_options,
         title='Googleマップの表示',
         tools='pan, reset')

# 設定を操作するためのウィジェット作成
options = [
    'roadmap',
    'satellite',
    'terrain',
    'hybrid'
]

map_type = Select(
    options=options,
    value='roadmap',
    width=100,
    title='map_type'  
)
# Selectのvalue属性とGMapOptionsのmap_type属性をリンク
map_type.js_link('value', map_options, 'map_type')  

zoom = Spinner(
    low=1,
    high=21,
    step=1,
    value=10,
    width=100,
    title='zoom'
) 
zoom.js_link('value', map_options, 'zoom')

tilt = Slider(
    start=0,
    end=45,
    step=45,
    value=45,
    width=100,
    title='tilt'
)
tilt.js_link('value', map_options, 'tilt')

show(row(p, column(map_type, zoom, tilt, margin=10), background='grey'))

gmapでグリフを重ねる場合、座標は緯度経度で指定する。

japan = world.loc[world['name'] == 'Japan']
source = GeoJSONDataSource(geojson=japan.to_json())

map_options = GMapOptions(lat=37.5, lng=135, zoom=5,)

p = gmap('GoogleAPIキー',
         map_options)

p.patches('xs', 'ys', source=source, alpha=0.5)
p.star(x=139.7673068, y=35.6809591, size=26, color='gold', alpha=0.9)

show(p)

メッシュデータを重ねて表示

これまでと同様に描画できるが、データが大きいと結構重たくなる。to_json()にも時間がかかるようになる。

国土数値情報ダウンロードサイト（国土交通省）の平年値メッシュデータを使用して描画してみる。
16個のメッシュデータをひとつにまとめて降水量と座標の列だけ抜き出す。データ数は約8万。

precipitation = gpd.read_file('データ/パス.shp', encoding='cp932')

precipitation.crs = 'epsg:4612'
precipitation.to_crs('epsg:3857', inplace=True)

precipitation.head(3)

	1月	2月	3月	4月	5月	6月	7月	8月	9月	10月	11月	12月	geometry
0	566	608	1233	1305	1683	1952	1491	1391	2386	1687	960	496	POLYGON ((15250770.239 4118674.237, 15250770.2...
1	563	608	1244	1304	1679	1952	1492	1380	2387	1661	959	495	POLYGON ((15250770.239 4119802.185, 15250770.2...
2	553	607	1243	1298	1673	1952	1494	1349	2388	1643	959	489	POLYGON ((15250770.239 4120930.246, 15250770.2...

from bokeh.models import Select, NumeralTickFormatter
from bokeh.transform import linear_cmap

def bkapp(doc):
    data = precipitation.loc[:, ['1月', 'geometry']]
    data = data.rename(columns={'1月': 'precipitation'})
    source = GeoJSONDataSource(geojson=data.to_json())
    
    months = [f'{i}月' for i in range(1, 13)]

    p = figure(
        width=800,
        aspect_ratio=5/4,
        x_range=(15180000, 15720000),
        y_range=(4100000, 4500000),
        x_axis_type='mercator',
        y_axis_type='mercator',
        title='月別降水量平年値（平成24年）',
    )

    tile_provider = get_provider('STAMEN_TERRAIN')
    p.add_tile(tile_provider)

    mapper = linear_cmap('precipitation', 'Inferno256', low=0, high=6000)
    p.patches('xs', 'ys',
              source=source,
              color=mapper,
              alpha=0.4)
    
    color_bar = ColorBar(color_mapper=mapper['transform'], width=15)
    color_bar.formatter = NumeralTickFormatter(format='0a')
    p.add_layout(color_bar, 'right')
    
    def callback(attr, old, new):
        data['precipitation'] = precipitation[new]
        source.geojson = data.to_json()
        
    select = Select(options=months, value='1月')
    select.on_change('value', callback)
    
    doc.add_root(column(select, p))
    
show(bkapp)

描画も動きも重たいしデータの入れ替えにも時間がかかる。

メッシュデータでデータの数が多い時はquadを使った方がいいかもしれない。quadは四隅を指定して四角形を描画するメソッドで、figure()の引数にoutput_backend='webgl'を渡すことで大きなデータの描画を改善できる。効果のあるグリフは限られているので公式
を参照。quadのサポートはバージョン2.4.3から。

from bokeh.models import ColumnDataSource

def bkapp(doc):
    # 座標データをquadに使える形に変換
    data = precipitation.loc[:, ['1月', 'geometry']]
    data = data.rename(columns={'1月': 'precipitation'})
    data[['left', 'bottom', 'right', 'top']] = data.bounds
    data = data.drop(columns=['geometry'])
    
    months = [f'{i}月' for i in range(1, 13)]
    
    # ColumnDataSourceを使う
    source = ColumnDataSource(data)

    p = figure(
        width=800,
        aspect_ratio=5/4,
        x_range=(15180000, 15720000),
        y_range=(4100000, 4500000),
        x_axis_type='mercator',
        y_axis_type='mercator',
        title='月別降水量平年値（平成24年）',
        output_backend='webgl',  # 足すだけ
    )

    tile_provider = get_provider('STAMEN_TERRAIN')
    p.add_tile(tile_provider)

    mapper = linear_cmap('precipitation', 'Inferno256', low=0, high=6000)
    
    patches = p.quad('left', 'right', 'top', 'bottom',
                     source=source,
                     color=mapper,
                     alpha=0.4)

    color_bar = ColorBar(color_mapper=mapper['transform'], width=15)
    color_bar.formatter = NumeralTickFormatter(format='0a')
    p.add_layout(color_bar, 'right')
    
    def callback(attr, old, new):
        data['precipitation'] = precipitation[new]
        source.data = data
        
    select = Select(options=months, value='1月')
    select.on_change('value', callback)
    
    doc.add_root(column(select, p))

show(bkapp)