カテゴリー: unclassified

BigQueryでTIMESTAMPのカラムに対し日付を指定する際に、自分が今まで誤って行った方法をまとめてみる。

課題設定

startTime というTIMESTAMP型のカラムがあるテーブルを対象として、

startTime
----------
2016-01-01 01:00:00 UTC
2013-04-21 11:33:01 UTC
2019-09-27 10:30:20 UTC
…

startTimeが2016年7月1日のレコードのみを取り出したいと想定する。

SELECT startTime, … 
FROM dataset.table
WHERE [  ?  ]  // 2016-07-01のデータのみ取り出したい場合に、どう書くのが適切なのか

（話を単純にするため、UTCからの変換は考慮に入れる必要がないと仮定する） ものすごく単純な課題にもかかわらず、色々踏み抜いてきたので以下一つずつ書いていく。

失敗例1：

日付である 2016-07-01 を指定するぞ、という意識だけで以下のクエリをまず書いた。

WHERE startTime = date(2016, 7, 1)

※ BigQueryのDATE型での特定の日時の書き方は以下を参照：
BigQueryリファレンス: Date functions 対象のstartTimeがTIMESTAMP型だという認識がそもそも持てていない時の失敗。（ No matching signature for operator = for argument types: TIMESTAMP, DATE. と怒られる）

失敗例2：

比較対象の型に合わせないといけないのだった、DATEからTIMESTAMPに変えよう…となって、次に以下のようなクエリを書いた。

WHERE startTime = timestamp('2016-07-01')

※ BigQueryでのTIMESTAMP型での特定の日時の書き方は以下を参照：
BigQueryリファレンス: Timestamp functions これは実行できてしまうのだけど、意図した形では動かない。 timestamp('2016-07-01') は 2016-07-01 00:00:00 UTCのことなので、00時ちょうどのデータしか合致しない。結果が1行も返ってこないか、想定よりとても少ない件数が返ってくることになる（count などで気づければよいのだけど、処理途中の部分だったりすると気づけないことがある）。

失敗例3：

では時刻まで含めた形で指定すれば良いのでは、となって以下のクエリを書いた。

WHERE startTime BETWEEN timestamp('2016-07-01 00:00:00') AND timestamp('2016-07-01 23:59:59')

2016/07/01 のレコードの大半がこれで取得できるけど、正確ではない（終了側の指定の仕方が適切ではない）。 TIMESTAMP の説明として以下の記載がある。
BigQueryリファレンス: Data-types

A TIMESTAMP object represents an absolute point in time, independent of any time zone or convention such as Daylight Savings Time with microsecond precision.

TIMESTAMPはマイクロ秒の精度であることがわかる。範囲としては 0001-01-01 00:00:00 to 9999-12-31 23:59:59.999999 UTCとある。 そのため、記載したクエリだと2016-07-01の23時59分59.5秒のデータなどが範囲に含まれないことになってしまう。

失敗例4：

ではTIMESTAMPの方をDATE型に変換してしまうか、となって以下のクエリを書いた。

WHERE date(startTime) = date(2016, 07, 01)

※ DATE(timestamp_expression)という形式でTIMESTAMPをDATEに変換できる BigQueryリファレンス: Date funcitons これでも2016/07/01のレコードを取得できる。ただ、これがパーティション分割テーブルだったりすると、全てスキャンされることになりパフォーマンスが下がる。 ※ パーティション分割テーブルのスキャン範囲についてはこのあたりに記載がある
BigQueryリファレンス: パーティション分割テーブルに対するクエリ

クエリでスキャンされるパーティションを制限するには、フィルタで定数式を使用します。クエリフィルタで動的式を使用すると、BigQuery はすべてのパーティションをスキャンする必要があります。

成功例：

数々の失敗を経て以下の形式にたどり着いた。

WHERE 
  startTime >= timestamp('2016-07-01 00:00:00') 
  AND startTime < timestamp('2016-07-02 00:00:00')

これだと、意図した範囲を指定できているし、定数式にもなっている。
（この条件を満たす他の書き方もあるのかもしれないけど…） 1つのネタでよくこれだけ失敗できたなという感があるけど、色々回り道したことで理解が深まったのでよかったと思うことにする。

前の記事では、Cloud Functionsを使い、Cloud StorageにCSVファイルが追加されたらBigQueryに自動で読み込むことをしてみた。この際、CSVは予め手動で整形した上でCloud Storageに置いていた。
この記事では、CSVをBigQueryに読み込み可能な状態に整形する過程を自動化する方法を探る。

概要

今回、Cloud Data Fusionを用いてCSV整形を実施してみた。
Cloud Data Fusionは、データパイプラインをGUIで作成できるサービスで、ETLのTransform部分も様々な機能が用意されているため、CSVの整形もGUIで行えてしまう。
今回は、以下の2段階で進めた手順をまとめた。
A. 定型的にCSVを整形し、BigQueryに読み込む手順をCloud Data Fusionで作成する
B. 任意のCSVファイル名を渡したらA.が実行されるようにする

事前準備

クイックスタート|Cloud Data Fusion ドキュメント を参照して、APIを有効にし、インスタンスを作成しておく。
下記のように、インスタンス名の前に緑のチェックマークが入ったら作成完了。

操作方法や概観をつかむため、事前に以下のチュートリアルをやった。
・ターゲティングキャンペーンパイプライン | Cloud Data Fusion ドキュメント 顧客一覧から、住所が特定の条件に合致する顧客情報のみを抜き出すという内容。今回の目的のうち、CSVの整形を行う部分で参考になった。
・再利用可能なパイプラインの作成 | Cloud Data Fusion ドキュメント 引数で情報を渡すよう設定し、再利用可能なパイプラインを構築するという内容。今回の目的のうち、CSVのファイル名を与えたら実行されるよう設定しておく部分で参考になった。

手順

A. 定型的にCSVを整形し、BigQueryに読み込む手順をCloud Data Fusionで作成する

スタートとゴール

前の記事同様、気象庁のデータを使用する。
気象庁のHPからダウンロードした段階では、以下のような形式になっている。このCSVをCloud Storageの特定のフォルダに置いた状態からスタートする。

ダウンロードした時刻：2021/11/21 16:31:04

,東京,東京,東京,東京,東京,東京
年月日,平均気温(℃),最高気温(℃),最低気温(℃),最大風速(m/s),最大風速(m/s),天気概況(昼：06時～18時)
,,,,,風向,
2021/1/1,4.4,10.5,-1.3,3.1,北北東,快晴
2021/1/2,4.8,10.8,0.1,4.7,北北東,快晴

これを下記の形式にData Fusion内で整形し、BigQueryに格納することを目指す。

date,ave_temp,max_temp,min_temp,max_wind_speed,wind_direction,weather
2021-01-01,4.4,10.5,-1.3,3.1,北北東,快晴
2021-01-02,4.8,10.8,0.1,4.7,北北東,快晴

① 読み込み先のBigQueryの準備をしておく

BigQueryに読み込み先のデータセットとテーブルを準備しておく。 スキーマはこれからCloud Data Fusionで整形する予定の内容に合わせておく。

② Cloud Data Fusion UIでCSVを読み込む

「事前準備」でインスタンスを作成してあるので、アクションの「インスタンスを表示」をクリックすると、UIに画面が遷移する。
Wrangleをクリックする。

Upload元としてCloud Storageを選択し、用意しておいたCSVを読み込むと、以下のような感じで表示される。行によってカンマの数が異なるデータのため、うまくカンマ区切りで読み込まれず、全て1列に入る形になっている。

③ データを整形する

ゴールの内容に合うようにデータを整形していく。

• 不要なヘッダーの削除
  今回は冒頭の5行を消してしまい、列名は後から付け直す方針で進める。 Filter > value matches regex で正規表現により日付から始まる列のみ残すように（冒頭の5行を削除するように）した。

• カンマ区切りで区切る Parse > CSV を選択し、区切り文字はCommaを選択する。

ここまでの処理でこのような形式になっている。

• 不要な列を削除する
  カラム区切りで分割した際、分割前の状態も列として残っていている。不要なため Delete columnより削除する。

