Uber H3のRインタフェース{h3forr}を試す

緯度・経度といった座標情報を集計して可視化するには総務省が定めた標準地域メッシュに加工することが多いが、Uberは矩形ではなく六角形グリッドを使用しているとのこと。

Uberはなぜ地図を六角形で埋めるのか - Qiita

ドコモのAIタクシーは矩形のメッシュだが中国の配車サービスDiDiでも六角形グリッドが使用されているようで六角形グリッドの方が見た目がカッコいい。

Uberの六角形グリッドシステムH3はオープンソースで公開されており、Rのインタフェースについてもいくつか公開されている。

ここでは六角形グリッドを{sf}のポリゴンに変換までしてくれる{h3forr}を試してみる。

メッシュとの比較

まずは東京タワー(139.745433, 35.658581)付近でH3グリッドと標準地域メッシュを比較してみる。緯度・経度の座標からH3のコードを取得するにはh3forr::geo_to_h3()を使用する。デフォルトでグリッドの解像度は7だが以下では7から10で取得している。 H3のコードから六角形の各座標を取得するにはh3forr::h3_to_geo_boundary()を使用する。sfオブジェクトにするにはh3forr::geo_boundary_to_sf()を使用する。

library(h3forr)
library(dplyr)
library(purrr)
library(tidyr)
library(sf)

grids <- tibble(lat = 35.658581,
                lng = 139.745433,
                res = seq(7, 10, 1)) %>%
  mutate(h3_index = pmap(.l = list(lat, lng, res),
                         function(x, y, z) {
                           geo_to_h3(coords = c(x, y), res = z)
                         })) %>%
  unnest() %>%
  mutate(boundary = h3_to_geo_boundary(h3_index))

grids <-
  st_sf(data.frame(grids, geom = geo_boundary_to_sf(grids$boundary)))

