大量のデータを画像化して描写するDatashader使ってみる


大規模なデータをラスタライズし、可視化・要約を行う[Holoviz](https://holoviz.org/)の[Datashader](https://datashader.org/#)で地図描写試してみたという話。

**この記事の要約**

- datashaderでは、ver0.16からgeopandasを読みこむことが出来る
- 画面のピクセルサイズを先に決めてしまい、データをピクセルに対応付け・集計する
- plotlyで利用する時は、画像が対応する緯度経度を取り出してmapboxレイヤー作る

# これは何するの？

下のようなことができます。

- **公式トップページ**

>With code just like the above, you can plot 300 million points of data (one per person in the USA) from the 2010 census without any parameter tuning:
（訳）上記のようなコードを使用するとパラメーターを調整せずに、2010年の国勢調査から3億ポイントのデータ(米国の一人ひとり) をプロットできます。
![](https://datashader.org/assets/images/usa_census.jpg)
https://datashader.org/#


- **[dash-world-cell-towers](https://dash.gallery/dash-world-cell-towers/)([github](https://github.com/plotly/dash-world-cell-towers))**

>This dataset contains over 40 million rows of data, each row including several cell tower features in addition to the latitude/longitude location of the tower. To handle this fairly large dataset, this app makes use of Dask for parallel processing, and Datashader for server side rendering.
（訳）このデータセットには4000万行を超えるデータが含まれており、各行には基地局の緯度/経度の位置に加えて、いくつかの付随する情報が含まれています。このかなり大きなデータセットを処理するために、このアプリは並列処理にDaskを使用し、サーバー側のレンダリングにDatashaderを使用します。
![](https://github.com/plotly/dash-world-cell-towers/raw/master/images/full_preview.png)
https://github.com/plotly/dash-world-cell-towers


やっていることとしては、画面のピクセルサイズを先に決めてしまい、データをピクセルに対応付け・集計を行っているみたいです。

集約関数は複数用意されており、目的に応じて利用できます。

↓この辺に書いてあります。

https://datashader.org/getting_started/Pipeline.html

https://datashader.org/api.html#reductions


# Plotlyのmapbox系プロットでとりあえず使ってみる

公式の説明の通り。

> Plotly and Datashader in PythonMapbox Map Layers in Python 
How to use datashader to rasterize large datasets, and visualize the generated raster data with plotly. 
https://plotly.com/python/datashader/

- [Plotly and Datashader in Python](https://plotly.com/python/datashader/)
- [xarray.Coordinates](https://docs.xarray.dev/en/stable/generated/xarray.Coordinates.html#xarray.Coordinates)

datashaderでは、ver0.16からgeopandasを読みこむことが出来る。

そのver0.16がpython 3.8などまだ割と使われているバージョンに対応していないことに留意。

https://datashader.org/user_guide/Geopandas.html

## kaggle Taxi Trajectory Dataを可視化してみる

https://www.kaggle.com/datasets/crailtap/taxi-trajectory

このデータは、ポルトガルを走るタクシーの軌跡データで、約170万！！の移動軌跡が収録されています。

描写対象として、軌跡の長さが、2より大きい166万行、点の数に直すと、8341万点に対して描写しています。

前処理で死にそうなほどメモリ食ってますが、datashaderの処理自体は割とサクサク動きます。

描写してしまえば、画像なので言わずもがなです。

### 前処理


~~~python

import polars as pl
import geopandas as gpd
import shapely

import datashader as ds
from colorcet import fire
import plotly.express as px

from IPython.display import Image
~~~
~~~python
df = pl.read_csv(
    'TaxiTrajectoryData.csv',
    schema=dict(
        TRIP_ID=pl.Int64,
        CALL_TYPE=pl.Utf8,
        ORIGIN_CALL=pl.Utf8,
        ORIGIN_STAND=pl.Int16,
        TAXI_ID=pl.Int32,
        TIMESTAMP=pl.Int64,
        DAY_TYPE=pl.Utf8,
        MISSING_DATA=pl.Boolean,
        POLYLINE=pl.Utf8
    )
)
df
~~~

| TRIP_ID             | CALL_TYPE     | ORIGIN_CALL  | ORIGIN_STAND | TAXI_ID  | TIMESTAMP  | DAY_TYPE      | MISSING_DATA | POLYLINE                              |
|---------------------|---------------|--------------|--------------|----------|------------|---------------|--------------|---------------------------------------|
| i64                 | str           | str          | i16          | i32      | i64        | str           | bool         | str                                   |
| 1372636858620000589 | &quot;C&quot; | &quot;&quot; | null         | 20000589 | 1372636858 | &quot;A&quot; | false        | &quot;[[-8.618643,41.141412],[-8.618… |
| 1372637303620000596 | &quot;B&quot; | &quot;&quot; | 7            | 20000596 | 1372637303 | &quot;A&quot; | false        | &quot;[[-8.639847,41.159826],[-8.640… |
| 1372636951620000320 | &quot;C&quot; | &quot;&quot; | null         | 20000320 | 1372636951 | &quot;A&quot; | false        | &quot;[[-8.612964,41.140359],[-8.613… |
| 1372636854620000520 | &quot;C&quot; | &quot;&quot; | null         | 20000520 | 1372636854 | &quot;A&quot; | false        | &quot;[[-8.574678,41.151951],[-8.574… |
| 1372637091620000337 | &quot;C&quot; | &quot;&quot; | null         | 20000337 | 1372637091 | &quot;A&quot; | false        | &quot;[[-8.645994,41.18049],[-8.6459… |
| &hellip;            | &hellip;      | &hellip;     | &hellip;     | &hellip; | &hellip;   | &hellip;      | &hellip;     | &hellip;                              |
| 1404171463620000698 | &quot;C&quot; | &quot;&quot; | null         | 20000698 | 1404171463 | &quot;A&quot; | false        | &quot;[[-8.612469,41.14602],[-8.6124… |
| 1404171367620000670 | &quot;C&quot; | &quot;&quot; | null         | 20000670 | 1404171367 | &quot;A&quot; | false        | &quot;[[-8.610138,41.140845],[-8.610… |
| 1388745716620000264 | &quot;C&quot; | &quot;&quot; | null         | 20000264 | 1388745716 | &quot;A&quot; | false        | &quot;[]&quot;                        |
| 1404141826620000248 | &quot;B&quot; | &quot;&quot; | 12           | 20000248 | 1404141826 | &quot;A&quot; | false        | &quot;[[-8.630712,41.154885],[-8.630… |
| 1404157147620000079 | &quot;B&quot; | &quot;&quot; | 34           | 20000079 | 1404157147 | &quot;A&quot; | false        | &quot;[[-8.615538,41.140629],[-8.615… |


~~~python
# 死ぬほどメモリ喰うので注意

# ポイント＝行として変換する例
#ddf = df\
#     .drop(['ORIGIN_STAND','DAY_TYPE'])\
#     .with_columns(
#         (pl.col('TIMESTAMP') * 10 ** 6).cast(pl.Datetime),
#         pl.col('POLYLINE').str.json_decode()
#     )\
#     .with_columns(
#         pl.int_ranges(0,pl.col('POLYLINE').list.len()).alias('SEC')
#     )\
#     .explode(['POLYLINE','SEC'])\
#     .with_columns(
#         pl.col('TIMESTAMP')+pl.duration(seconds=(pl.col('SEC') * 15))
#     )

# ラインストリングにとして扱うための下準備
ddf = df\
    .drop(['ORIGIN_STAND','DAY_TYPE','MISSING_DATA'])\
    .with_columns(
        (pl.col('TIMESTAMP') * 10 ** 6).cast(pl.Datetime),
        pl.col('POLYLINE').str.json_decode()
    )\
    .filter(pl.col('POLYLINE').list.len() > 2)
~~~

~~~python
%%time
# 死ぬほどメモリ喰うので注意

t=ddf.to_pandas()
t['POLYLINE'] = t['POLYLINE'].apply(lambda x: shapely.LineString(x))

gdf = gpd.GeoDataFrame(t,geometry='POLYLINE')

# ポイントとして整形する例
# a=ddf.to_pandas()
# gdf = gpd.GeoDataFrame(
#     a,
#     geometry=gpd.points_from_xy(
#             a['POLYLINE'].apply(lambda x: x[0]),
#             a['POLYLINE'].apply(lambda x: x[1])
#     )
# )

~~~

CPU times: user 2min 42s, sys: 8.94 s, total: 2min 51s
Wall time: 2min 51s

~~~python
gdf
~~~

|         | TRIP_ID             | CALL_TYPE | ORIGIN_CALL | TAXI_ID  | TIMESTAMP           | POLYLINE                                          |
|---------|---------------------|-----------|-------------|----------|---------------------|---------------------------------------------------|
| 0       | 1372636858620000589 | C         |             | 20000589 | 2013-07-01 00:00:58 | LINESTRING (-8.61864 41.14141, -8.6185 41.1413... |
| 1       | 1372637303620000596 | B         |             | 20000596 | 2013-07-01 00:08:23 | LINESTRING (-8.63985 41.15983, -8.64035 41.159... |
| 2       | 1372636951620000320 | C         |             | 20000320 | 2013-07-01 00:02:31 | LINESTRING (-8.61296 41.14036, -8.61338 41.140... |
| 3       | 1372636854620000520 | C         |             | 20000520 | 2013-07-01 00:00:54 | LINESTRING (-8.57468 41.15195, -8.5747 41.1519... |
| 4       | 1372637091620000337 | C         |             | 20000337 | 2013-07-01 00:04:51 | LINESTRING (-8.64599 41.18049, -8.64595 41.180... |
| ...     | ...                 | ...       | ...         | ...      | ...                 | ...                                               |
| 1666761 | 1388660427620000585 | C         |             | 20000585 | 2014-01-02 11:00:27 | LINESTRING (-8.60697 41.16228, -8.60697 41.162... |
| 1666762 | 1404171463620000698 | C         |             | 20000698 | 2014-06-30 23:37:43 | LINESTRING (-8.61247 41.14602, -8.61249 41.145... |
| 1666763 | 1404171367620000670 | C         |             | 20000670 | 2014-06-30 23:36:07 | LINESTRING (-8.61014 41.14084, -8.61017 41.140... |
| 1666764 | 1404141826620000248 | B         |             | 20000248 | 2014-06-30 15:23:46 | LINESTRING (-8.63071 41.15488, -8.63073 41.154... |
| 1666765 | 1404157147620000079 | B         |             | 20000079 | 2014-06-30 19:39:07 | LINESTRING (-8.61554 41.14063, -8.61542 41.140... |

### 描写

重要なのは、`agg = canvas.line(source=gdf, geometry="POLYLINE",agg=ds.count())`と`img = ds.tf.shade(agg, cmap=fire)[::-1].to_pil()`の部分。
`agg=ds.count()`でどうゆう集約をするか？どの列を参照するか？決定する。この例だと、「何回ピクセル上にマッピングされたか？」のカウント。

 参考: [Datashader - api #reductions](https://datashader.org/api.html#reductions)

~~~python
# From https://plotly.com/python/datashader/
canvas = ds.Canvas(plot_width=2000, plot_height=2000)
agg = canvas.line(source=gdf, geometry="POLYLINE",agg=ds.count())

# agg is an xarray object, see http://xarray.pydata.org/en/stable/ for more details
coords_lat, coords_lon = agg.coords['y'].values, agg.coords['x'].values

# Corners of the image, which need to be passed to mapbox
# 要するに上の行と合わせ、作成した画像がどのに位置するのか？エリア（lon,lat）と画像の頂点（角）の対応を取るために必要
coordinates = [[coords_lon[0], coords_lat[0]],
               [coords_lon[-1], coords_lat[0]],
               [coords_lon[-1], coords_lat[-1]],
               [coords_lon[0], coords_lat[-1]]]

# import datashader.transfer_functions as tf
img = ds.tf.shade(agg, cmap=fire)[::-1].to_pil()

# Trick to create rapidly a figure with mapbox axes
# トリックを使用しないとこうなるらしい https://github.com/plotly/dash-world-cell-towers/blob/18fe0d3e2966a3e0b8b8fb03fe40e8daae565687/dash_opencellid/app.py#L720
x,y = gdf.POLYLINE.iloc[0].xy
fig = px.scatter_mapbox(
    gdf.iloc[0],
    lon=x,
    lat=y,
    zoom=7,
    # height=800,
    # width=800
    opacity=0,

    )
# Add the datashader image as a mapbox layer image

fig.update_layout(
                 mapbox_style="carto-darkmatter",
                # mapbox_style="carto-positron",
                 mapbox_layers = [{
                    "sourcetype": "image",
                    "source": img,
                    "coordinates": coordinates
                }]
)
fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0})
# fig.show()
Image(fig.to_image())
~~~

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3472791/4dc9ec99-3fd8-b99a-d777-d8d4959c05c5.png)

https://gist.github.com/matsu2254/d2497b8e016bb7c4fb6e1c8161ea81e1


## 町丁ポリゴンを塗り分けてみる

山形県の全域（！！）の町丁目ポリゴンを適当に塗り分けてプロットします。

4428のポリゴンがあります。数としては少ないですが、複雑な形状も含まれます。普通に描写したらplotlyだとまともに動かんです…

ちなみに、乱数で塗り分けているので、隣り合うポリゴンが違う色になる保証はないです。

(N100・12GBメモリのノートでも難なく動いたので、試すにはこちらのほうが向いているかも…)

データの出典： [e-stat - 境界データダウンロード 山形県全域](https://www.e-stat.go.jp/gis/statmap-search?page=1&type=2&aggregateUnitForBoundary=A&toukeiCode=00200521&toukeiYear=2020&serveyId=A002005212020&prefCode=06&coordsys=1&format=shape&datum=2000)

ほとんど同じなのでGistのみ貼っておきます。（こっちは、`agg=ds.by('r')`のように集約しています）

↓こんな具合に描写できる

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3472791/dd3df971-f9ff-9b90-3f58-5130d4587fbc.png)

https://gist.github.com/matsu2254/48f1e2aea1a9115b7ed6893d7f9fd0d1#file-datashader_plotly-ipynb

## HoloViewsでの描写

holoviewsで可視化する場合、データの座標系をWEBメルカトル法にする必要があるみたいです。

座標系の変換にはややこい計算が必要ですが、~~pyprojとか使うと簡単です~~geopandas.to_crsしたほうが簡単だとあとから気づいた人

（この例ではデータ量が多すぎるので1万件に絞っています）

その他は[チュートリアル](https://datashader.org/getting_started/Introduction.html)通り。

https://gist.github.com/matsu2254/506e1444dc3d1b5424ae3daf262203f8


## 参考

- [Holoviz - Datashader](https://datashader.org/#)
    - [Datashader - getting started Introduction](https://datashader.org/getting_started/Introduction.html)
    - [Datashader - getting started Pipeline](https://datashader.org/getting_started/Pipeline.html)
    - [Datashader - api #reductions](https://datashader.org/api.html#reductions)
    - [Qiita - 大量データで散布図を描いた時のオーバーラップをなんとかする（Matplotlib, Pandas, Datashader）](https://qiita.com/khigashi02/items/5197796711a316ae801f)
- 座標系関連
    - [pyproj Documentation](https://pyproj4.github.io/pyproj/stable/)
    - [Convert GPS coordinates to Web Mercator EPSG:3857 using python/pyproj](https://gis.stackexchange.com/questions/247871/convert-gps-coordinates-to-web-mercator-epsg3857-using-python-pyproj )
  - [Transformation Grids](https://pyproj4.github.io/pyproj/stable/transformation_grids.html)
  - [座標変換](http://hydro.iis.u-tokyo.ac.jp/~akira/page/Python/contents/pyproj/transform.html)
  - [OSM - Mercator](https://wiki.openstreetmap.org/wiki/Mercator#JavaScript)
- マップタイルの変更
  - [HoloViews - Tiles](https://holoviews.org/reference/elements/bokeh/Tiles.html)
  - [Plotly - Mapbox Map Layers in Python ](https://plotly.com/python/mapbox-layers/)
- [HoloViz - Colorcet](https://colorcet.holoviz.org/)
    - [Qiita - colorcetを用いて均一なカラーバーを作成する](https://qiita.com/shu32/items/e10e974220db3d59695c)