• 列名をつける
  各列打ち替えていく。

• 型変換
  • 日付：読み込んだ時点では1行全体で1つの文字列として判断されているため、現状日付も文字列となっている。Parse > Simple date より、今回2021/1/1形式だったのでCustom formatを選択し、yyyy/m/dで実行した。

すると、不要な時間までついてきてしまった。

公式のドキュメントでは見つけられなかったのだけど、stack overflowを参照して、Custom transform > date.toLocalDate()をかけたら、Date型になってくれた。

• 数値：Change data type > Floatを選択して変換する。

Wranglerで行った各変換は画面右で表示されており、不要な処理を削除することが可能。

④パイプライン作成

Create a Pipeline をクリックし、Batch pipelineを選択すると

データ元のGCSと設定したWranglerがパイプラインでつながった状態で表示される。

⑤ データの書き出し先を設定する（BigQuery）

Sink > BigQueryを選択（Sinkはデータの出力先）してBigQueryのノードを追加し、Wranglerと接続する。

ポインタをBigQuryノードの上に置くとPropertiesが表示されるので、クリックする。BigQuery側のデータセットやテーブル等の指定を行う。

⑥ デプロイと実行

Runをクリックすると、処理が進んでいく（けっこう時間がかかった 20分程度）。Statusの部分の表示で処理の段階がわかり、無事成功すると最終的にSucceededになる。

B. 任意のCSVファイル名を渡したらA.が実行されるようにする

A. では、元となるCSVファイル名をベタ打ちで指定していた。これだとファイル名が変わるたびに毎回作り直す必要があるので、任意のファイル名を渡して処理を実行できるようにしたい。
A. の処理をベースにし、一部書き換える。

① データの読み込み先の変更

Cloud Storageからの読み込んでいる最初のノードのPropertiesをクリックし、読み込み先のパスの指定部分末尾のCSVファイル名が入る部分を${csv_file_name}に変更する。
※ 下図のFormat部分をcsvとしているが、ここは今回textにする必要があった（csvだと1行目しか読み込まれなかった。おそらく、今回のCSVファイルがCSVと判断される形式でなかったのが要因）

② デプロイ

変数を入れてRunしてみると、無事成功した。

未解決な部分

今回文字列として読み込まれた小数点の数値列について、文字列からFloatに変換すると元々存在しなかった不要な桁が出てしまったが、これをうまく処理することができなかった。

round関数があるので、wranglerでround(ave_temp, 1)などと試してみたが、エラーになり、roundの桁数を指定する方法が見つけられなかった。
この記事の中で、

