この記事は、BigQuery の Advent Calendar 2020 の7日目の記事です。 https://qiita.com/advent-calendar/2020/bigquery

私はGIS関係のベンチャー企業に就いて、そろそろ1年が立ちました。その中で学んできた事を含めてBigQueryで試してみようと思います。

GIS

GISとは、Geographic Information Systemの略で地理情報システムサービスです。 BigQueryGIS については、他の人の記事を読んで頂いた方が良いと思います。 https://qiita.com/shin_ishiguro/items/2429038b2c4c99c10837

地図にどのような情報を載せるかを定義するものです。GoogleMapですね。

どんなデータでも当然ですが、データ量が多くなるとそのままWebサービスに表示しても、重すぎてWebサービスには簡単にはできません。 GISの場合、ズームインすると小さいエリアの情報を出すことでデータ量は小さくなるのですが、ズームアウトしてエリアが広くなるとそれだけ多くの情報を必要になります。 GoogleMapのマーカーを大量に表示しても、利用者は扱いにくいし重すぎて表示が出来なくなると思います。

その場合に、データをクラスタ化してデータを削減する必要があります。

クラスタの手法

Geo情報の中でも今回はポイントを利用します。

ポイントが大量にあるので、あるエリアの情報に絞って、集計する必要があります。エリアをどう定義するのかなどは考える必要があります。

1つ目は、Indexingです。位置情報として国単位・県単位・市単位などある一定のエリアを決めてエリアによる集計があります。しかし、大きい場所はデータが大きくなるなどばらつきがでてしまい、良いクラスタとは言えません。無いよりはマシですが。

Girdのように、全体を張り巡らせるようなIndex化するクラスタがとても容易です。 Indexには、geohash や UberのH3 があります。

2つ目は、位置情報以外のデータから、アルゴリズムを通じてクラスタリングをする方法です。DBScan や k-means によるクラスタリングが有名です。

以下の3つを記載してます。

DBScan
k-means
H3

対象データ

利用するデータは、BigqueryのPublicデータを利用します。データロードなど不要で、利用するだけという神ですね。これだけでBigQuery最高。ニューヨークのタクシーデータを使います。ゾーンを定義されたデータセットもありますが、今回はそのような定義が無いデータをクラスタ化する為、使用しません。

利用データ: bigquery-public-data.new_york_taxi_trips.tlc_yellow_trips_2016
全体データ量: およそ 18 GB
件数: 131,165,043

少し多すぎるので、5月のみに絞ります。なぜか7月以降からPickupした位置情報がありませんでした。マスク化されたのでしょうかね。

SELECT count(1)  FROM `bigquery-public-data.new_york_taxi_trips.tlc_yellow_trips_2016` 
WHERE '2016-05-01'<= pickup_datetime and pickup_datetime  <  '2016-06-01' and pickup_longitude is not null and not pickup_longitude = 0
LIMIT 1000

11,689,146

お題

Pickupしたエリアの情報元に、支払われた平均額を求めます。これによって、儲けやすいエリアを割り出せるので、タクシー運転手は嬉しい情報ですね。

BigQuery Geo Viz によって確認していきます。

まずは、クラスタ化しない生データのまま表示します。全部表示出来ないので、10万件にサンプリングしています。また、小さい料金が多すぎて面白みがないので、20以上にしています。これだけでも十分に何かは得られるかもしれません。

SELECT total_amount, ST_GEOGPOINT(pickup_longitude, pickup_latitude) as geom FROM `bigquery-public-data.new_york_taxi_trips.tlc_yellow_trips_2016` 
WHERE '2016-05-01'<= pickup_datetime and pickup_datetime  <  '2016-06-01' and pickup_longitude is not null and not pickup_longitude = 0 and total_amount > 20
LIMIT 100000

f:id:aquarickn:20201130094140p:plain — 1

DBScan

DBScanは、既にGIS関数で用意されてます。まさに神。

しかし、1000万件のデータではDBScanは実行出来なかった。そのため、今回は100万件に絞りました。残念です。 100万件から1200クラスタまで減らせました。ただ、DBScanは密集した距離が近いもの全てを巻き込む性質があるので、パラメータによっては密集エリアで稼げる場所が割り出せなくなってしまいます。

create table `aaa.dbscan_new_york_taxi` as
WITH
  raw AS (
  SELECT 
    total_amount,
    ST_GEOGPOINT(pickup_longitude, pickup_latitude) as geom
  FROM
    `bigquery-public-data.new_york_taxi_trips.tlc_yellow_trips_2016`
   WHERE '2016-05-01'<= pickup_datetime and pickup_datetime  <  '2016-06-01' and pickup_longitude is not null and not pickup_longitude = 0
   LIMIT 1000000
  ),
  clusters AS (
  SELECT
    total_amount,
    ST_CLUSTERDBSCAN(geom, 100, 5) OVER() AS cluster_id,
    geom
  FROM raw
 )