glimpse(grids)
Observations: 4
Variables: 6
$ lat      <dbl> 35.65858, 35.65858, 35.65858, 35.65858
$ lng      <dbl> 139.7454, 139.7454, 139.7454, 139.7454
$ res      <dbl> 7, 8, 9, 10
$ h3_index <chr> "872f5aad8ffffff", "882f5aadc7fffff", "892f5aadc6fffff", "8a2f5aadc6dffff"
$ boundary <list> [<matrix[7 x 2]>, <matrix[7 x 2]>, <matrix[7 x 2]>, <matrix[7 x 2]>]
$ geometry <POLYGON [°]> POLYGON ((139.7487 35.6602,..., POLYGON ((139.7529 35.65271..., POLYGO…

大きさを比較するために標準地域メッシュも用意する。

library(jpmesh)
meshes <- tibble(
  lat = 35.658581,
  lng = 139.745433,
  size = c("1km", "500m", "250m", "125m")
) %>%
  mutate(code = pmap(.l = list(lng, lat, size),
                     function(x, y, z) {
                       coords_to_mesh(x, y, mesh_size = z)
                     })) %>%
  unnest()

meshes <- export_meshes(meshes$code)
glimpse(meshes)
Observations: 4
Variables: 2
$ meshcode <chr> "53393599", "533935992", "5339359921", "53393599212"
$ geometry <POLYGON [°]> POLYGON ((139.7375 35.65833..., POLYGON ((139.7438 35…

{leaflet}で可視化してみると以下のような感じになる。

library(leaflet)
grids %>%
  leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(weight = 2, fillColor = "Transparent") %>%
  addPolygons(
    data = meshes,
    weight = 2,
    fillColor = "Transparent",
    color = "red"
  )

f:id:tak95:20190501181757p:plain:w600 — compare grid to mesh

青線がH3のグリッド、赤線が標準地域メッシュを表していて125mメッシュと解像度10のH3グリッドが大体同じくらいの大きさとなる。デフォルトの解像度7だとかなり広く見える。面積で比較すると解像度7のH3グリッドは1kmメッシュの4倍以上も大きい。

library(lwgeom)
meshes %>% st_area() %>% knitr::kable()

x
1047333.14 [m²]
259745.28 [m²]
65092.88 [m²]
16717.90 [m²]

grids %>% st_area() %>% knitr::kable()

x
4743067.20 [m²]
677739.26 [m²]
96815.85 [m²]
13830.71 [m²]

NYC OpenData 2016 Green Taxi Trip Dataの可視化

ある点だけをグリッドにしても面白くないので以下のデータで可視化してみる。

https://data.cityofnewyork.us/Transportation/2016-Green-Taxi-Trip-Data/hvrh-b6nb

データの取得

上記データセットはCSVにエクスポートできるようなのだがファイルサイズが大きく、自宅環境では途中でキャンセルとなってしまう。一方でAPIも用意されており、R用のパッケージ{RSocrata}まで提供されている。以下では2016-02-01の1日分を取得している。

library(RSocrata)
library(lubridate)
df <- read.socrata(
  "https://data.cityofnewyork.us/resource/pqfs-mqru.json?$where=date_trunc_ymd(lpep_pickup_datetime)='2016-02-01'"
) %>%
  mutate(
    pickup_longitude = as.numeric(pickup_longitude),
    pickup_latitude = as.numeric(pickup_latitude),
    hour = hour(lpep_pickup_datetime),
    h3_index = pmap(.l = list(pickup_latitude, pickup_longitude),
                    function(x, y) {
                      geo_to_h3(coords = c(x, y), res = 8)
                    })
  ) %>%
  select(lpep_pickup_datetime,
         hour,
         pickup_longitude,
         pickup_latitude,
         h3_index) %>% 
  # 緯度・経度が正しくなさそうなデータを除外
  filter(pickup_latitude >= 40.48333 & pickup_latitude <= 45)

glimpse(df)
Observations: 42,366
Variables: 5
$ lpep_pickup_datetime <dttm> 2016-02-01 00:00:01, 2016-02-01 00:01:33, 2016-0…
$ hour                 <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ pickup_longitude     <dbl> -73.93902, -73.89149, -73.98378, -73.80752, -73.9…
$ pickup_latitude      <dbl> 40.80521, 40.74665, 40.67613, 40.70037, 40.74493,…
$ h3_index             <list> ["882a1008d9fffff", "882a100c69fffff", "882a1077…

集計 & ポリゴンの作成

時間とグリッドごとに乗車が何回発生したか集計を行う。

df.aggr <- df %>% 
  unnest() %>% 
  group_by(hour, h3_index) %>% 
  count() %>% 
  ungroup() %>% 
  mutate(boundary = h3_to_geo_boundary(h3_index))

sf <- st_sf(data.frame(df.aggr, geom = geo_boundary_to_sf(df.aggr$boundary)))

glimpse(sf)
Observations: 5,208
Variables: 5
$ hour     <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ h3_index <chr> "882a100101fffff", "882a10011bfffff", "882a100127fffff", "882…
$ n        <int> 1, 1, 3, 5, 1, 2, 1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 9, 6, 4, 12, …
$ boundary <list> [<matrix[7 x 2]>, <matrix[7 x 2]>, <matrix[7 x 2]>, <matrix[…
$ geometry <POLYGON [°]> POLYGON ((-73.85562 40.8794..., POLYGON ((-73.82913 4…

可視化

{leaflet}で18時を可視化してみる。マーカーは自由の女神。

sf %>% 
  filter(hour == 18) %>% 
  leaflet() %>% 
  setView(lng = -73.94506, lat = 40.74829, zoom = 11) %>% 
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(weight = 1, fillColor = ~pal(n), fillOpacity = 0.5) %>% 
  addMarkers(lng = -74.0445, lat = 40.6892) %>% 
  addLegend("bottomright", pal = pal, values = ~n,
            title = "number of pickup",
            opacity = 1
  )