LoginSignup
4
2

More than 1 year has passed since last update.

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

Last updated at Posted at 2021-07-18

bokeh学習の記録。

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

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

最初に

使うもの

  • bokeh
  • geopandas

参考サイト

ページ内のコード

  • 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)

map_01.png

見栄えを整える

形が崩れてしまっているので、ぱっと見違和感のないようにする。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)

map_02.png

さらに調整

形は整ったので、色を変えたり、サイズを変えたりしてみる。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)

map_03.png

点を重ねて表示する

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

map_04.png

表示範囲を指定する

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

形を崩さないようにするには自分で縦横の長さを設定する必要がある。
中央のフレーム部分の縦横を指定するframe_heightframe_widthを使うのが簡単だがsizing_modeとの相性が悪い。sizing_mode'fixed'以外にすると思った通りの形にならない可能性がある。plot_width, plot_heightaspect_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))

map_05.png

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

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

ソースに色の値を持つ列を加えるか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'))

map_06.png

値に応じて塗り分ける

linear_cmap()またはlog_cmap()を使う。第一引数から適用するソースの列名、色のリスト。lowhighで範囲を指定する。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)

map_07.png

ホバーツールを使う

HoverToolクラスを使う。tooltipsがホバーツールの中身になる。渡すのはliststringlistの値は(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', ...)]

map_08.png

ウィジェットを使う

jupyter上でウィジェットを使う場合、ipywidgetsbokeh.models.widgetsの2つが選択肢になる。基本的にはどちらか一方を選んで使う。ipywidgets_bokehjupyter_bokehを使うと、jupyterとbokehアプリケーションの両方で2つとも使えるらしいが、試していないのでここでは書かない。
ipywidgetsbokeh.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)

jupyter_dropdown.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()がある。詳細は公式参照。作ったオブジェクトは.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_dropdown.gif

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

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

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)

map_09.png

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)

tile_map_01.jpg

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

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

グリフを重ねて表示

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

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

参考

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)

map_10.png

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

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

map_12.png

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)

map_13.png

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

これまでと同様に描画できるが、データが大きいと結構重たくなる。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)

map_14_01.png

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

メッシュデータでデータの数が多い時は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)

mesh_quad.gif

4
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
4
2