SELECT
  cluster_id,
  ST_CENTROID_AGG(geom) AS center_geom,
  MAX(total_amount) AS max_ta,
  MIN(total_amount) AS min_ta,
  AVG(total_amount) AS avg_tda
  COUNT(1) AS cnt
FROM
  clusters
GROUP BY
  cluster_id
ORDER BY
  cluster_id DESC

人が多いほど circleRadiusを大きくしています。

f:id:aquarickn:20201130094315p:plain

k-means

k-meansはBQの関数に存在しません。その代わりBQではMLの機能があるので、組み合わせることでk-meansを実装出来ます。ただし、100クラスタまでしか作れませんでした。

DBScanのように近しい人々のクラスタリングのモデル

geograpy型は対応していないので、FLOATで行う

-- Train
CREATE OR REPLACE MODEL
  aaa.near_people_new_york_taxi_model_k_means
OPTIONS
  (model_type='kmeans',
    num_clusters=100,
    standardize_features = TRUE) AS
SELECT ST_X(ST_GEOGPOINT(pickup_longitude, pickup_latitude)) as lng, ST_Y(ST_GEOGPOINT(pickup_longitude, pickup_latitude)) as lat FROM `bigquery-public-data.new_york_taxi_trips.tlc_yellow_trips_2016`
WHERE '2016-05-01'<= pickup_datetime and pickup_datetime  <  '2016-06-01' and pickup_longitude is not null and not pickup_longitude = 0 LIMIT 1000000

f:id:aquarickn:20201130094316p:plain

評価

f:id:aquarickn:20201130094420p:plain

予測テーブル

-- predict
CREATE TABLE
  `aaa.near_people_new_york_taxi_model_k_means_predict` AS
SELECT
  * EXCEPT(nearest_centroids_distance)
FROM
  ML.PREDICT(MODEL aaa.near_people_new_york_taxi_model_k_means,
    ( SELECT
    total_amount,
    ST_GEOGPOINT(pickup_longitude, pickup_latitude) as geog,
    ST_X(ST_GEOGPOINT(pickup_longitude, pickup_latitude)) as lng,
    ST_Y(ST_GEOGPOINT(pickup_longitude, pickup_latitude)) as lat
  FROM
    `bigquery-public-data.new_york_taxi_trips.tlc_yellow_trips_2016`
  WHERE
    '2016-05-01'<= pickup_datetime
    AND pickup_datetime < '2016-06-01'
    AND pickup_longitude IS NOT NULL
    AND NOT pickup_longitude = 0 ))

表示 100クラスタに分けてみましたが、今回のお題では使い方を誤ったなと。一応クラスタ化されているようです。

f:id:aquarickn:20201130094457p:plain

金額でのクラスタリング

-- Train
CREATE OR REPLACE MODEL
  aaa.total_amount_new_york_taxi_model_k_means
OPTIONS
  (model_type='kmeans',
    num_clusters=5,
    standardize_features = TRUE) AS
SELECT total_amount FROM `bigquery-public-data.new_york_taxi_trips.tlc_yellow_trips_2016`
WHERE '2016-05-01'<= pickup_datetime and pickup_datetime  <  '2016-06-01' and pickup_longitude is not null and not pickup_longitude = 0 LIMIT 1000000

f:id:aquarickn:20201130094515p:plain

評価

f:id:aquarickn:20201130094533p:plain

予測テーブル

-- predict
CREATE TABLE
  `aaa.total_amount_new_york_taxi_model_k_means_predict` AS
SELECT
  * EXCEPT(nearest_centroids_distance)
FROM
  ML.PREDICT(MODEL aaa.near_people_new_york_taxi_model_k_means,
    ( SELECT
    total_amount,
    ST_GEOGPOINT(pickup_longitude, pickup_latitude) as geog,
    ST_X(ST_GEOGPOINT(pickup_longitude, pickup_latitude)) as lng,
    ST_Y(ST_GEOGPOINT(pickup_longitude, pickup_latitude)) as lat
  FROM
    `bigquery-public-data.new_york_taxi_trips.tlc_yellow_trips_2016`
  WHERE
    '2016-05-01'<= pickup_datetime
    AND pickup_datetime < '2016-06-01'
    AND pickup_longitude IS NOT NULL
    AND NOT pickup_longitude = 0 ))

表示1

SELECT CENTROID_ID, total_amount as ta,geog FROM aaa.total_amount_new_york_taxi_model_k_means_predict limit 100000

f:id:aquarickn:20201130094614p:plain

表示2

SELECT CENTROID_ID, avg(total_amount) as ta, ST_CENTROID_AGG(geog) as geog  FROM aaa.total_amount_new_york_taxi_model_k_means_predict 
group by CENTROID_ID

センター位置でまとめてみました。

f:id:aquarickn:20201130094634p:plain

H3