We could write the recipe for the transformation directly with Directives using JEXL syntax(https://commons.apache.org/proper/commons-jexl/reference/syntax.html)

とあったので、該当ページを参照したが、今回行いたい操作のヒントは得ることができなかった。

まとめ

今回、Data Fusionを利用してCSVファイルの整形を行ってみた。GUIで様々な処理が行えるのはとっつきやすくはあったが、少し込み入った処理をしようとするとやりづらかったり、結局学習コストがかかるので、コードを書いて管理する方が楽かもしれないと感じた。
次は、ワークフローを管理できるCloud Composerを触ってみたい。

前記事で示したOpenCVを用いた矩形検出の改善案のうち、ここでは下記2案の内容について記載する。
・改善案3 エッジ検出を使う
・改善案4 適応的閾値処理を使う

環境

Python 3.7.8
OpenCV-Python 4.5.1.48

改善案3 エッジ検出を使う

案の概要

画像の輝度が急激に変化している箇所をエッジ（今回だと輪郭）として捉える処理を行い、エッジのみになった画像に対して輪郭検知を行うという案。
エッジとみなすのは周辺と比べ値が急激に変化している箇所なので、微分した値が大きい箇所をエッジとみなす形で処理が行われる。調べると、エッジ検出としてはSobelフィルタ、Laplacianフィルタ、Cannyフィルタなどがよく使われているようだったが、ここではCannyフィルタを利用した場合の結果を示す。

案の実践

エッジ検出した結果に対して輪郭の検知を行った。

• 

• 

一見うまくいくように見えるが、実際にはその後矩形のみに絞る処理がうまくいかず、1枚もレシートを検出できなかった。これは、レシートの外形を輪郭が囲っているように見えるがそれらはひとつながりの線になっておらず、別々の矩形と捉えられているからと思われる。
そこで輪郭の検知前にノイズ処理（ここではモルフォロジー変換）を追加し、検出したエッジを単純化（膨張処理）することにした。

• 

• 

概ねレシートの輪郭を捉えられているように見えるが、一部途切れているなどしていて、その後レシートの矩形のみに絞る処理をすると1枚のレシートしか検出できなかった。
そこで、エッジの膨張させる程度をもう少し強くしてみる。

• 

• 

これだと3枚ともレシートを検出できた。
検出結果を見ると、レシートを囲む2重線のうち、内側のみ矩形と判断されていた。これは、外側のレシートを囲む線は背景のノイズを拾った線とつながるなどしていて、矩形と判定されなかったことによる。
一方、内側のレシートを囲む線はレシートの印字の輪郭線とつながっていないため、矩形と判断できた。今回の膨張の程度だと問題ないが、もう少し大きくすると今度はレシートの印字の検知の範囲とレシートを囲む内側の線がくっついてしまい、うまく矩形検出されなくなることが予想される。

※ モルフォロジー変換を膨張でなくクロージング（エッジをいったん膨張させて他の細かいエッジと結合させた後、膨張を戻す）にすれば上記の問題が解決できるのではと考え試してみた。しかし、以下の結果となり、最終的にレシートの矩形はうまくいかなかった。

• 

• 

案の評価

レシートの矩形検出自体は行えるが、上記のようにモルフォロジー変換の強さを各画像に合わせ調整する必要がありそうなので、任意の画像に対応するのは難しいだろう。

コード

案3の一連の処理を行う際に書いたコードを以下に示す。

※ 案3に特有の部分
エッジ処理（Canny法）は cv2.Canny で行っている。

Pythonbinary_img = cv2.Canny(gray_img, 100, 200) # 100はminVal、200はmaxVal
# 画素値の微分値が maxVal 以上であればエッジとみなす  
# 画素値の微分値が minVal 以下であればエッジではないとみなし除外する
# 画素値の微分値が二つの閾値の間の場合、エッジと区別された画素（maxVal以上）につながっていればエッジとみなし，そうでなければエッジではないとみなし除外する

モルフォロジー変換（膨張）は cv2.dilate で行っている。

kernel = np.ones((30,30), np.uint8) # 処理の際参照する領域のサイズ
dilation = cv2.dilate(img, kernel, iterations=1) # iterationsは処理回数

一連の処理を行い、検出結果を画像で返すところまでのコード

import cv2
import numpy as np
from matplotlib import pyplot as plt


def binarize(img):
    """画像を2値化する
    """
    gray_img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    binary_img = cv2.Canny(gray_img, 100, 200)
    plot_img(binary_img, 'binary_img')
    binary_img = binary_img.astype('uint8')
    return binary_img


def noise_reduction(img):
    """ノイズ処理(膨張)を行う
    """
    kernel = np.ones((30,30), np.uint8)
    dilation = cv2.dilate(img, kernel, iterations=1)
    plot_img(dilation, 'dilation')
    return dilation


def find_contours(img):
    """輪郭の一覧を得る
    """
    contours, _ = cv2.findContours(img, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    return contours


def approximate_contours(img, contours):
    """輪郭を条件で絞り込んで矩形のみにする
    """
    height, width, _ = img.shape
    img_size = height * width
    approx_contours = []
    for i, cnt in enumerate(contours):
        arclen = cv2.arcLength(cnt, True)
        area = cv2.contourArea(cnt)
        if arclen != 0 and img_size*0.02 < area < img_size*0.9:
            approx_contour = cv2.approxPolyDP(cnt, epsilon=0.01*arclen, closed=True)
            if len(approx_contour) == 4:
                approx_contours.append(approx_contour)
    return approx_contours


def draw_contours(img, contours, file_name):
    """輪郭を画像に書き込む
    """
    draw_contours_file = cv2.drawContours(img.copy(), contours, -1, (0, 0, 255, 255), 10)
    plot_img(draw_contours_file, file_name)


def plot_img(img, file_name):
    """画像の書き出し
    """
    plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
    plt.title(file_name)
    plt.show()
    cv2.imwrite('./{}.png'.format(file_name), img)


def get_receipt_contours(img):
    """矩形検出までの一連の処理を行う
    """
    binary_img = binarize(img)
    noise_reduction_binary_img = noise_reduction(binary_img)
    contours = find_contours(noise_reduction_binary_img)
    approx_contours = approximate_contours(img, contours)
    draw_contours(img, contours, 'draw_all_contours')
    draw_contours(img, approx_contours, 'draw_rectangle_contours')


input_file = cv2.imread('/path/to/example.jpg')
get_receipt_contours(input_file)

改善案4 適応的閾値処理を使う

案の概要

前記事で示した現在の検出手順では、画像全体の画素値を対象として計算し2値化の閾値を決めていた。適応的閾値処理は画像全体ではなく、画像中の小領域ごとに閾値を計算する方法になる。そのため、領域ごとに光源環境が変わるような画像であっても限られた領域内の画素を対象とすることで、画像全体を対象とした場合よりも良い結果が得られる。

案の実践

OpenCVで用意されている適応的閾値処理の関数として cv2.adaptiveThreshold がある。以下の2引数を動かして変化を確認した。
・Block Size： 閾値計算時に対象にする小領域の大きさ（奇数とする）
・C： 計算された閾値から引く定数

Block Sizeを動かす（Cは2に固定）

• 

• 

• 

Cを動かす（Block Sizeは101に固定）

• 

• 

• 

上記結果をみて、Bloce Size=255、C＝2でやってみた。

• 

• 

• 

レシート3枚拾えているが、余計な背景も矩形として拾ってしまっている。
輪郭検知の段階でかなりノイズが多いことが要因と考え、ノイズ処理（中央値フィルタ）を追加してみた。

• 

• 

• 

今度はレシートのみ過不足なく拾えている。

案の評価

上記のように結果を見ながら閾値を調整したので、画像によってはレシートの矩形検出がうまく行えない恐れがある。

コード

案4の一連の処理を行う際に書いたコードを以下に示す。

※ 案4に特有の部分
適応的閾値処理の関数として cv2.adaptiveThreshold を利用している。

binary_img = cv2.adaptiveThreshold(gray_img, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY, 255, 2)
# 第2引数：輝度値の最大値（今回255）
# 第3引数：閾値計算の方法（今回のADAPTIVE_THRESHOLD_GAUSSIAN_Cだと小領域で閾値を計算する方法にガウス分布による重み付けをした平均値を使うことになる）
# 第4引数：閾値処理の種類（今回のTHRESH_BINARYだと閾値より小さい範囲は黒大きい範囲は白に変換する）
# 第5引数：閾値計算時に対象にする小領域の大きさ（今回255としている）
# 第6引数：計算された閾値から引く定数

ノイズ処理（中央値フィルタ）として cv2.medianBlur を利用している。

median = cv2.medianBlur(img, 9) # 9はカーネルサイズ（中央値を計算する対象とする範囲）

一連の処理を行い、検出結果を画像で返すところまでのコード

import cv2
import numpy as np
from matplotlib import pyplot as plt


def binarize(img):
    """画像を2値化する
    """
    gray_img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    binary_img = cv2.adaptiveThreshold(gray_img, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY, 255, 2)
    plot_img(binary_img, 'binary_img')
    return binary_img


def noise_reduction(img):
    """ノイズ処理(中央値フィルタ)を行う
    """
    median = cv2.medianBlur(img, 9)
    plot_img(median, 'median')
    return median


def find_contours(img):
    """輪郭の一覧を得る
    """
    contours, _ = cv2.findContours(img, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    return contours


def approximate_contours(img, contours):
    """輪郭を条件で絞り込んで矩形のみにする
    """
    height, width, _ = img.shape
    img_size = height * width
    approx_contours = []
    for i, cnt in enumerate(contours):
        arclen = cv2.arcLength(cnt, True)
        area = cv2.contourArea(cnt)
        if arclen != 0 and img_size*0.02 < area < img_size*0.9:
            approx_contour = cv2.approxPolyDP(cnt, epsilon=0.05*arclen, closed=True)
            if len(approx_contour) == 4:
                approx_contours.append(approx_contour)
    return approx_contours


def draw_contours(img, contours, file_name):
    """輪郭を画像に書き込む
    """
    draw_contours_file = cv2.drawContours(img.copy(), contours, -1, (0, 0, 255, 255), 10)
    plot_img(draw_contours_file, file_name)


def plot_img(img, file_name):
    """画像の書き出し
    """
    plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
    plt.title(file_name)
    plt.show()
    cv2.imwrite('./{}.png'.format(file_name), img)


def get_receipt_contours(img):
    """矩形検出までの一連の処理を行う
    """
    binary_img = binarize(img)
    binary_img = noise_reduction(binary_img)
    contours = find_contours(binary_img)
    approx_contours = approximate_contours(img, contours)
    draw_contours(img, contours, 'draw_all_contours')
    draw_contours(img, approx_contours, 'draw_rectangle_contours')


input_file = cv2.imread('/path/to/example.jpg')
get_receipt_contours(input_file)

前記事で示したOpenCVを用いた矩形検出の改善案のうち、ここでは下記2案の内容について記載する。レシートを白・背景を黒に2値化できるような変換を目指す。
・改善案1 減色する
・改善案2 HSV色空間にする

環境

Python 3.7.8
OpenCV-Python 4.5.1.48

改善案1 減色する

案の概要

前記事で示した現在の検出手順では、2値化に大津の2値化を利用している。大津の2値化は画素値の分散を用いて閾値が決められるため、画素値をヒストグラムで表した時に双峰性を持つような分布になる画像だとうまく閾値を決めることができる（と理解している）。
今回の画像のヒストグラムを見てみると、山がいくつもある。また、大津の2値化による閾値は107だったが、背景で白飛びしている箇所（画像右下）が閾値より右（白色側）に入っているためうまくレシートと背景を2値化で分離できていない状態になっている。

単純に考えると、2値化する際に背景の机が黒と判定されれば（ヒストグラムで背景が白色側の山に入らず、黒色側の山に入るように調整できれば）うまくいくはず。
単純化すれば扱う画像は以下の3色で構成されている。
・レシート：白
・背景：任意の単一色（グレースケールだと灰）
・レシートの印字：黒
この3色に減色できている画像をグレースケールに変換すれば、ヒストグラムの山が3つになる。その状態で2値化し、レシート（白）と背景（灰）の間で閾値が引かれればうまくいくという案となる（仮に期待通りに減色できたとしても、背景（灰）と印字（黒）の間で閾値が引かれてしまうという問題は残っているが、背景が1色になれば少なくとも背景が白黒両方に分布することはなくなるはず）。

案の実践

減色はk-meansを用いて行うことができる。（考えてみれば当然だけど、色は3つの数字の組で表されているので、3次元空間でクラスタリングするのと同じことだった）
概要に記載した3色に減色を試してみると以下のようになった。

• 

• 

• 

そう思い通りにはいかず、背景が複数色に分かれ、かつ背景の一部はレシートと同色になる結果だった。
減色後のヒストグラムを確認すると、背景は1色にならず3色いずれにも分布している状態だった。

レシートと背景が同化しないよう、色数（クラスタリング数）を少し増やして5色にすると以下の結果だった。

• 

• 

• 

5色に増やすと、3色の時にはレシートと同色に分類された画像右下部分に関してもレシートと別の色に分けられている。しかし、閾値が期待した位置で引けていないため、結局2値化した際にレシートと背景の一部が同化している。 そこで、2値化の閾値を人間が与える形に変えてみる。何色に減色したかによるが、x色に減色したうちレシートは白側上位1色か2色に属することが多いだろうから、2値化の閾値を白側から2色と3色の間の位置とするルールにしてみた。
※ 以下は7色に減色し、閾値は白側から2色と3色の間に引いた場合

• 

• 

• 

これだと、一応レシートが全て無事検出できた。

案の評価

「レシートは減色したx色のうち白側上位1色か2色になる」という仮定の元、閾値を任意で定めている。そのため、レシートが白側上位3色になっている場合、あるいは白側上位1、2位に背景も含まれてしまっている場合、この案は役に立たなくなる。特に背景が白色系だと背景とレシートをうまく分離できないだろう。

コード

案1の一連の処理を行う際に書いたコードを以下に示す。

※ 案1に特有の部分
　減色は cv2.kmeans を利用している。

pixels = img.reshape(-1, 3).astype(np.float32) # 画像の変換（np.float32型で渡す必要がある）
criteria = cv2.TERM_CRITERIA_MAX_ITER + cv2.TERM_CRITERIA_EPS, 10, 1.0 # 繰り返しの終了条件 これだと精度が1に達するor繰り返し10回いずれかに達したら終了する
attempts = 10 # k-meansの初期値の試行回数
flags = cv2.KMEANS_RANDOM_CENTERS # k-meansの重心の初期値の決め方
_, labels, centers = cv2.kmeans(pixels, K, None, criteria, attempts, flags)

一連の処理を行い、検出結果を画像で返すところまでのコード

import cv2
import numpy as np
from matplotlib import pyplot as plt


def sub_color(img, K):
    """色数を指定して減色する
    """
    pixels = img.reshape(-1, 3).astype(np.float32)
    criteria = cv2.TERM_CRITERIA_MAX_ITER + cv2.TERM_CRITERIA_EPS, 10, 1.0
    attempts = 10
    flags = cv2.KMEANS_RANDOM_CENTERS
    _, labels, centers = cv2.kmeans(pixels, K, None, criteria, attempts, flags)
    sub_color_img = centers[labels].reshape(img.shape).astype(np.uint8)
    plot_img(sub_color_img, 'sub_color_img')
    return sub_color_img


def plot_histgram(img):
    """画像の画素値の分布をヒストグラムにする
    """
    hist = cv2.calcHist([img], [0], None, [256], [0,256])
    plt.bar([i for i in range(0,256)], hist.ravel())
    plt.show()


def binarize(img):
    """画像を2値化する
    """
    gray_img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    plot_img(gray_img, 'gray_img')
    threshold = np.unique(np.array(gray_img).ravel())[-2] -1  # 白側から2色と3色の間の位置を閾値とする
    _, binary_img = cv2.threshold(gray_img, threshold, 255, cv2.THRESH_BINARY)
    plot_img(binary_img, 'binary_img')
    return gray_img, binary_img


def find_contours(img):
    """輪郭の一覧を得る
    """
    contours, _ = cv2.findContours(img, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    return contours


def approximate_contours(img, contours):
    """輪郭を条件で絞り込んで矩形のみにする
    """
    height, width, _ = img.shape
    img_size = height * width
    approx_contours = []
    for i, cnt in enumerate(contours):
        arclen = cv2.arcLength(cnt, True)
        area = cv2.contourArea(cnt)
        if arclen != 0 and img_size*0.02 < area < img_size*0.9:
            approx_contour = cv2.approxPolyDP(cnt, epsilon=0.01*arclen, closed=True)
            if len(approx_contour) == 4:
                approx_contours.append(approx_contour)
    return approx_contours


def draw_contours(img, contours, file_name):
    """輪郭を画像に書き込む
    """
    draw_contours_file = cv2.drawContours(img.copy(), contours, -1, (0, 0, 255, 255), 10)
    plot_img(draw_contours_file, file_name)


def plot_img(img, file_name):
    """画像の書き出し
    """
    plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
    plt.title(file_name)
    plt.show()
    cv2.imwrite('./{}.png'.format(file_name), img)


def get_receipt_contours(img, K):
    """矩形検出までの一連の処理を引数の色数で行う
    """
    sub_color_img = sub_color(img, K)
    gray_img, binary_img = binarize(sub_color_img)
    contours = find_contours(binary_img)
    approx_contours = approximate_contours(img, contours)
    draw_contours(img, contours, 'draw_all_contours')
    draw_contours(img, approx_contours, 'draw_rectangle_contours')
    plot_histgram(gray_img)


# 7色に減色して矩形検出を試す
input_file = cv2.imread('/path/to/example.jpg')
get_receipt_contours(input_file, 7)

改善案2 HSV色空間にする

案の概要

ほとんどのレシートは白色なので、色を条件として背景と分離できないかと考えた（背景も白色系だと使えなくなるが）。

色を条件とした検出の場合、RGB色空間よりHSV色空間を使ったほうが検出しやすいらしい。
・RGBは、赤(R)・緑(G)・青(B)の各要素がどれだけ含まれているか（3つの色の混色の割合）で表される。同一の色でも明度や彩度の違いによってRGB3つのパラメータが変動し、範囲を指定しづらい。
・HSVは、色相(H)・彩度(S)・明度(V)のパラメータを使って表される。色相（色合い）を単独で指定できる（= 特定の色を指定しやすい）。

白色はHSV色空間だと色相は関係なくなる（0°~360°全て）ので、その長所が生かせないような気もするが、色相：制限なし / 彩度：小さめ / 明度：大きめという範囲を条件として処理してみる。

案の実践

以下の手順で処理を行う。
1. 画像をHSV色空間に変換する
2. レシートの白色とみなす範囲をHSVで指定して、それ以外はマスク（黒に変換）する
3. 2値化できた状態になるので、輪郭の検出を行う

• 

• 

• 

2値化が期待した形でできているので、3枚とも検出できている。

案の評価

任意の画像だと、背景色やレシートの白色度合いによってはレシートの色（白色）とみなす指定範囲を調整する必要がある。今回の画像では偶然うまくいったが、レシートの色（白色）とみなす範囲からレシートの一部が外れた場合、その箇所は黒に分類されるため矩形検出がうまくいかなくなる。

• 

• 

• 

コード

案2の一連の処理を行う際に書いたコードを以下に示す。

※ 案2に特有の部分
HSV色空間への変換後、白色部分のみにするマスク処理は cv2.inRange で行っている。

hsv_img = cv2.cvtColor(img, cv2.COLOR_BGR2HSV) # HSV色空間への変換
lower_white = np.array([0,0,100]) # 白色とみなすHSVの各値の下限
upper_white = np.array([180,25,255]) # 白色とみなすHSVの各値の上限
binary_img = cv2.inRange(hsv_img, lower_white, upper_white) # 上限と下限を指定してマスク処理する

一連の処理を行い、検出結果を画像で返すところまでのコード

import cv2
import numpy as np
from matplotlib import pyplot as plt


def binarize(img):
    """画像を2値化する
    """
    hsv_img = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)
    lower_white = np.array([0,0,100])  # 白色とみなすHSVの各値の下限
    upper_white = np.array([180,25,255])  # 白色とみなすHSVの各値の上限
    binary_img = cv2.inRange(hsv_img, lower_white, upper_white)
    plot_img(binary_img, 'binary_img')
    return binary_img


def find_contours(img):
    """輪郭の一覧を得る
    """
    contours, _ = cv2.findContours(img, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    return contours


def approximate_contours(img, contours):
    """輪郭を条件で絞り込んで矩形のみにする
    """
    height, width, _ = img.shape
    img_size = height * width
    approx_contours = []
    for i, cnt in enumerate(contours):
        arclen = cv2.arcLength(cnt, True)
        area = cv2.contourArea(cnt)
        if arclen != 0 and img_size*0.02 < area < img_size*0.9:
            approx_contour = cv2.approxPolyDP(cnt, epsilon=0.01*arclen, closed=True)
            if len(approx_contour) == 4:
                approx_contours.append(approx_contour)
    return approx_contours


def draw_contours(img, contours, file_name):
    """輪郭を画像に書き込む
    """
    draw_contours_file = cv2.drawContours(img.copy(), contours, -1, (0, 0, 255, 255), 10)
    plot_img(draw_contours_file, file_name)


def plot_img(img, file_name):
    """画像の書き出し
    """
    plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
    plt.title(file_name)
    plt.show()
    cv2.imwrite('./{}.png'.format(file_name), img)


def get_receipt_contours(img):
    """矩形検出までの一連の処理を行う
    """
    binary_img = binarize(img)
    contours = find_contours(binary_img)
    approx_contours = approximate_contours(img, contours)
    draw_contours(img, contours, 'draw_all_contours')
    draw_contours(img, approx_contours, 'draw_rectangle_contours')


input_file = cv2.imread('/path/to/example.jpg')
get_receipt_contours(input_file)

以前作ったレシートのOCRアプリを改善するため、OpenCVを利用した矩形検出について試行錯誤を行った。

環境

Python 3.7.8
OpenCV-Python 4.5.1.48

課題設定

レシートOCRアプリにおいて、一部のレシートが正しく検出されないことがそれなりの頻度で発生している（人間の目で見る分には背景との差異がある程度あるように思える場合にも）。レシートの矩形を過不足なく検出できるようにしたい。

現在の検出手順

どこに改善点があるか調べるため、現在の処理手順を順に追ってみる。
[ 手順 ]
1. 読み込んだ画像をグレースケールに変換
2. 白と黒の2値化（大津の2値化を利用）
3. 輪郭の検知
4. 検知した輪郭を条件で取捨（内の面積が一定以上ある輪郭に絞った上で、輪郭の形状を近似し、頂点が4点の輪郭のみ選択）

• 

• 

• 

• 

手順2.の2値化の段階でレシートの地が白、背景とレシート印字が黒になるイメージだったが、背景の一部も白になってしまっている。そのため、レシートと背景の境界が一部消え、手順3.で輪郭が検知されない結果となっている。

改善ポイント

2値化の段階で、レシートと背景をうまく分離できれば、その後の輪郭の検知等はうまくいくと思われる。そのため、レシートを白・背景を黒に変換する精度を上げることを目指し、2値化の前処理について下記のような複数の手法を試してみた。
改善案1 減色する
改善案2 HSV色空間にする
改善案3 エッジ検出を使う
改善案4 適応的閾値処理を使う
改善案5 ハフ変換を使う ※ レシート外形を検知できなかったので省略

先に結果

上記改善案1~4を試した結果を先に示しておく。

詳細は別記事に記載しているが、この画像を元にパラメータを調整したのでどの手法でも全てのレシートが検出できている。
どの程度汎用性があるか見るため、異なる条件を付加した2パターンの画像で試してみた。
まず、背景の明るさの変化が大きい画像で試した。

背景の明るさの変化というより、各手法の中で指定しているパラメータの範囲外になることで検出できない部分が生じている。
次に部分的に影がかかっている画像で試した。

影の影響を受けやすい案1、2において、検出できないレシートが生じている。

各改善案の詳細

各改善案の詳細・コードは下記に分割して記載した。
OpenCVを利用した矩形検出の試行錯誤_減色・色空間の変更
・改善案1 減色する
・改善案2 HSV色空間にする
OpenCVを利用した矩形検出の試行錯誤_エッジ検出・適応的閾値処理
・改善案3 エッジ検出を使う
・改善案4 適応的閾値処理を使う

まとめ

今回試した3つの画像に限ると、適応的閾値処理が一番検出率が高い結果となった。しかし、この方法も万能ではなく、設定したパラメータから外れる画像に対しては検出できない結果となる。
より良い手法や調整の仕方がありそうだが、今の知識でわかるのはここまでなので、いったんレシートOCRアプリを適応的閾値処理を利用する形で修正を行った。

これはGMOペパボ ディレクター Advent Calendar 2020の18日目の記事です。

 

 

自由研究の楽しいところ

テーマの選び方

• • 精神や体調は、何によって決まるのか（2018）
    前職で毎日ぎりぎりの闘いをしていた中で、切実に「風邪をひきたくない」「少しでも心健やかにいたい」と感じたことがきっかけ。活動量や食事内容、睡眠時間等の様々な値を記録し、それと体調等の関係性を検証した。

• •  なるべく日陰を通るには、どのルートを選ぶべきか（2019）
    2018夏、初めて福岡で迎える夏があまりに暑く「もしやこれから毎年この暑さを乗り切らねばならぬのか」と恐れおののいたことがきっかけ。出発点と到着点を指定した時、どのルートを通るとどの程度日陰になるか計算を行った。（なお、2018夏は特異な年で、翌年以降はそれほど暑くなかった）

• • レシートの家計簿入力を少しでも楽したい（2020）
    スプレッドシートに手入力していたが、「外税か内税か軽減税率適用商品か目で見て判断するのほんとめんどい」「同じ品目名何回も入力するの人生の無駄だな」と思ったのがきっかけ。レシートをOCRした後、結果を調整してCSVに吐き出す仕組みを作った。

自由研究の進め方

 

 

 

 

 

自由研究のメリットデメリット

・メリット：心の安定に寄与する
自分が動くことで、今までなかったものが生み出される。たとえ、仕事で価値をあまり生み出せていない時期も、自由研究で成果が出れば、自分の中で多少心の支えになる（業務面の課題は何も解決していないのだが、別方面で成果を得ると脳が錯覚するのである）。
また、業務後に仕事とは異なる内容に頭を使わないといけなくなり、業務のあれこれを考える余裕がなくなるので、結果的に気持ちの切り替えもしやすい。
※ なお、自由研究でも進捗が出なくて二重に苦しくなることもままある

・デメリット：可処分時間の減少
自由研究に取り組むと、思った以上に時間を食うことに気づく。しかし、可処分時間の大半を注がないとなかなか進捗が出ないので、楽しくならないというジレンマがある。遊ぶ時間や他の勉強をする時間に影響するので地味につらい。

結論

 

オンライン教育サービスであるCouseraで、機械学習の基礎的な学習講座として有名なMachine Learning講座を受講した。

前々から気になりつつ、「自分がついていけるレベルだろうか」「英語わかんないしな」（講義には日本語訳がついているけどテストや課題は英語）、「Pythonじゃないし」（プログラミング課題はOctaveで提出する）となかなか踏み出せずにいた。
在宅勤務になり家にいる時間が増えたことも後押しになって、受講することができたので、感想等記録しておく。

扱われる内容

下記のキーワードに関連する内容について、11週にわたり講義が行われる。

Week 1　機械学習の概要
線形単回帰、最小二乗法、最急降下法

Week 2　線形重回帰
特徴量のスケーリング、正則化

Week 3　分類　ロジスティック回帰
過学習、正則化、One-vs-All

Week 4・5　ニューラルネットワーク
隠れ層、論理ゲート、backpropagation、gradient checking、ランダム初期化

Week 6　機械学習の評価
交差検証、high bias、high variance、学習曲線、適合率、再現率、F値

Week 7　サポートベクタマシン(SVM)
マージン、決定境界、カーネル法

Week 8　クラスタリング・主成分分析(PCA)
K平均法、局所最適、エルボー法
次元削減、射影誤差、共分散行列

Week 9　異常検知・レコメンドシステム
正規分布、多変量正規分布
協調フィルタリング、類似度

Week 10　大規模データ
確率的勾配降下法、ミニバッチ勾配降下法、逐次学習、並列化

Week 11　Photo OCR
パイプライン、スライディングウインドウ、データ合成、ceiling analysis

学習の進め方

各週、動画による講義 + テスト + プログラミング課題から構成されている。
私は講義を平日にざっと見て（わからないところがあってもあまり気にせず最後までいったん見る）、休日にもう1回見直しながら話の流れをノートにメモしていった。平日見たときは「何言ってるか全然わからん」と思った箇所も、2回目に休日に見た時は「あれ、なんで理解できなかったんだろう」となることも多かった（何回見直してもやっぱりわからん、となる箇所ももちろんあった）。
講師のAndrew先生の説明がわかりやすかったのはもちろんだけど、ノートにまとめようとする過程で理解が曖昧な箇所をつぶし、流れをしっかり捉えられるのが自分にはよかったのかなと思う。

かかった時間

ネット上の体験談だと1ヶ月未満で完了したという方も割と見かけたけど（みんな天才なのかな？と思った）、私は設定されているのと同じペースで進めていったので、約2ヶ月かかった。週により内容の重さの差が大きくて、3時間で終わった週もあれば15時間以上かかった週もあった。

感想

講座全体を通じて感じたことを記録しておく。

• 勉強を続けて納得感を高めていきたいと思った

とても勉強になったので受講してよかったなと思う。
説明が非常にわかりやすいので、そのアルゴリズムの意図や、数式の意味するところの大枠について理解することができた。全体的に納得しながら進めることができたのだけど、「この数式がなぜこうなるかは、ここでは説明しない」という説明で次に進む部分もしばしばあったので、数式の導出などは少しもやもやが残った。また、課題は一から自分でコードを書くのではなく、重要な数式部分のみ自分で書くという形式だった。そのため、ふんわりとした理解になっている自覚がある。自分でPythonで書き直してみるとか、数式の導出をしてみるとかするべきなのだろう。勉強を続けていって、理解を深めていきたいと思う。

• 英語をもうちょっと読めるようになると色々楽だなと思った

講義には有志の方がつけてくれた訳文がついているけど、テストや課題は英語で書かれている。最初の数週は頑張って単語を調べながら読んでいたけど、途中から面倒になってDeepLに突っ込んで訳文を読むようになった。俄然課題の進みが速くなって、自分の英語力の低さを改めて実感した。DeepLは訳が自然で、訳文だと意味が不明で結局原文を読むみたいなことはほとんど発生しなかった（訳文を読んで「解くのに必要な前提条件が足りない…」となって、原文見たら訳されていない箇所があった、ということは数度あった）。

• ペース配分を自分でしなくて済むのって楽だなと思った

オンライン講座では当たり前なのかもしれないけど、見終わった動画にはチェックマークがついたり、プログラミング課題も途中提出の度にどこまで終わったか明確に表してくれる。「今週中にあと3個動画見て課題を解くから、x時間くらいで終わる」とわかると、残りの時間は別の学習に充てようなど見通しが立てやすい。ペース配分を自分でしなくて済むこと、どこまで進んでいるか示してくれることがこんなに楽なのかというのは意外な発見だった。

今後機械学習の勉強をしていく上で、非常に役に立つ講座だった。継続して学習していきたい。

こちらでふりかえりを書いたWebアプリについて、デプロイでうまくいかなかった部分を記録しておく。

概要

ローカルでは一応意図した形で動いてくれていたDjangoをHerokuにデプロイしたところ、エラーが出てうまく動かなかった。そして、そのエラーを解決することができなかった。
自分の理解では、利用しているライブラリが参照しているファイルをうまく参照させることができなかったのでエラーとなったと認識している。
コードはここにいったん置いた。

詳細

1- ログの確認

Herokuでデプロイして、該当サイトにアクセスしたところ「Application error」という表示になっていた。

その下に表示されていた文に従い、ログを確認した。

$ heroku logs --tail

(前略)
2019-12-30T05:46:04.064263+00:00 app[web.1]: import osmnx as ox
2019-12-30T05:46:04.064272+00:00 app[web.1]: File "/app/.heroku/python/lib/python3.7/site-packages/osmnx/__init__.py", line 9, in
2019-12-30T05:46:04.064274+00:00 app[web.1]: from .core import *
2019-12-30T05:46:04.064276+00:00 app[web.1]: File "/app/.heroku/python/lib/python3.7/site-packages/osmnx/core.py", line 10, in
2019-12-30T05:46:04.064278+00:00 app[web.1]: import geopandas as gpd
2019-12-30T05:46:04.064288+00:00 app[web.1]: File "/app/.heroku/python/lib/python3.7/site-packages/geopandas/__init__.py", line 1, in
2019-12-30T05:46:04.064290+00:00 app[web.1]: from geopandas.geoseries import GeoSeries # noqa
2019-12-30T05:46:04.064292+00:00 app[web.1]: File "/app/.heroku/python/lib/python3.7/site-packages/geopandas/geoseries.py", line 15, in
2019-12-30T05:46:04.064294+00:00 app[web.1]: from geopandas.base import GeoPandasBase, _delegate_property
2019-12-30T05:46:04.064296+00:00 app[web.1]: File "/app/.heroku/python/lib/python3.7/site-packages/geopandas/base.py", line 16, in
2019-12-30T05:46:04.064298+00:00 app[web.1]: from rtree.core import RTreeError
2019-12-30T05:46:04.064300+00:00 app[web.1]: File "/app/.heroku/python/lib/python3.7/site-packages/rtree/__init__.py", line 1, in
2019-12-30T05:46:04.064302+00:00 app[web.1]: from .index import Rtree
2019-12-30T05:46:04.064310+00:00 app[web.1]: File "/app/.heroku/python/lib/python3.7/site-packages/rtree/index.py", line 5, in
2019-12-30T05:46:04.064312+00:00 app[web.1]: from . import core
2019-12-30T05:46:04.064314+00:00 app[web.1]: File "/app/.heroku/python/lib/python3.7/site-packages/rtree/core.py", line 125, in
2019-12-30T05:46:04.064316+00:00 app[web.1]: raise OSError("Could not find libspatialindex_c library file")

計算して求めたルートを地図にプロットするのに、osmnxというライブラリを利用している。そして、osmnxではRtreeというライブラリが使用されている。Rtreeのコードの中で、libspatialindex_c library fileが見つからなかったことによりエラーになっていると理解した。

2- ローカルではエラーが起きなかったのは何が違うのか確認する

上記のエラー文にFile "/app/.heroku/python/lib/python3.7/site-packages/rtree/core.py", line 125 とあり、.venv/lib/python3.7/site-packages/rtree/core.pyの125行目付近を確認した。
（/app/.heroku/~というのがどのように確認できるかわからなかったのでvenvでローカルでインストールしたものと同じだろうと考えて上記を確認した）

elif os.name == 'posix':
    if 'SPATIALINDEX_C_LIBRARY' in os.environ:
        lib_name = os.environ['SPATIALINDEX_C_LIBRARY']
    else:
        lib_name = find_library('spatialindex_c')

    if lib_name is None:
        raise OSError("Could not find libspatialindex_c library file")

    rt = ctypes.CDLL(lib_name)
else:
    raise RTreeError('Unsupported OS "%s"' % os.name)

lib_nameが定義できていなかった（環境変数の中にSPATIALINDEX_C_LIBRARYがなく、またspatialindex_cというライブラリも見つけられなかった）ためにエラーとなっていたことがわかった。

ターミナルでPythonを立ち上げて、ローカル環境だとどうなっているか確認してみた。

>>> import os
>>> os.name
'posix'
>>> os.environ
environ({'TERM_PROGRAM': 'Apple_Terminal',(後略)

>>> from ctypes.util import find_library
>>> find_library('spatialindex_c')
'/usr/local/lib/libspatialindex_c.dylib'

os.environでSPATIALINDEX_C_LIBRARYがなかったことから、ローカル環境ではlib_name = find_library('spatialindex_c')の方が実行されたのだろうこと、またspatialindex_cはローカル環境だと/usr/local/lib/libspatialindex_c.dylibにあるのだということがわかった。

検索してみると、ローカル環境で同じような問題にぶつかっている人がいた。Homebrewでspatialindexをインストールしたら解決したとされており、$ brew listで確認してみるとインストール済みだった。そのため、ローカルでは問題が起きなかったのだろうと思った。

3- どうしたら本番環境でも参照できるようになるか考える

ローカル環境では該当ファイルが存在していたため参照できたが、本番環境には該当のファイルがないためエラーになる。
pipでインストールできるファイルならrequirements.txtに書いておけるけど、そうではなさそうだ。それなら本番環境にもそのファイルを置いて、参照できればとりあえずは解決できるのではと考えた。

/usr/local/lib/libspatialindex_c.dylibをコピーして、本番環境のプロジェクト直下にlibディレクトリを作り、置いてみた。

これが正しい行為なのかわからないけど、環境変数にパスを設定してみた。
$ heroku configで、登録できていることが確認できた。

(前略)
SPATIALINDEX_C_LIBRARY: lib/libspatialindex_c.dylib

この修正をしても同様のエラー画面だったため、再度、$ heroku logs --tailでログを確認した。

(前略)
2019-12-30T06:38:08.950338+00:00 app[web.1]: File "/app/.heroku/python/lib/python3.7/site-packages/rtree/core.py", line 127, in
2019-12-30T06:38:08.950339+00:00 app[web.1]: rt = ctypes.CDLL(lib_name)
2019-12-30T06:38:08.950340+00:00 app[web.1]: File "/app/.heroku/python/lib/python3.7/ctypes/__init__.py", line 364, in __init__
2019-12-30T06:38:08.950342+00:00 app[web.1]: self._handle = _dlopen(self._name, mode)
2019-12-30T06:38:08.950344+00:00 app[web.1]: OSError: lib/libspatialindex_c.dylib: invalid ELF header

表示が前回と変わっている。
最初のエラーログで出ていたraise OSError("Could not find libspatialindex_c library file")はなく、次の行のrt = ctypes.CDLL(lib_name)が表示されているので、libspatialindex_cを見つけること自体はできたと考えてよいのだろうか。

エラーの最後の行で、invalid ELF headerと表示されている。検索すると、違うOSでのビルドだから読み込めない場合に表示されるという記載があった。

「2- ローカルではエラーが起きなかったのは何が違うのか確認する」に記載した.venv/lib/python3.7/site-packages/rtree/core.pyで、もう少し上の方に以下のような記載箇所がある。

if 'SPATIALINDEX_C_LIBRARY' in os.environ:
  lib_path, lib_name = os.path.split(os.environ['SPATIALINDEX_C_LIBRARY'])
  rt = _load_library(lib_name, ctypes.cdll.LoadLibrary, (lib_path,))
else:
  rt = _load_library('spatialindex_c.dll', ctypes.cdll.LoadLibrary)
if not rt:
  raise OSError("could not find or load spatialindex_c.dll")

このspatialindex_c.dllの方が必要なのかなと思い（多分違う）、探したが見つけられなかった。

4- いったん諦める

いったん諦める。もう少し知識がついたらまた考えてみる。

目的

昨年の夏、外出のたびに「暑い…少しでも日陰を歩きたい」と思っていた。出かける前に、道路の方位と太陽の位置から影の出来かたを脳内でなんとなくシミュレートするんだけど、だいたい実際の影とずれていた。
夏を迎える前に、時間帯ごとの影の位置を求めておきたい。何時に出かけるのがベストなのかを知った上で、覚悟を持って日差しの元に出て行きたい。
そう思ったので計算した。

作ったもの

まず、与えた建物と道路のモデルに対し、道路の何割が日陰になるかを求めるコードを書いた。
– 入力：求める地点の緯度経度、日付、建物と道路の平面図、建物の高さ
– 出力：建物と道路の平面図に影を重ねた図、道路の方位角ごとの影の割合を求めたグラフ

次に、この結果を使って経路の何割が日陰になるかを求めるコードを書いた。
– 入力：出発地点と到着地点の緯度経度
– 出力：時間帯ごとの影の割合を求めたグラフ

結果

先に結果を示す。

影の計算は、下図のような単純化した建物と道路のモデルを対象とした。4.0ｍの道路から1.0ｍ後退した位置に、幅・奥行・高さが外法で6,000mmの建物が、1.0mの離れで立ち並んでいる状態とした。

道路の方位を22.5°ずつ変え、日の出〜日の入りまで、15分ごとに道路のうち影のかかる割合を計算した。（道路の方位は北を0°とした。南北に通る道路なら0°、北西-南東に通る道路なら45°となる） 地点は福岡市（ここでは福岡県庁の緯度経度を代表として利用）とした。

まずは、夏至(2019/6/22)の結果を示す。

• 南中時刻には、どの方位でもほとんど影ができていない（太陽高度が約80°とかなり天頂に近い位置に太陽があるので、どの方位であっても影ができにくい）。
• 方位によってかなり影の割合が異なることがわかる。例えば90°（東西の道路）だと、8時過ぎには影がなくなり、16時をすぎるまでそれが続く。一方、0°（南北の道路）だと、10時でもまだ半分以上影ができている。

• 特徴的なのは67.5°と112.5°の時。前者を取り出して見てみる。
  

太陽方位と道路の方位が一致、あるいは180°となる（つまり影のできる方位と道路が平行になる）のが日が出ている間に2回あり、その前後は道路に全く影ができない状態になる。
南中後、他の方位だとどんどん影の割合が増えるのに比べ、16時頃を頂点として影の割合が増えた後は太陽方位が道路の方位に近づくのに比例してまた影が減っていっている。

次に、真夏（2019/8/1）の場合の結果を見てみる。

1ヶ月強で、けっこう変化があることがわかる。
方位角と時間帯によって影の割合が違うので、すごく頑張って毎日方違えとかすれば、日差しを避けて生きていけそう。

方位角と時間帯ごとの影の割合がわかったところで、目的地まで至る経路上の影の割合を時間ごとに求めてみる。
今回は緯度経度の代表として利用した福岡県庁から、天神駅までの経路を対象とした。最短経路の方位と距離を求めて、方位については22.5°ずつの方位で近似した。

この方位ごとの延長距離から、経路上の影の割合を時間帯ごとに求めてみた。

方位角に偏りのある経路を選んだので、一番長い方位角(135°)の形状を概ね踏襲した結果となった。

計算方法

Pythonで画像処理ができるOpenCVと地理情報を得られるosmnxを利用して求めた。
コードはJupyter Notebookの形式でGitHubに置いた。

大まかな計算の流れは以下のとおり
1. モデルとなる画像を読み込み、建物の高さを画像のピクセル数に変換しておく
2. 道路の方位ごとの計算を行うため画像を回転させる
3. 建物の輪郭を検出、間引いてから座標化
4. 緯度経度・日付・時刻から太陽高度・方位角を求める
5. 建物の輪郭の座標ごとに影の位置を求め、それを図形化する
6. （建物と影を重ねた図を書き出し）
7. 影と道路の重なる割合を求め、グラフ化
8. ある経路上における道路の方位とその距離を調べる
9. 方位の割合に応じた時間ごとの影の割合を算出する

反省点/改善点

• モデルは適当なの？
  正直十分ではないと思う。道路からの後退距離や隣の建物との離れなど、もっと色々なパターンを用意したモデルの平均値を取るなどした方がよかったかもしれない。例えば、今回のモデルでは隣の建物との離れの位置を道路の両側で同じにしているけど、これを互い違いにしたモデルで計算したら影の割合が10%程度違う時間帯もあった。

• 建物の高さが均一というのはどうなの？
  道路の幅と建物の高さは、ある程度一定（幅の広い道路沿いほど建物の高さは高い）となる傾向にあるのでは、と思う。そうすると、高層地域と低層地域で時間帯間で比べた影の割合は似た傾向になるのでは…という仮説になった。（経路上に高層地域と低層地域両方がある場合は上記の仮説は成り立たなくなるので、現実的ではない）
  今回、GWの10連休中に作り切る、と決めていたので、かなり乱暴な近似をしているけど、影のでき方の傾向性を見る程度には使えるのでは、と考えている。

• 道路全面に対する影の割合を計算しているけど、真ん中を歩くことは実際にはないよね？
  ある程度幅の広い道路なら、中心でなく端の方を歩くので、道路の端1.0m程度での影の割合を計算する方が実際的だったと思う。ただ、そうすると左右の端ごとの計算が必要になり、結果が複雑になるので今回はやらなかった。

• ある程度幅の広い道路だと、街路樹があるよね？
  街路樹は今回考慮に入れることができなかった。建物は垂直方向の遮蔽物だけど、樹木は水平方向の遮蔽物で、影のできる面積がかなり大きいので、本当は考慮に入れたかった。幅の広い道路に対しては、x[m]おきにy[m]の高さの円錐状の木のモデルが立っている、というようにすれば良いかもしれない。

• 隣の建物にかかった影の形が再現できていない
  隣の建物の壁に影がかかると影の形が変形するけど、それは再現できていない。ただ、影を遮った隣の建物による影もあることを考えると、道路にかかる影の量は結果的に変わらないはずと考えている。

• わざわざ画像として読み込まなくても、もっと簡単にモデルを作れたんじゃない？
  たぶんそうだと思う。当初は実際の地図から建物を読み込んでそれに対する影を計算しようとしていたので、そういう回りくどい形式になってしまった。1回計算するのに20分くらいかかる。

• コードが整理できていない
  そう思う。後日整理する。

その他全般的な反省は別記事に書いた。
→　「『夏には日陰を歩きたい』という欲望を満たすために計算する」のふりかえり

前職にいた時、必要以上に真面目に取り組むことで許されようとする悪癖があった（今もちょっと残っている）。
期限に間に合わせるため、土日や年末年始にも非公式に職場にやってきて働くとか、そういう行動のことです。

この行動は、前職においておおむね批判されなかった。仕事量も人員数も調整がきかない場だったので、便利な存在という側面もあったと思う（そうやって働く人も多かった）¹。
一方で、その姿勢を諌めてくれる先輩もいた。「お前のまじめさは美徳なんかじゃねーからな」と。
だって、やるしかないじゃん、と思っていた。さらには、この日々は何らかの形で報われると勝手に思っていた。

3年目のある日、気づいた。
部長や課長、あるいは組織がこの頑張りに報いてくれるわけでないし、私の人生に責任を持ってくれるわけではない（当然、持つ必要もない）。
頑張りによる成果は彼らを喜ばせ、それは単純に嬉しいけど、それを自分の最終目的にしても、（頑張っている内容からして）自分に何かが残るわけではない。

この人生を生きているのは誰？
私だ。
まじか、そうかと思った。

目前の仕事だけでいっぱいな日々は、その間思考停止できて、実は楽だ。「でも今はこれをしなきゃ」が汎用的な言い訳になる。
でも、それを何年も続けていると、自分の価値は目前の仕事だけだから、その出来に精神状態が100%左右されるようになる。そして、頑張る見返りとして、自分の人生の責任を取るのは自分であるということを放棄したい心持ちになっていった。

（しつこいですが、ここまで前職での話です）

転職して数ヶ月、楽しんで働くってなんだ、とぼんやり考え続けていた²。

まだ考えがまとまっていないけど、一つ考えたことは、真面目であることは、真剣なのと似ているけど違うということだ。真面目さは義務感から発している一方、真剣さはその対象を自分ごととして捉えることから発している。真面目さは深刻さへつながっていく一方、真剣さはユーモアと同居できる。
「どうしたらそれができるって思う？」「今の自分・状況で取りうる、最善の行動ってなんだろうね？」と自分に問いかけてやって、困難³にも、にいっと笑って、楽しんでやっていきたい。今、自分は事業にしっかり貢献できるだけの実力が持てていないけれど、その中でも、考えることはできるはずなのだ。

真面目さで評価されようとするよりも、真剣に取り組んで成果を出す方が厳しい世界なのは気づいている。
無意識に慣れた手法を取ってしまうけど、頑張り方も変えていきたい。1歩進んで、1歩下がって、の繰り返しだろうけど、それでも振り返れば進んでいる。

ホワイトボード化までの思考回路

きっかけは、思いついたことをふせんにちまちま書いて壁に貼っていたことだった。
ふせんは狭い。ぼんやりと考えたこととかは切り捨てて内容をまとめ、収まるように調整してしまう。とりあえず書き出したい。余計なことを考えず、自由に書きたいのだ（幼児のお絵描きが画用紙を飛び出して壁や床に行われる時と同じ思考回路である）。

Amazonに壁に貼れるホワイトボードのシートがあったので、はじめはそれを机の横の壁に設置することを考えた。
でも、側面だと45度回転しないと書けない。新しいものを導入した時にちょっとした手間があると、結局活用しなくなってしまう。がんばらなくても使える状態にする必要がある。
どこならそのままの状態で目に入り、書けるか。

机の上だ。

作業内容

机の上にはキーボードや本を置いていたけど、物があるとインクがつくことを気にして自由に書けない。
そこで棚とキーボードの収納スペースを作って避難させた上で、机にホワイトボードのシートを設置した。

→
元の状態（見た目的にはこっちの方がすっきりしていて好き）
→ 改造後（なお、ホワイトボードに書いた計算は間違えている）

コストとしては、ホワイトボードシート・木材・金具類等で、7,000〜8,000円くらい。ホワイトボードシートだけだったら、2,000円しない。
散発的に進めていたけど、累計で、サイズを決めたり計画するのに3〜4時間、材料をネットやホームセンターで集めるのに5〜6時間（キーボードのスライド棚に使ったスライドレールについて調べるのに手間取った）、作業に4〜5時間くらいかかった。

やってみてわかったこと

• 太いペンで、自由に書けるのは、気持ちいい。いらない紙に書いてその後捨てるのと原理的に変わらないはずなのに、それよりも書くことへのハードルが下がる。アイデアを考えるとか、何となく考えてることを発展させるとかいう用途に特に向いてるなと感じた。
• キーボードは避難させたけど、ノートや本を広げる時は結局ホワイトボードの上になるので、それが不便。ノートを使うときには一旦消す必要があるな、とか、この書き込みはしばらく残しておきたいから端の方に書いておこう、とか本質と違うことを色々考えてしまう。そういう点では、壁設置に軍配が上がるなと思った。
• 机が白いと、結構目が疲れる。天板が白い机って結構あるけど、白色度はあまり高くないのは、そういう理由もあるのかもしれない。
• 線が太いので、文字や図が大きくなる。それが書くことの気持ち良さにつながっている部分もあるんだけど、この広さでも狭いなと感じてしまう。
• 職場の机（白い）にナチュラルにペンで書こうとして、一人ではっとする。

まだ設置して間もないので、しばらく使ってみたら、新たな発見があるかもしれない。

新卒で入った地方公共団体で数年間働いて、今回、GMOペパボ株式会社に仲間に入れてもらった（この経緯もいつかまとめられたらと思っている）。
ここでは、公務員という立場から転職して感じた違いについて書いてみる。

一言でいうと、全体的にワンダーランドだった。
共通点の方がむしろ少なくて、いちいち新鮮で面白かった。
入る前の想像と違うということでなく、頭で理解しているのに体がついていかないというような感覚だった。単に仕事内容や進め方の違いだけでなく、何を大切にするかという根本的なところが違っていたので、そう感じたのだと思う。
たくさんの違いの中でも、象徴的だなと感じたのは以下の2点だった。

• 楽しむということ：
  当初、言われた言葉「楽しんでやりましょう」。
  楽しむって何だと思った。今まで、楽しそうにしていると通報されることのある立場だったから（改めて考えるとすごい）、無意識に自制するようになっていた。
  仕事を楽しむってどういうことかわからなかったけれど、まずは「見当違いかなと思ってもとりあえず手を動かしてやってみる」ことをしようと考えた。
  前職で難しい対応をする時、相手の課題をより良く解決しようという視点でなく、訴えられない等の自分を守る視点になっていることが時々あって、そういう時は本当に楽しくないし、その守るという結果以外、何も新しいものが生まれなかった。¹
  おそらく、楽しむというのはそれと逆で、フラットな状態で物事に接し、前向きに取り組んで、結果としてもっとおもしろいものが生まれる、という状態なのかなと現時点では考えている。

• 組織への愛情：
  最も違いを感じた部分かもしれない。
  新卒生の研修終了の会を見ていて、ペパボへのある種純粋すぎるほどの愛を感じて、すごくまぶしかったし、その一員となれたことを嬉しく思った。
  同時に、前職の若手の仲間の顔を思い出していた。ペパボと同じくらい、前職の皆も住民のために頑張っていて、概ね人間関係も良好なのに、なぜ組織への誇りを持てない雰囲気だったのか ² 、そういう空気しか作れなかったことを後輩に申し訳なく思った。
  ペパボの雰囲気は、各人が醸成している部分もあるし、そうなるように意識的に作り維持されている部分もあるのだろうと思う。私はそれをいいなと思って入ってきたので、どうやって作られたものなのか観察して、前職との違いを生んでいる原因は何か考えたいと思っている。