2020-12-06

BigqueryGISと機械学習やってみた

この記事は、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サービスへ出すことも可能です。地図サービスを作っていく上で参考になるかもしれません。

2018-05-21

日付の羅列からSQLで継続日数を割り出してみた

SQL

日の列があれば、継続日数 を動的に割り出したい。データはだいたいデータベース内に入っている為、SQLで出したい。どうやる？を書きました。あんまりどこにも書いてないんだよねぇ。ほぼSQLパズル。

※ Redshiftを対象としています。

元テーブルと参考データ

create

create table continued ( date Date );
insert into continued values ('2017-01-01'),('2017-01-02'),('2017-01-03'),('2017-01-05'), ('2017-01-06'),('2017-01-07'), ('2017-01-09'), ('2017-01-11');

ここから割り出す

    date    
------------
 2017-01-01
 2017-01-02
 2017-01-03
 2017-01-05
 2017-01-06
 2017-01-07
 2017-01-09
 2017-01-11

開始日を割り出す

開始日を割り出すには、対象日の前の日に居なかったら開始とみなせる。その為、開始日対象日をFで継続中をTと出す。

create table continued_tf as 
select a.date, case when ( select 'T' as TF  from continued b where b.date = (a.date - 1)) = 'T' then 'T' else 'F' end  as TF from continued a;

select 内でサブクエリってパワー使いそうであまり使いたくないけど、仕方ないね。この時点で地獄感ある。

    date    | tf 
------------+----
 2017-01-01 | F
 2017-01-02 | T
 2017-01-03 | T
 2017-01-05 | F
 2017-01-06 | T
 2017-01-07 | T
 2017-01-09 | F
 2017-01-11 | F

最も近い開始日を各日付に割り当てる

あんまりこういういSQLは好きじゃないけど仕方ないかな・・。ユーザー毎に作りたかったらユーザー毎のIDを指定すればよろし。

with aaa as (
  select a.date, a.tf,  
  (
    select
      b.date start_date 
    from continued_tf b
    where b.date <= a.date 
      and b.tf = 'F'
    order by b.date desc
    limit 1
  ) start_date
  from continued_tf a
)
select f.start_date, f.date, f.tf, datediff('days', f.start_date, f.date ) as continue 
from aaa f

postgresql の場合はcontinueを以下に書き換え。

( f.date - f.start_date ) as continue

↓ こんな感じで完了

 start_date |    date    | tf | continue 
------------+------------+----+----------
 2017-01-01 | 2017-01-01 | F  |        0
 2017-01-01 | 2017-01-02 | T  |        1
 2017-01-01 | 2017-01-03 | T  |        2
 2017-01-05 | 2017-01-05 | F  |        0
 2017-01-05 | 2017-01-06 | T  |        1
 2017-01-05 | 2017-01-07 | T  |        2
 2017-01-09 | 2017-01-09 | F  |        0
 2017-01-11 | 2017-01-11 | F  |        0

♡SQL

2018-05-15

自営業4年目にして株式会社として起業した

自営業として屋号Labbit を掲げていたが、5月1日にて株式会社Labbitを設立した。

時間が取れないのと顧問になって頂く税理士の事務所紹介と言うこともあり、設立には司法書士の方にご依頼する形で設立をお願いした。

費用的には、株式会社設立で必ず必要となる20万+数万程度。今の時代、ぐぐると「起業 0円」や「会社設立freee」などのサービスをよく目にするが、法務局に行くという時間が勿体ないのと数万程度ならばお願いしたほうがいいなと感じた。

2017-12-08

フリーランスになって3年間で試したこと

この3年間歩んできて、フリーランスの方（特に新たになろうと思う人）にはお得感があってやった方が良いと思った事を3つ書きます。

社会人から個人事業になったら出来る免除
- 国民年金を全額免除出来る！
- 健康保険を減らせる！
  - 任意継続がオススメ。
  - 文芸美術保険組合に入る！
個人事業税支払事業にならない方法
信用について
- カードローンの審査
  - ネット銀行で審査
  - 大手銀行で審査
- 家のローン審査
  - フリーランス1年目なりたて
  - フリーランス2年目
  - フリーランス3年目
  - 金融機関から信用を得る方法まとめ
さいごに

2017-11-19

Scrapy shell で画像を表示するデバッグ方法

scrapy shell を使ってスクレイピングしていると、実際にアクセスした画像URLをコピペして開くのが面倒だったので、コンソールから今どの画像を見ているのかを知りたくなった。

ダウンロードして開く

当たり前かもしれないが、ローカルで画像を開く場合、一度ダウンロードする必要があった。

$ scrapy shell
fetch('http://yahoo.co.jp')
import requests
import io
from PIL import Image
Image.open(io.BytesIO(requests.get(response.css('img::attr(src)')[0].extract()).content)).show()

Webブラウザから開く

こっちの方が簡単だった。

$ scrapy shell
fetch('http://yahoo.co.jp')
import webbrowser
webbrowser.open(response.css('img::attr(src)')[0].extract())

スクレイピング♡Python

2017-11-09

Terminal.app から Vim でマウススクロールが固まる件について解決した

解決策

set visualbell の設定を削除する。

たったこれだけ

経緯と雑感

Terminal.app を利用していて、マウススクロールをした際、行末および行初まで思いっきり移動した際、固まる問題が解決した。僕はなんちゃって自称Vim使いであり、 .vimrc にそこまでこだわりがなかった。

今回も新しいMac環境となった際、「vim 設定おすすめ」で一番上に来た記事の vimrc を意味もわからず設定している。同じような人は居ないだろうか？絶対いる。

こいつを反映したせいで、私のMacでは、 ｍouse=a を指定するはめになった。なんちゃって Vimmer な私はこの mouse=a は好きじゃない。

しかし、やらなければ、マウススクロールが勝手に固まるのだ。

たまたま Vagrant を利用した環境を整えていたところ、マウススクロールをしてもデフォルトのvimでは固まらない。しかし使いにくいからいつものように設定を入れる。固まる！「そうか、.vimrc の設定を入れることで止まっていたのか。」とやっと気づいた。

1行ずつ試し、ついに見つけた。ど真ん中においてある set visualbell 。こいつだった。

全ての環境の vimrc へ制裁という名の削除を行った事により、私の vim 環境では mouse=a も同時に葬ることに成功した。繰り返す。

成功した。

2017-11-08

ネットワーク速度をコマンドラインから見る方法

Webブラウザから、「速度」と検索して計測している皆さん、こんにちわ。私はついに脱しました。