LuigiとBQで計算するVariogram

# 概要

これは [Python, LuigiでPipeline管理の基本を学ぶ](https://qiita.com/yuusei/items/6ba669a781b9f8ec7f63)の続編です。

今回はpandas gbqを使って[(参考)](https://qiita.com/i_am_miko/items/68cb516ad2be61d59554)、bigqueryを操作するpipelineを作成しようと思います。

ただ、単純に操作するだけでは物足りないので、最近勉強している空間統計学に絡めて[Variogram](https://en.wikipedia.org/wiki/Variogram)を計算する事を考えます。

### 環境
macOS Mojave (バージョン 10.14.6)
Python バージョン 3.7.4
pandas-gbq 0.13.2

### やる事

[Variogram](https://en.wikipedia.org/wiki/Variogram)とは、空間統計学における「距離と観測値の相関」です。
(詳しくは[GDS:バリオグラムとは](http://www.math.keio.ac.jp/~kei/GDS/2nd/vario.html)や[空間統計学への招待：3.1 バリオグラム(Variogram)の推定](https://www.msi.co.jp/splus/learning/spatial/3-1.html)等を見てください)

今回は[前回](https://qiita.com/yuusei/items/6ba669a781b9f8ec7f63)で学んだ内容を踏まえ、仕事で頻繁に使うBigQueryと連携するために[pandas-gbq](https://pandas-gbq.readthedocs.io/en/latest/)を使っていこうと思います。
([BigQueryとPython(pandas)を連携させてみた](https://qiita.com/i_am_miko/items/68cb516ad2be61d59554)を参考にしています)

### PipeLine構造

| Task名       |  概要 |
|:-----------------|:------------------|
| TaskA             |  BQのテーブルからSQLで一次加工したデータをローカルにcsv保存 |
| TaskB           |  TaskAで保存したcsvにテーブル名をつけてBQへアップロード |
| TaskC             | TaskBの結果を元に2地点間のvariogramを計算してcsv保存 |

となります。

### データ

COVID19の感染者を日別に、検査都道府県別に記録した公開資料(7/16まで分)からVariogramを求めるための緯度経度情報等をを取り出して事前にBQに格納しています。

csvでの見た目はこんな感じです。

<img width="691" alt="キャプチャ.PNG" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/138270/bb1f6195-4d55-032d-09c4-bdfd16a9d70d.png">

# コード

それではコードを見ていきましょう

フォルダ構造は作業フォルダ下に出力したcsvを入れるoutputフォルダとTarget格納用のmanagementフォルダがあるだけのシンプル構造です。
各タスクをモジュール化して流用しやすくしました。

```
.
└output
└management
```

### PipeLine

TaskA->TaskB->TaskCの順に実行するだけのシンプルコードです。

```{python}
import os
import luigi

import taskA
import taskB
import taskC

path_manage = './management/'
path_output = './output/'

class TaskA(luigi.Task):
    """
    TaskA class
    """
    def requires(self):
        return []         # In case of no pre-running task, an empty list.

    def output(self):
        return luigi.LocalTarget(path_manage + 'TaskA.txt')  # Where the managiment file is.
        # return luigi.LocalTarget(path_manage + 'TaskA')  # Where the managiment file is.

    def run(self):
        taskA.task_a()

      # If the line does not exist no managed file created
        with self.output().open('w') as out_file:
            out_file.write("")

class TaskB(luigi.Task):
    """
    TaskB class
    """
    def requires(self):
        return [TaskA(), ]  # TaskB depends on TaskA.

    def output(self):
        return luigi.LocalTarget(path_manage + 'TaskB.txt')  # Where the managiment file is.

    def run(self):
        taskB.task_b()

      # If the line does not exist no managed file created
        with self.output().open('w') as out_file:
            out_file.write("")

class TaskC(luigi.Task):
    """
    TaskB class
    """
    def requires(self):
        return [TaskB(), ]  # TaskB depends on TaskA.

    def output(self):
        return luigi.LocalTarget(path_manage + 'TaskC.txt')  # Where the managiment file is.

    def run(self):
        taskC.task_c()

      # If the line does not exist no managed file created
        with self.output().open('w') as out_file:
            out_file.write("")

if __name__ == '__main__':
    luigi.run()
```


### TaskA

簡単な前処理をして結果をcsvに出力します。

```{python}
import pandas as pd
import itertools

project_id = '<your_pjt_id>'
dataset_id = '<your_dataset_id>'
table_id = '<your_table_id>'

start_dt = '2020-02-01'
end_dt = '2020-08-01'

def task_a():

  query = '''

  #standardSQL

  WITH tmp AS  -- 簡単な条件抽出
   (
     SELECT Hospital_Pref
               , lon
               , lat
               , yyyymmdd
        FROM {dataset}.{table}
     WHERE Hospital_Pref != 'Unknown'
          AND ( lon > 0 OR lat > 0 )
  )

  SELECT Hospital_Pref  -- 抽出用の集計
             , avg(lon) as lon
             , avg(lat) as lat
             , count(*) as CNT
    FROM tmp
   WHERE   '{start_date}' <= yyyymmdd and yyyymmdd < '{end_date}'
  GROUP BY Hospital_Pref
  '''

     
  # strベースのクエリと、project_idが必要
  df = pd.read_gbq(query.format(dataset=dataset_id
    , table=table_id
    , start_date = start_dt
    , end_date = end_dt
    )
  , project_id)

  df.to_csv('./output/{p0}to{p1}.csv'.format(p0=start_dt, p1=end_dt), index=False)

# main関数
if __name__ == "__main__":
  task_a()
```


### TaskB

TaskAで出力したcsvをBQにアップロードします。
pandas-gbqを事前にインストールしないといけませんが、
インストールしてあればimportは不要です。

```{python}
import pandas as pd
import itertools

project_id = '<your_pjt_id>'
dataset_id = '<your_dataset_id>'

start_dt = '2020-02-01'
end_dt = '2020-08-01'

# strベースのクエリと、project_idが必要
def task_b():  # 書き込む
  # 読み込む
  df = pd.read_csv('./output/{p0}to{p1}.csv'.format(p0=start_dt, p1=end_dt))

  # BQに書き出す
  df.to_gbq('{dataset_id}.variogram{p0}to{p1}'.format(
  	          dataset_id=dataset_id
  	        , p0=start_dt.replace(u'-',u'')
  	        , p1=end_dt.replace(u'-', u'')
  	          )
            , project_id)

# main関数
if __name__ == "__main__":
  task_b()
```


### TaskC

TaskBで作成した中間テーブルをCROSS JOINして、BQ組み込みの地理計算関数 ST_DISTANCE で距離計算、また、二点間のVariogram(差の二乗平均)を求めています。

```{python}
import pandas as pd
import itertools

project_id = '<your_pjt_id>'
dataset_id = '<your_dataset_id>'
table_id = 'variogram'
start_dt = '2020-02-01'.replace(u'-',u'')
end_dt = '2020-08-01'.replace(u'-',u'')

def task_c():
  query = '''

  #standardSQL
  WITH tmp AS -- CROSS JOIN用1
   (SELECT * FROM {dataset}.{table}{start_date}to{end_date}),
          tmp2 AS   -- CROSS JOIN用2
  (SELECT * FROM {dataset}.{table}{start_date}to{end_date}),

  tmp3 AS(  -- CROSS JOIN で全組み合わせを採用
   SELECT t1.Hospital_Pref AS t1_Pref
              , t1.lon AS t1_lon
              , t1.lat AS t1_lat
               ,t1.CNT AS t1_CNT
               ,t2.Hospital_Pref AS t2_Pref
              , t2.lon AS t2_lon
              , t2.lat AS t2_lat
              ,t2.CNT AS t2_CNT
              , ST_DISTANCE(ST_GEOGPOINT(t1.lon, t1.lat)
                                    , ST_GEOGPOINT(t2.lon, t2.lat)
                                    ) / 1000 AS dist_km  -- 距離をkm変換
             , POWER(t1.CNT - t2.CNT, 2)/2 AS variogram
     FROM tmp AS t1
     CROSS JOIN tmp2 AS t2
   )
   
   SELECT * FROM tmp3
   '''

     
  # strベースのクエリと、project_idが必要

  df = pd.read_gbq(query.format(dataset=dataset_id
    , table=table_id
    , start_date=start_dt
    , end_date=end_dt
    )
  , project_id)

  df.to_csv('./output/variogram_result{p0}to{p1}.csv'.format(
    p0=start_dt, p1=end_dt
    ), index = False)

# main関数
if __name__ == "__main__":
  task_c()
```


### 出力結果

都道府県の緯度経度から求めた相互の距離をkm換算し、
そこからVariogramを求めました。

csvの出力はこんな感じです。

<img width="577" alt="キャプチャ2.PNG" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/138270/3ad9c578-9133-737b-30a1-ef0ba375ddef.png">

一番右がVariogramの値。右から2番目がkm距離になります。

# 次回予告

空間統計学の話としては、上記を応用して距離を離散化し、経験VariogramやKrigingを求めていければなと思っています。

PipeLineの方の話では2点。

まず、TaskBが冗長なので SQL に create_table を入れることを考えています。
ただし、この場合 pd.read_gbq() の出力結果がEmptyになり Alartが表示されるので、そこの処理をどうしようか考え中です。

もう一つは、典型的なLuigiのPipelineでは、Targetファイルを作成するために以下のコードが必須の様なのですが、

```
# If the line does not exist no managed file created
    with self.output().open('w') as out_file:
        out_file.write("")
```

今回みたいにモジュール化した自作関数を使い、PipeLine構造をシンプルにしたい時(まぁ、実務上はだいたいこれ)はこの項目が削れないかと悩ましくなりました。

そうしたらエムスリーさんのBlogで[機械学習で便利なluigiフレームワークの紹介](https://www.m3tech.blog/entry/2018/10/17/105115)の「他のタスクを実行するだけのタスク」によると、このあたりが省略できそうな感じになってます。

次回は上記の内容を踏まえてブラッシュアップした上で、
空間統計学の内容である経験バリオグラムとクリギング[参考](https://qiita.com/Keisuke08aw/items/7cae37b223ee6eaaa498)まで手を出していけたらと思います。

以上