H3を実装していきます。H3はUDFを利用します。ここからGCSへアップロード

2つ H3変換のUDFを作成

-- Geo情報からH3Indexを取得
CREATE FUNCTION `aaa.geoToH3`(lat FLOAT64, lng FLOAT64, res INT64)
RETURNS STRING
LANGUAGE js AS """
  return h3.geoToH3(lat, lng, res);
"""
OPTIONS (
  library="gs://{bucket}/h3-js.umd.js"
);

-- H3Index から H3Polygon のGeo情報を取得
CREATE FUNCTION `aaa.h3ToGeoBoundary`(h3Index STRING)
RETURNS STRING
LANGUAGE js AS """
  var formatAsGeoJson = true;
  var coords = h3.h3ToGeoBoundary(h3Index, formatAsGeoJson);
  return JSON.stringify({"type": "Polygon", "coordinates": [coords]});
"""
OPTIONS (
  library="gs://{bucket}/h3-js.umd.js"
);

加工・集計

--- リアルタイムに計算すると遅いので、一旦保存した
create table `aaa.h3_new_york_taxi` as
WITH
  raw AS (
  SELECT
    ST_GEOGPOINT(pickup_longitude, pickup_latitude) AS geom,
    total_amount
  FROM
    `bigquery-public-data.new_york_taxi_trips.tlc_yellow_trips_2016`
  WHERE '2016-05-01'<= pickup_datetime and pickup_datetime  <  '2016-06-01' and pickup_longitude is not null and not pickup_longitude = 0
)
SELECT
  geom,
  total_amount,
  `aaa.geoToH3`(ST_Y(geom), ST_X(geom),0)  as h3_0,
  `aaa.geoToH3`(ST_Y(geom), ST_X(geom),1)  as h3_1,
  `aaa.geoToH3`(ST_Y(geom), ST_X(geom),2)  as h3_2,
  `aaa.geoToH3`(ST_Y(geom), ST_X(geom),3)  as h3_3,
  `aaa.geoToH3`(ST_Y(geom), ST_X(geom),4)  as h3_4,
  `aaa.geoToH3`(ST_Y(geom), ST_X(geom),5)  as h3_5,
  `aaa.geoToH3`(ST_Y(geom), ST_X(geom),6)  as h3_6,
  `aaa.geoToH3`(ST_Y(geom), ST_X(geom),7)  as h3_7,
  `aaa.geoToH3`(ST_Y(geom), ST_X(geom),8)  as h3_8,
  `aaa.geoToH3`(ST_Y(geom), ST_X(geom),9)  as h3_9,
  `aaa.geoToH3`(ST_Y(geom), ST_X(geom),10) as h3_10,
  `aaa.geoToH3`(ST_Y(geom), ST_X(geom),11) as h3_11,
  `aaa.geoToH3`(ST_Y(geom), ST_X(geom),12) as h3_12,
  `aaa.geoToH3`(ST_Y(geom), ST_X(geom),13) as h3_13,
  `aaa.geoToH3`(ST_Y(geom), ST_X(geom),14) as h3_14,
  `aaa.geoToH3`(ST_Y(geom), ST_X(geom),15) as h3_15,
FROM
  raw

表示

with raw as (
SELECT 
  -- 0 ~ 15 resolution
    h3_8 as h3,
    max(total_amount) as max_ta,
    min(total_amount) as min_ta,
    avg(total_amount) as avg_ta,
    count(1) as cnt_ta
FROM `aaa.h3_new_york_taxi`
GROUP BY 
  h3
)
select 
  h3,  
  ST_GEOGFROMGEOJSON(aaa.h3ToGeoBoundary(h3)) as h3_geom,
  max_ta,
  min_ta,
  avg_ta,
  cnt_ta
 from raw

100万件のデータが、2800件だけになりました。サンプリング無しで全てを集計してます。

fillColor は、avg_ta で平均の額を元に色濃くしています。 strokeColor は、cnt_ta で乗車人数が多いほど、色濃くしています。

f:id:aquarickn:20201130094730p:plain

全データで表示できたので、Resolution 4くらいだと、位置情報がバグってる？ようなものもあって面白い。海の上・・。

f:id:aquarickn:20201130094801p:plain

まとめ

このデータであれば、タクシーの運ちゃん大儲け！という訳までには行かなそうです。目安ですね。

DBScanやk-meansによるクラスタリングは、今回のように結果をマップに表示する方式が合いませんでした。データ量の削減もポイントをどの位置にするかを考える必要があります。Centerなら簡単だけど、数が多すぎて微妙。

今回の例だと、Indexing形式に変換することで、クラスタ化するのが一番効果的でした。表示のデータ容量も多くが削減でき、このままMapBoxなどを利用してWebサービスへ出すことも可能です。地図サービスを作っていく上で参考になるかもしれません。

DailyDiaryCron

かきため

BigqueryGISと機械学習やってみた

GIS