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

bokeh学習の記録。

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

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

# 最初に

## 使うもの
- bokeh
- geopandas

## 参考サイト

- [bokeh](https://docs.bokeh.org/en/2.4.3/index.html)
- [GeoPandas](https://geopandas.org/index.html)
- [Shapely](https://shapely.readthedocs.io/en/stable/manual.html)

* [ipywidgetsとBokeh使ってインタラクティブな可視化をする](https://qiita.com/driller/items/0730325bf5c1cd689979)

## ページ内のコード

- 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の機能を使ってマップタイルを表示する。


```python
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()
```

<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>pop_est</th>
      <th>continent</th>
      <th>name</th>
      <th>iso_a3</th>
      <th>gdp_md_est</th>
      <th>geometry</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>920938</td>
      <td>Oceania</td>
      <td>Fiji</td>
      <td>FJI</td>
      <td>8374.0</td>
      <td>MULTIPOLYGON (((180.00000 -16.06713, 180.00000...</td>
    </tr>
    <tr>
      <th>1</th>
      <td>53950935</td>
      <td>Africa</td>
      <td>Tanzania</td>
      <td>TZA</td>
      <td>150600.0</td>
      <td>POLYGON ((33.90371 -0.95000, 34.07262 -1.05982...</td>
    </tr>
    <tr>
      <th>2</th>
      <td>603253</td>
      <td>Africa</td>
      <td>W. Sahara</td>
      <td>ESH</td>
      <td>906.5</td>
      <td>POLYGON ((-8.66559 27.65643, -8.66512 27.58948...</td>
    </tr>
    <tr>
      <th>3</th>
      <td>35623680</td>
      <td>North America</td>
      <td>Canada</td>
      <td>CAN</td>
      <td>1674000.0</td>
      <td>MULTIPOLYGON (((-122.84000 49.00000, -122.9742...</td>
    </tr>
    <tr>
      <th>4</th>
      <td>326625791</td>
      <td>North America</td>
      <td>United States of America</td>
      <td>USA</td>
      <td>18560000.0</td>
      <td>MULTIPOLYGON (((-122.84000 49.00000, -120.0000...</td>
    </tr>
  </tbody>
</table>
</div>

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

# 地図の表示や設定

## とにかく地図を表示する

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

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

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

show(p)
```

![map_01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/f8d19692-2a91-2247-9423-0626a81fa254.png)

## 見栄えを整える

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

```python
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)
```

![map_02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/5239ceb1-b54f-eace-1dec-93b74323d873.png)

# さらに調整

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

```python
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)
```

![map_03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/0c194c67-14eb-ec01-698c-67447e7973b7.png)

# 点を重ねて表示する

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

```python
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
    
![map_04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/4c1ed20e-9036-beba-5372-7ba5618be303.png)

# 表示範囲を指定する

`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を使う意味がないと気付いたのでコードを修正。

```python
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))
```

![map_05.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/d5f8ff56-7ae8-4719-ecb4-519950c20c25.png)

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

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

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

```python
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'))
```

![map_06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/054f6338-64d8-7331-c5a9-d8dbb93de8f6.png)

## 値に応じて塗り分ける

`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'は描画範囲内、それ以外は縁の位置。  


```python
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)
```

![map_07.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/ebf29c38-6ef0-02b1-03f7-ad5f1e1f9eb5.png)

# ホバーツールを使う

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

## ツールについて

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

```python
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', ...)]
    
![map_08.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/35f7dc39-1a0c-861d-2989-46ecd93f9bf5.png)

# ウィジェットを使う

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

```python
# 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)
```

![jupyter_dropdown.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/b8d0fd2c-7179-a4fb-76c0-1e3bdb24da62.gif)

## 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()`がある。詳細は[公式](https://docs.bokeh.org/en/2.3.3/docs/user_guide/layout.html)参照。作ったオブジェクトは`.children`で確認できる。  
`show()`に定義した関数を引数で渡して実行すると、プロットが表示される。jupyterのポートがデフォルトの8888でない時は、`notebook_url='http://localhost:8889'`のように現在のポート番号を引数で渡す。

```python
# 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_dropdown.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/5991006a-2ec1-ce04-1ba7-97e5c86e4e1d.gif)

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

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

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


```python
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)
```

![map_09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/9c3e46fc-e189-fc94-7f5b-67a358feb0c5.png)

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

add_tile()を使う。

引数は文字列か[xyzservices](https://xyzservices.readthedocs.io/en/stable/)ライブラリのTileProviderインスタンス。文字列は[xyzservicesが認識できる名前](https://xyzservices.readthedocs.io/en/stable/introduction.html)。TileProviderはxyzservices.providersで定義されているものか自分で作成して使う。

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


```python
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)
```

![tile_map_01.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/96ab3f10-0c0a-a643-db49-b8e1e4c607b5.jpeg)

```python
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)
```

![tile_map_02.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/bd3ccfa4-39bb-57f2-2682-4f121af9d45c.jpeg)

```python
# 国土地理院の地図タイルを使う
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)
```

![tile_map_03.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/6d736ace-530c-2ca9-0f25-ad0b95270185.jpeg)


## グリフを重ねて表示

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

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

参考 
- [stackoverflowの記事](https://stackoverflow.com/questions/61091738/how-to-remove-shapes-that-cause-a-problem-upon-reprojection-from-a-geopandas-dat)
- [メルカトル図法｜特徴とWebメルカトルについて解説](https://club.informatix.co.jp/?p=1183)


```python
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

```python
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)
```

![map_10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/2a91e137-8395-d8e1-a37c-d4470a54fa2c.png)

```python
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)
```

![map_11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/8bbae512-0152-0699-2422-42f106cb86a4.png)

## Googleマップの表示

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

```python
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'))
```

![map_12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/b1992794-aec5-2c0d-65fd-13e5ff500941.png)

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

```python
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)
```

![map_13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/cdb5b979-9aad-b227-cfd8-66c17e13b373.png)

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

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

[国土数値情報ダウンロードサイト（国土交通省）](https://nlftp.mlit.go.jp/ksj/index.html)の[平年値メッシュデータ](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-G02.html)を使用して描画してみる。
16個のメッシュデータをひとつにまとめて降水量と座標の列だけ抜き出す。データ数は約8万。

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

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

precipitation.head(3)
```

<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>1月</th>
      <th>2月</th>
      <th>3月</th>
      <th>4月</th>
      <th>5月</th>
      <th>6月</th>
      <th>7月</th>
      <th>8月</th>
      <th>9月</th>
      <th>10月</th>
      <th>11月</th>
      <th>12月</th>
      <th>geometry</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>566</td>
      <td>608</td>
      <td>1233</td>
      <td>1305</td>
      <td>1683</td>
      <td>1952</td>
      <td>1491</td>
      <td>1391</td>
      <td>2386</td>
      <td>1687</td>
      <td>960</td>
      <td>496</td>
      <td>POLYGON ((15250770.239 4118674.237, 15250770.2...</td>
    </tr>
    <tr>
      <th>1</th>
      <td>563</td>
      <td>608</td>
      <td>1244</td>
      <td>1304</td>
      <td>1679</td>
      <td>1952</td>
      <td>1492</td>
      <td>1380</td>
      <td>2387</td>
      <td>1661</td>
      <td>959</td>
      <td>495</td>
      <td>POLYGON ((15250770.239 4119802.185, 15250770.2...</td>
    </tr>
    <tr>
      <th>2</th>
      <td>553</td>
      <td>607</td>
      <td>1243</td>
      <td>1298</td>
      <td>1673</td>
      <td>1952</td>
      <td>1494</td>
      <td>1349</td>
      <td>2388</td>
      <td>1643</td>
      <td>959</td>
      <td>489</td>
      <td>POLYGON ((15250770.239 4120930.246, 15250770.2...</td>
    </tr>
  </tbody>
</table>
</div>

```python
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)
```

![map_14_01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/bbbbad0b-3be4-62c7-7dcb-2779fa84668d.png)

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

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

```python
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)
```

![mesh_quad.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1714033/efb90df1-03ca-3079-172d-d116c8417b6d.gif)