GIS奮闘記

現役GISエンジニアの技術紹介ブログ。主にPythonを使用。

スポンサーリンク

Python で2点の緯度経度から距離を計測する方法

さて、本日は2点の緯度経度から距離を計測する方法について紹介しようと思います。平面上の単純な2点間の距離計測などは簡単にできるかと思いますが、地球は回転楕円体なため(実は地球は完全な球体というわけではありません)、その丸みを考慮した計算が必要になります。

地球の形については国土地理院さんが公開しているサイトをご参照ください。

距離計算方法

ヒュベニの公式というものがあり、こちらであれば簡易的に距離計測をすることができますので、今回はこちらを使ってみます。 公式については以下をご参照ください。

ヒュベニの公式

2点間の距離計測式 D = SQRT((Ay * M)^2 + (Ax * N * cos(P))^2)

Ax:2点の経度の差
Ay:2点の緯度の差
P:2点の緯度の平均

M = Rx (1 - e^2 )/W3:子午線曲率半径
N = Rx/W:卯酉線曲率半径
W = SQRT(1 - E^2*sin(P)2):子午線・卯酉線曲率半径の分母

e = SQRT((Rx2 - Ry2)/Rx2):第1離心率
Rx:長半径(赤道半径)
Ry:短半径(極半径)
E2 = e2:第2離心率

長半径(赤道半径)、短半径(極半径) という見慣れない単語がありますが、こちらに関しては Wikipediaを参照してください。そういう固定値を使うという認識でまずは大丈夫だと思います。

f:id:sanvarie:20210116195907p:plain

色々調べてみると他には以下のような計算方法があるようです。国土地理院さんが公開している計算方法は複雑すぎて私の頭では理解できないですね・・・

距離計測する区間

「蒲田駅」から「横須賀中央駅」の距離を計測してみようと思います。それぞれの緯度経度は以下です。

  • 蒲田駅
    • 緯度 35.562479
    • 経度 139.716073

f:id:sanvarie:20210116143305p:plain

  • 横須賀中央駅
    • 緯度 35.278699
    • 経度 139.670040

f:id:sanvarie:20210116143410p:plain

国土地理院さんのサイトで確認したら「31.761253km」という計測結果が出ました。今回の結果がこの数値にどれだけ近づけるかのでしょうか。楽しみです。

f:id:sanvarie:20210116151537p:plain

実行環境

Windows 10 64bit
Python 3.6.6

サンプルコード

# -*- coding: utf-8 -*-
import math

pole_radius = 6356752.314245                  # 極半径
equator_radius = 6378137.0                    # 赤道半径
latlon_kamata = (35.562479, 139.716073)       # 蒲田駅の緯度経度
latlon_yokosukachuo = (35.278699, 139.670040) # 横須賀中央駅の緯度経度

def cal_distance():

    # 緯度経度をラジアンに変換
    lat_kamata = math.radians(latlon_kamata[0])
    lon_kamata = math.radians(latlon_kamata[1])
    lat_yokosukachuo = math.radians(latlon_yokosukachuo[0])
    lon_yokosukachuo = math.radians(latlon_yokosukachuo[1])

    lat_difference = lat_kamata - lat_yokosukachuo       # 緯度差
    lon_difference = lon_kamata - lon_yokosukachuo       # 経度差
    lat_average = (lat_kamata + lat_yokosukachuo) / 2    # 平均緯度

    e2 = (math.pow(equator_radius, 2) - math.pow(pole_radius, 2)) \
            / math.pow(equator_radius, 2)  # 第一離心率^2

    w = math.sqrt(1- e2 * math.pow(math.sin(lat_average), 2))

    m = equator_radius * (1 - e2) / math.pow(w, 3) # 子午線曲率半径

    n = equator_radius / w                         # 卯酉線曲半径

    distance = math.sqrt(math.pow(m * lat_difference, 2) \
                   + math.pow(n * lon_difference * math.cos(lat_average), 2)) # 距離計測

    print(distance / 1000)

if __name__ == '__main__':
    cal_distance()

以下のような結果になりました。

f:id:sanvarie:20210116151941p:plain

31.7612541505034km」という計測結果が出ました。国土地理院さんのサイトで計測した「31.761253km」とほぼ同じ値になりました。

さいごに

少し複雑な計算式でしたが、蒲田駅~横須賀中央駅間の距離計測に成功しました。今度は球面三角法など他の計算方法も試してみたいですね。

また、ライブラリなどを使用すればこのような計算をしなくてもいいのですが、今回のエントリーを書くにあたって色々調査をする中で勉強になることがたくさんありました。私は GIS 歴はけっこう長いのですが地理的な知識に関してはあまりなく、GIS に携わるものとして地理の基礎的な勉強をもっとしなければいけないなと実感しました。本日は以上です。

参考サイト

https://www.trail-note.net/tech/calc_distance/
http://hp.vector.co.jp/authors/VA002244/yacht/geo.htm

Python と OpenCV で画像の幾何学変換をしてみよう

さて、本日は画像の幾何学変換について書いてみようと思います。以下エントリーと同様に OpenCV を使ってみようと思います。

www.gis-py.com

幾何学変換とは

幾何学変換と聞くだけで頭が痛くなりそうな気がしますがそんなに難しく考える必要はありません。まずは画像を拡大・縮小、回転、スキュー(画像を傾けるような処理)するもの、くらいに捉えておけば大丈夫かと思います(もちろん専門的に勉強するならこれだと足りないと思いますが)。

今回試すこと

  1. 拡大・縮小
  2. 移動
  3. 回転
  4. アフィン変換
  5. 射影変換

使用する画像

大田区公式PRキャラクターの「はねぴょん」の画像を使用します。

f:id:sanvarie:20210111155935p:plain

実行環境

Windows 10 64bit
Python 3.6.6

サンプルコード

拡大・縮小

はねぴょんを横に1.5倍、縦に3倍拡大させてみます。cv2.resize の引数の fx、fyを設定することで拡大縮小の調整ができます。

# -*- coding: utf-8 -*-
import cv2

def main():
    img = cv2.imread(r"D:\blog\data\opencv\hane.jpg")

    res = cv2.resize(img, None, fx=1.5, fy=3,         # 横に1.5倍、縦に3倍にする
                     interpolation = cv2.INTER_CUBIC) # 補間方法(縮小する場合はcv2.INTER_AREA)

    cv2.imshow("res", res)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

if __name__ == '__main__':
    main()

このようにはねぴょんが拡大されたことがわかります。

f:id:sanvarie:20210113150749p:plain

移動

はねぴょんを移動させます。

# -*- coding: utf-8 -*-
import cv2
import numpy as np

def main():
    img = cv2.imread(r"D:\blog\data\opencv\hane.jpg")

    rows, cols, _ = img.shape

    M = np.float32([[1,0,100],[0,1,50]]) # X移動100,Y移動50
    dst = cv2.warpAffine(img, M, (cols,rows))

    cv2.imshow("dst", dst)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

if __name__ == '__main__':
    main()

このように画像が移動したことがわかります。

f:id:sanvarie:20210113151952p:plain

回転

はねぴょんを180度回転させます。

# -*- coding: utf-8 -*-
import cv2
import numpy as np

def main():
    img = cv2.imread(r"D:\blog\data\opencv\hane.jpg")
    rows, cols, _ = img.shape

    M = cv2.getRotationMatrix2D((cols / 2, rows / 2), 180, 1)
    dst = cv2.warpAffine(img, M, (cols, rows))

    cv2.imshow("dst", dst)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

if __name__ == '__main__':
    main()

このように画像が180度回転したことがわかります。

f:id:sanvarie:20210116090245p:plain

アフィン変換

アフィン変換は線形変換(拡大や縮小,回転など)と平行移動を組み合わせた変換です。アフィン変換を行うには入力画像と出力画像の対応点の座標が少なくとも3組必要です。今回は以下のように対応点を設定してみました。

f:id:sanvarie:20210116092836p:plain

# -*- coding: utf-8 -*-
import cv2
import numpy as np

def main():
    img = cv2.imread(r"D:\blog\data\opencv\hane.jpg")
    rows, cols, _ = img.shape

    pts1 = np.float32([[16,36],[138,36],[16, 130]])
    pts2 = np.float32([[34,67],[154,84],[63,175]])

    M = cv2.getAffineTransform(pts1,pts2)

    dst = cv2.warpAffine(img, M, (cols, rows))

    cv2.imshow("dst", dst)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

if __name__ == '__main__':
    main()

出力画像の対応点にアフィン変換された結果になりました。

f:id:sanvarie:20210116093405p:plain

射影変換

任意の四角形を任意の四角形に変形することができる変換が射影変換です。はねぴょんの画像では少しわかりにくいので、射影変換に関しては以下の画像を使用します。

f:id:sanvarie:20210116094427p:plain

射影変換を使うことで雑誌を真上から見たように補正することができます。射影変換を行うには入力画像と出力画像の対応点の座標が少なくとも4組必要です。今回は以下のように対応点を設定してみました。

f:id:sanvarie:20210116095309p:plain

# -*- coding: utf-8 -*-
import cv2
import numpy as np

def main():
    img = cv2.imread(r"D:\blog\data\opencv\book.jpg")
    rows, cols, _ = img.shape

    pts1 = np.float32([[765,521],[3135,551],[167, 2023],[4025, 1849]])
    pts2 = np.float32([[0,0],[cols,0],[0,rows],[cols, rows]])

    M = cv2.getPerspectiveTransform(pts1,pts2)

    dst = cv2.warpPerspective(img, M, (cols, rows))
    cv2.imwrite(r"D:\blog\data\opencv\book2.jpg", dst) # 保存

if __name__ == '__main__':
    main()

文字はぼやけてしまっていますが、雑誌を真上から見たように補正することができました。

f:id:sanvarie:20210116103858p:plain

さいごに

実は意外と GIS をつかっていると画像を幾何学変換する機会はあるかと思います。私の場合は地図上に水道の図面を幾何学変換させていました。その図面を下図に水道管などを作図するイメージですかね。まだ基本的なことしか紹介していませんが、OpenCV を使えば簡単に画像の幾何学変換を行うことができます。興味のある方はぜひ試してみてください。

Python の認証プロキシ設定について

さて、本日は Python の認証プロキシ設定について書いてみようと思います。社内環境で WebAPI などを叩こうと思ったときに認証プロキシにひっかかることがあるかと思います。今回は認証プロキシの突破方法について書いてみようと思います。

実行環境

Windows 10 64bit
Python 3.6.6

認証プロキシ設定方法

  1. 環境変数を設定する
  2. プログラム内で環境変数を定義

環境変数を設定する

環境変数を設定することで認証プロキシを通過することができます。

set HTTP_PROXY=http://user:pass@server:port

例えば、以下のような設定の場合、環境変数に設定する値は「http://gis:gis2021@hoge:8080」のようになります。

  • ユーザ:gis
  • パスワード:gis2021
  • サーバー名:hoge
  • ポート:8080

プログラム内で環境変数を定義

環境変数を設定したくない場合、以下のようにプログラム内で環境変数を定義すれば認証プロキシを通過することができます。

import os
os.environ["http_proxy"] = "http://user:pass@server:port"

さいごに

認証プロキシにお困りの方はぜひ参考にしてみてください。短いですが、本日は以上です。

Python と OpenCV で画像処理をしてみよう!

さて、本日は OpenCV を使って画像処理をしてみようと思います。私は最近画像を扱うことが多く、OpenCV を使ってもう少し画像関連の知識の幅を広げようと思っています。

今回のテーマは GIS とは直接は関係ないのですが、GIS の世界では航空写真など多くの画像を扱います。GIS のお仕事に興味がある方は画像の扱いにも慣れておくことをおすすめします。

画像関連のエントリーをいくつか書いていますので、興味のある方はぜひ読んでみてください。

www.gis-py.com

www.gis-py.com

www.gis-py.com

OpenCV とは?

オープンソースのコンピュータビジョン向けのライブラリですね。要するに画像や動画用のライブラリですね。OpenCV は元々 C/C++ で使えるライブラリでしたが、Python でも動作するようになりました。ありがたい話です。

今回試すこと

まずは以下の基本的な機能を試してみたいと思います。

  1. 画像を読み込む
  2. 興味領域(ROI)の抽出
  3. 画像のブレンド
  4. 画像の重ね合わせ
  5. 画像から顔と目を検出

実行環境

Windows 10 64bit
Python 3.6.6

使用する画像

大田区公式PRキャラクターの「はねぴょん」の画像を使用します。

f:id:sanvarie:20210111155935p:plain

サンプルコード

画像を読み込む
# -*- coding: utf-8 -*-
import cv2

def main():
    img = cv2.imread(r"D:\blog\data\opencv\hane.jpg", 1) # カラーで読込
    cv2.imshow("test", img)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

if __name__ == '__main__':
    main()

このようにはねぴょんの画像が読み込まれたことが確認できました。

f:id:sanvarie:20210111160141p:plain

興味領域(ROI)の抽出

ROIとは、Region of Interestの略で、処理対象とする領域のことですね。以下を ROI としてみます。「ペイント」アプリを使えば ROI とする座標がわかります。

f:id:sanvarie:20210111162318p:plain

# -*- coding: utf-8 -*-
import cv2

def main():
    img = cv2.imread(r"D:\blog\data\opencv\hane.jpg", 1) # カラーで読込
    roi = img[68:128, 35:143] # ROI領域
    cv2.imwrite(r"D:\blog\data\opencv\hane_roi.jpg", roi) # 保存

if __name__ == '__main__':
    main()

このように ROI を抽出することができました。 f:id:sanvarie:20210111163959p:plain

画像のブレンド

以下の画像をはねぴょんの画像にブレンディングします(はねぴょんの画像の重みを0.7、上の画像の重みを0.3に設定)。

f:id:sanvarie:20210111171544p:plain

# -*- coding: utf-8 -*-
import cv2

def main():
    img1 = cv2.imread(r"D:\blog\data\opencv\hane.jpg")
    img2 = cv2.imread(r"D:\blog\data\opencv\opencv_logo.png")

    # 画像の大きさをあわせる
    img1 =cv2.resize(img1,(200,200))
    img2 =cv2.resize(img2,(200,200))

    img1 = cv2.cvtColor(img1, cv2.COLOR_BGR2RGB)
    img2 = cv2.cvtColor(img2, cv2.COLOR_BGR2RGB)

    blended = cv2.addWeighted(src1=img1,alpha=0.7,src2=img2,beta=0.3,gamma=0)

    cv2.imshow("blended",blended)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

if __name__ == '__main__':
    main()

このように画像がブレンドされたことがわかります。

f:id:sanvarie:20210111171929p:plain

画像の重ね合わせ

今度はブレンドではなく、画像を重ね合わせてみようと思います。

# -*- coding: utf-8 -*-
import cv2

def main():
    img1 = cv2.imread(r"D:\blog\data\opencv\hane.jpg")
    img2 = cv2.imread(r"D:\blog\data\opencv\opencv_logo.png")

    # 片方の画像を小さくする
    img2 =cv2.resize(img2,(100,100))

    #img1 = cv2.cvtColor(img1, cv2.COLOR_BGR2RGB)
    #img2 = cv2.cvtColor(img2, cv2.COLOR_BGR2RGB)

    large_img = img1
    small_img = img2

    large_img[0:small_img.shape[0], 0:small_img.shape[1]] = small_img

    cv2.imshow("large_img", large_img)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

if __name__ == '__main__':
    main()

このように画像が重なったことがわかります。 f:id:sanvarie:20210111172544p:plain

画像から顔と目を検出

OpenCV では顔,目,笑顔などを検出することができます(検出させるものを学習させることも可能です)。今回は顔と目を検出してみようと思いますが、はねぴょんでは検出ができないので以下の画像を使用します。

f:id:sanvarie:20210111173628j:plain

# -*- coding: utf-8 -*-
import os
import cv2

# OpenCVのディレクトリ
opencv_path = r"C:\Python\Lib\site-packages\cv2\data"

def main():

    # 顔検出用データ
    face_cascade = cv2.CascadeClassifier(os.path.join(opencv_path,
                                         "haarcascade_frontalface_default.xml"))

    # 目検出用データ
    eye_cascade = cv2.CascadeClassifier(os.path.join(opencv_path,
                                        "haarcascade_eye.xml"))

    # 画像読込
    img = cv2.imread(r"D:\blog\data\opencv\CR7.jpg")
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

    # 顔検出
    faces = face_cascade.detectMultiScale(gray,
                                          scaleFactor=1.3,
                                          minNeighbors = 5,
                                          minSize=(30, 30))
    for (x,y,w,h) in faces:
        # 顔検出した部分に枠を描画
        img = cv2.rectangle(img,(x,y),(x+w,y+h),(255,0,0),2)
        roi_gray = gray[y:y+h, x:x+w]
        roi_color = img[y:y+h, x:x+w]

        # 目検出
        eyes = eye_cascade.detectMultiScale(roi_gray)
        for (ex,ey,ew,eh) in eyes:
            # 目検出した部分に枠を描画
            cv2.rectangle(roi_color,(ex,ey),(ex+ew,ey+eh),(0,255,0),2)

    cv2.imshow("img", img)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

if __name__ == '__main__':
    main()

このように顔と目を検出することができることがわかります。 f:id:sanvarie:20210111181036p:plain

さいごに

まだまだ OpenCV の基本機能を触れただけですが、簡単なコードで画像から顔と目を検出できたり非常に便利なライブラリということがわかりました。OpenCV はすごく多機能なライブラリのため、次回以降でもOpenCV の紹介をしていこうと思います。本日は以上です。

Python で緯度経度を地域メッシュコードに変換する方法

明けましておめでとうございます。

昨年は一昨年に続き変化の多い年でした。新しい分野に挑戦することになり四苦八苦していますが、頑張ろうと思います。

さて、本日は Python で緯度経度を地域メッシュコードに変換してみようと思います。以下のエントリーで地域メッシュコードを緯度経度に変換する方法を紹介していますが、今回についてその逆を紹介します。

www.gis-py.com

地域メッシュとは

前回エントリーに記載しましたので、こちらを参照してください。

www.gis-py.com

やり方

地域メッシュの考え方が以下のサイトで公開されています。これを参考にしてコードを作ってみようと思います。

https://www.stat.go.jp/data/mesh/pdf/gaiyo1.pdf

以下にサイトの画面ショットをのせておきます。

f:id:sanvarie:20210107211517p:plain

f:id:sanvarie:20210107211554p:plain

使用する緯度経度

こちらで取得したJR蒲田駅の緯度経度を使用したいと思います。なぜ蒲田かは秘密です。

緯度:35.562479
経度:139.716073

f:id:sanvarie:20210107211845p:plain

蒲田駅の地域メッシュコードをESRIジャパンさんが公開しているデータで確認してみます。

1次メッシュ:5339
2次メッシュ:533925
3次メッシュ:53392577

ということがわかりました。

f:id:sanvarie:20210107212512p:plain

実行環境

Windows 10
Python 3.6.6

サンプルコード

1~3次メッシュまで対応させたサンプルコードを作成しました。

# -*- coding: utf-8 -*-
def latlon2mesh(lat, lon):
    #1次メッシュ上2けた
    quotient_lat, remainder_lat = divmod(lat * 60, 40)
    first2digits = str(quotient_lat)[0:2]

    #1次メッシュ下2けた
    last2digits = str(lon - 100)[0:2]
    remainder_lon = lon - int(last2digits) - 100

    #1次メッシュ
    first_mesh = first2digits + last2digits

    #2次メッシュ上1けた
    first1digits, remainder_lat = divmod(remainder_lat, 5)

    #2次メッシュ下1けた
    last1digits, remainder_lon = divmod(remainder_lon * 60, 7.5)

    #2次メッシュ
    second_mesh = first_mesh + str(first1digits)[0:1] + str(last1digits)[0:1]

    #3次メッシュ上1けた
    first1digits, remainder_lat = divmod(remainder_lat * 60, 30)

    #3次メッシュ下1けた
    last1digits, remainder_lon = divmod(remainder_lon * 60, 45)

    #3次メッシュ
    third_mesh = second_mesh + str(first1digits)[0:1] + str(last1digits)[0:1]

    print("1次メッシュ:" + first_mesh)
    print("2次メッシュ:" + second_mesh)
    print("3次メッシュ:" + third_mesh)

if __name__ == '__main__':
    latlon2mesh(35.7007777, 139.71475)

結果を確認すると想定通りの地域メッシュコードを取得することができました。

1次メッシュ:5339
2次メッシュ:533925
3次メッシュ:53392577

f:id:sanvarie:20210107213306p:plain

さいごに

「緯度経度を地域メッシュコードに変換」と聞くとなんとなく難解そうな気がしますが、コードにするととても短くシンプルなプロセスということがわかるかと思います。地域メッシュコードは GIS の世界ではよく使用するものなので、今まで聞いたことなかったという方はぜひ色々調べてみてください。本年もよろしくお願いします。

Python で交通事故統計データを可視化してみよう!

さて、本日は交通事故統計データの可視化をしてみようと思います。本エントリーですが、以下二つのエントリーで行った作業を前提として進めていきます。興味がありましたらぜひ以下エントリーも読んでみてください。

www.gis-py.com

www.gis-py.com

交通事故統計データとは

警察庁が公開しているオープンデータです。その名の通り交通事故の統計データですね。

www.npa.go.jp

交通事故統計データ概要

2019年の交通事故統計データが公開されています。こちらに記載されていますが、以下三つのデータが CSV 形式で存在します。

  • 本票

    • 交通事故1件につき1レコード
    • 記録事項は、交通事故の内容に関する事項及び交通事故に関与した者(当事者A、当事者B)に関する事項
  • 補充票

    • 本票以外の交通事故に関与した者1人について1レコード
    • 記録事項は、本票に記録されない交通事故に関与した者のうち、死亡若しくは負傷した者又は車両等の運転者で死傷がなく死亡事故に関与した者に関する事項
  • 高速票

    • 高速自動車国道及び道路交通法第110条第1項により国家公安委員会が指定する自動車専用道路における交通事故1件について1レコード
    • 記録事項は、交通事故の発生地点、道路構造等に関する事項

要するに本票がメインのレコードで補充票、高速票がそれに紐づくレコードといった感じでしょうか。今回は「本票」「高速票」のデータを使用して可視化を行いたいと思います。補充票に関しては可視化はしませんが、一応データの格納だけはしようと思います。

手順

  1. 本票、補充票、高速票を SQLite に格納する
  2. 本票のデータを加工
  3. 加工したデータからシェープファイルに変換するデータを取得
  4. PyShp を使ってポイントのシェープファイルを作成
  5. 作成したシェープファイルを ArcGIS Online で可視化する

PyShp と ArcGIS Online に関しては以下のエントリーで紹介しているので、興味のある方がいましたらぜひ読んでみてください。

www.gis-py.com

www.gis-py.com

実行環境

Windows10 64bit
Python3.6.6
ArcGIS Online

手順1(本票、補充票、高速票を SQLite に格納) サンプルコード

まずは データを SQLite に格納します。前回エントリーと同じ要領ですね。

# -*- coding: utf-8 -*-
import sqlite3
from os.path import join, split
from glob import glob
import pandas as pd

# 本票、補充票、高速票格納先
csv_path = r"D:\blog\data\csv"

# SQLite 格納先
db_path = r"D:\blog\data\sqlite\test.db"

def insert_csv():
    # CSV ファイルのディレクトリ取得
    csv_files = glob(join(csv_path, "*.csv"))

    for csv in csv_files:

        table_name = split(csv)[1] # テーブル名作成

        table_name = table_name[0:-4] # テーブル名から拡張子を削除

        df = pd.read_csv(csv, dtype=object, encoding="SHIFT-JIS") # CSV 読込

        with sqlite3.connect(db_path) as conn:
            df.to_sql(table_name, con=conn) # SQLiteにCSVをインポート

if __name__ == '__main__':
    insert_csv()

前回エントリーで作成したDBに「hojuhyo_2019」「honhyo_2019」「kosokuhyo_2019」というテーブルが追加されました。

f:id:sanvarie:20201230143715p:plain

f:id:sanvarie:20201230144550p:plain

手順2(手順「1」で格納した本票のデータを加工) サンプルコード

オープンデータにありがちですが、公開されているデータがあまりきれいな状態ではなく、少しデータの整形をしたいと思います。ここのパートはややこしく、説明がわかりづらいかもしれませんがご容赦ください。以下3つのステップを踏んでデータ加工を行います。

  1. 「honhyo_2019」テーブルの「死者数」、「負傷者数」カラムに入っている不要な0を削除
  2. 一意になるキーを作成
  3. 「honhyo_2019」テーブルの「路線コード」カラムを名称に変換
1.「honhyo_2019」テーブルの「死者数」、「負傷者数」カラムに入っている不要な0を削除

負傷者が一人だったら、001のような値が入っているのですが、このような形でデータが格納されているのは一体なぜなのでしょうか。ちょっとよくわからないですね。。。

 update honhyo_2019
    set 死者数 = case
                    when substr(死者数,1,2) = '00' then
                        substr(死者数,3,1)
                    when substr(死者数,1,1) = '0' and substr(死者数,2,1) <> '0' then
                        substr(死者数,2,2)
                end,
        負傷者数 = case
                      when substr(負傷者数,1,2) = '00' then
                          substr(負傷者数,3,1)
                      when substr(負傷者数,1,1) = '0' and substr(負傷者数,2,1) <> '0' then
                          substr(負傷者数,2,2)
                  end

加工前

f:id:sanvarie:20201230153528p:plain

加工後

f:id:sanvarie:20201230153647p:plain

2.一意になるキーを作成

本票番号というカラムがキーになると思ったのですが、色々調べてると違うことがわかりました。。。「honhyo_2019」「kosokuhyo_2019」テーブルともに「都道府県コード」、「警察署等コード」、「本票番号」カラムを結合させてキーを作りたいと思います。

事故番号」というカラムを作成します。

alter table honhyo_2019 add column 事故番号[text];
alter table kosokuhyo_2019 add column 事故番号[text];

以下で当該カラムに「都道府県コード」、「警察署等コード」、「本票番号」カラムの結合した文字列を格納します。

update honhyo_2019 
   set 事故番号 = 都道府県コード || 警察署等コード || 本票番号

update kosokuhyo_2019
   set 事故番号 = 都道府県コード || 警察署等コード || 本票番号

これをキーとして使用したいと思います。

3.「honhyo_2019」テーブルの「路線コード」カラムを名称に変換

最初は「本票_路線コード」テーブルと結合して路線名を取得しようと思ったのですが、ちょっと困ったことに当該テーブルのレコードが以下のようになっていて結合ができないので、カラム追加をしてそこに路線名を格納しようと思います。

f:id:sanvarie:20201230172219p:plain

まずは「路線名」カラム追加します。

alter table honhyo_2019 add column 路線名[text];

次に追加したカラムに路線名を格納します。

update honhyo_2019
   set 路線名 = case 
                  when substr(路線コード,1,4) between '0001' and '0999' then '一般国道(国道番号)'
                  when substr(路線コード,1,4) between '1000' and '1499' then '主要地方道-都道府県道'
                  when substr(路線コード,1,4) between '1500' and '1999' then '主要地方道-市道'
                  when substr(路線コード,1,4) between '2000' and '2999' then '一般都道府県道'
                  when substr(路線コード,1,4) between '3000' and '3999' then '一般市町村道'
                  when substr(路線コード,1,4) between '4000' and '4999' then '高速自動車国道'
                  when substr(路線コード,1,4) between '5000' and '5499' then '自動車専用道-指定'
                  when substr(路線コード,1,4) between '5500' and '5999' then '自動車専用道-その他'
                  when substr(路線コード,1,4) between '6000' and '6999' then '道路運送法上の道路'
                  when substr(路線コード,1,4) between '7000' and '7999' then '農(免)道'
                  when substr(路線コード,1,4) between '8000' and '8499' then '林道'
                  when substr(路線コード,1,4) between '8500' and '8999' then '港湾道'
                  when substr(路線コード,1,4) between '9000' and '9499' then '私道'
                  when substr(路線コード,1,4) = '9500' then 'その他'
                  when substr(路線コード,1,4) = '9900' then '一般の交通の用に供するその他の道路'
             end

これでデータの加工は終了です。

手順3(手順「2」で加工したデータをもとにシェープファイルに変換するデータを取得) サンプルコード

「honhyo_2019」「kosokuhyo_2019」テーブルから必要なカラムを選定し、コードを名称に変換します。

以下のSQLで取得するデータをシェープファイルに変換します(「手順4」参照)。カラム数がかなり多いので、必要最低限のカラムのみを取得しています。

select a.事故番号   as 事故番号,
       c.都道府県名 as 都道府県,
       e.区分       as 事故内容,
       a.死者数     as 死者数,
       a.負傷者数   as 負傷者数,
       a.路線名     as 路線名,
       a.発生日時  年 || '/' || a.発生日時  月 || '/' || a.発生日時  月 || ' ' || a.発生日時  時 || ':' ||a.発生日時  分 as 発生日時, 
       f.区分       as 天候,
       g.区分       as 高速道路事故発生場所,
       h.区分       as 道路構造,
       a.地点 緯度(北緯) as 緯度,
       a.地点 経度(東経) as 経度
  from honhyo_2019 a
    left outer join kosokuhyo_2019 b
      on a.事故番号 = b.事故番号
    left outer join 本票_補充票_高速票_都道府県コード c
      on a.都道府県コード = c.コード
    left outer join 本票_事故内容 e
      on a.事故内容 = e.コード
    left outer join 本票_天候 f
      on a.天候 = f.コード
    left outer join 高速票_道路区分 g
      on b.道路区分 = g.コード
    left outer join 高速票_道路構造 h
      on b.道路構造 = h.コード

こんな感じのデータを取得することができました。

f:id:sanvarie:20201230175502p:plain

手順4(PyShp を使ってポイントのシェープファイルを作成) サンプルコード

「手順3」で使用している SQL を使用してポイントのシェープファイルを作成を作成します。PyShp に関しては以下のエントリーを参考にしてみてください。

www.gis-py.com

# -*- coding: utf-8 -*-
import math
from decimal import Decimal, ROUND_HALF_UP
from os.path import join
import sqlite3
import shapefile
import pandas as pd

# シェープファイル作成先
shape_path = r"D:\blog\data\shape"

# SQLite 格納先
db_path = r"D:\blog\data\sqlite\test.db"

# epsg
epsg = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]'

def dms2deg(dms):
    """緯度経度を度分秒から度へ変換します"""

    h = dms[0]
    m = dms[1]
    s = dms[2]
    deg = Decimal(str(int(h) + (int(m) / 60) + (float(s) / 3600))).quantize(Decimal('0.0001'), rounding=ROUND_HALF_UP)
    return deg

def get_data():
    """SQLiteからデータを取得します"""

    conn = sqlite3.connect(db_path)
    sql = """select a.事故番号  as 事故番号,
                   c.都道府県名 as 都道府県,
                 e.区分       as 事故内容,
                 a.死者数     as 死者数,
                 a.負傷者数   as 負傷者数,
                 a.路線名     as 路線名,
                 a.発生日時  年 || '/' || a.発生日時  月 || '/' || a.発生日時  月 || ' ' || a.発生日時  時 || ':' ||a.発生日時  分 as 発生日時,
                 f.区分       as 天候,
                   g.区分       as 高速道路事故発生場所,
                   h.区分       as 道路構造,
                 a.地点 緯度(北緯) as 緯度,
                 a.地点 経度(東経) as 経度
              from honhyo_2019 a
                left outer join kosokuhyo_2019 b
                on a.事故番号 = b.事故番号
              left outer join 本票_補充票_高速票_都道府県コード c
                on a.都道府県コード = c.コード
                left outer join 本票_事故内容 e
                on a.事故内容 = e.コード
              left outer join 本票_天候 f
                on a.天候 = f.コード
              left outer join 高速票_道路区分 g
                on b.道路区分 = g.コード
              left outer join 高速票_道路構造 h
                on b.道路構造 = h.コード
            """

    return pd.read_sql_query(sql=sql, con=conn)

def create_shape():
    """シェープファイルを作成します"""

    # SQLite からデータを取得
    df = get_data()

    # データフレームのカラム取得
    column_list = df.columns.values

    with shapefile.Writer(join(shape_path, "traffic_data.shp")) as w:
        # シェープファイルのフィールドを作成
        for column in column_list:
            if column == "緯度":
               w.field(column, "F", 18, 8)
            elif column == "経度":
                w.field(column, "F", 18, 8)
            else:
                w.field(column, "C", 100)

        # シェープファイル作成処理
        for _, row in df.iterrows():

            lat = 0
            lon = 0
            attributes = []

            for column in column_list:

                if column == "緯度":
                    lat = row[column]
                    lat = dms2deg((lat[0:2], lat[2:4], lat[4:6] + "." + lat[6:]))
                    attributes.append(lat)

                elif column == "経度":
                    lon = row[column]
                    lon = dms2deg((lon[0:3], lon[3:5], lon[5:7] + "." + lon[7:]))
                    attributes.append(lon)

                else:
                    if row[column] == None:
                        attributes.append("")
                    else:
                        attributes.append(row[column])


            # ポイントを作成
            w.point(float(lon), float(lat))
            w.record(attributes[0],attributes[1],attributes[2],attributes[3],
                     attributes[4],attributes[5],attributes[6],attributes[7],
                     attributes[8],attributes[9],attributes[10],attributes[11])

    # プロジェクションファイル作成
    with open("%s.prj" % join(shape_path, "traffic_data"), "w") as prj:
        prj.write(epsg)

if __name__ == '__main__':
    create_shape()

このようにシェープファイルが作成されていることが確認できました。やはり件数が多いだけあって容量が大きいですね。

f:id:sanvarie:20201230202558p:plain

手順5(ArcGIS Online で可視化)

シェープファイル を ArcGIS Online にアップロードしてみます。

ArcGIS Online の使い方に関してはESRIジャパンさんが主催した「ArcGIS 開発者のための最新アプリ開発塾 2020」の「Web GIS 基礎 ~ArcGIS Online を使ってみよう!~」を参考にしてみてください。

community.esri.com

路線ごとに色分けした結果を表示してみました。正しくデータがプロットされていることがわかります。

f:id:sanvarie:20201231092538p:plain

f:id:sanvarie:20201231092620p:plain

天候ごとにも色分けしてみました。

f:id:sanvarie:20201231092449p:plain

事故ごとの負傷者数で色分けなんかもできます。

f:id:sanvarie:20201231123750p:plain

データをより細かく分析するにはどうしたらいいのか?

以下のようなことがしたい場合は ArcGIS Pro などのデスクトップアプリケーションを導入することをお勧めします。

  1. 事故多発地帯の分析
  2. 同じ地点に短期間限定で事故が多発している箇所の分析
  3. 天候ごとに事故が発生しやすい箇所の分析
  4. etc

さいごに

オープンデータにありがちですが、データが素直に取り込めばいいという状態になっておらず色々苦労しました。交通事故統計データだけでなく、日本のオープンデータすべてが利用者にとってもう少し優しい形になってほしいですね。

手順が多いエントリーになってしまいましたが、何とか交通事故統計データの可視化までもっていくことができました。ただ、GIS的にはここからが本番でデータ分析してこそ、そのデータに価値が生まれるかと思います。もし ArcGIS Pro でこのデータを使った分析方法を知りたいという方がいらっしゃいましたら本エントリーにコメントを残していただければと思います。

本エントリーが今年最後です。今年もありがとうございました。来年もよろしくお願いします。それではよいお年をお迎えください。

Python で CSV ファイルを SQLite にインポートする方法

さて、本日は CSV ファイルを SQLite にインポートしてみようと思います。前回のエントリーで「交通事故統計情報のオープンデータ」のコード表を CSV に変換してみました。本エントリーではその CSV を SQLite にインポートしてみようと思います。

PDF を CSV に変換する方法について興味のある方はぜひ前回のエントリーを読んでみてください。

www.gis-py.com

また、以下のエントリーで SQLite に関して書いていますので、こちらも興味がありましたらぜひ読んでみてください。

www.gis-py.com

対象データ

前回のエントリーで出力した CSV です。 これらをファイルごとにインポートしていこうと思います。

f:id:sanvarie:20201229171900p:plain

f:id:sanvarie:20201229172351p:plain

やり方

  1. Pandas で データフレームとして CSV を読み込む
  2. 「1」のデータフレームをSQLiteにインポート

実行環境

Windows10 64bit
Python3.6.6

サンプルコード

CSV ファイルを SQLite にインポートするサンプルです。すごく短いコードで済みましたね。

# -*- coding: utf-8 -*-
import sqlite3
from os.path import join, split
from glob import glob
import pandas as pd

csv_path = r"D:\blog\data\pdf2csv\csv"
db_path = r"D:\blog\data\sqlite\test.db"

def insert_csv():
    # CSV ファイルのディレクトリ取得
    csv_files = glob(join(csv_path, "*.csv"))

    for csv in csv_files:

        table_name = split(csv)[1] # テーブル名作成

        table_name = table_name[0:-4] # テーブル名から拡張子を削除

        df = pd.read_csv(csv, dtype=object) # CSV 読込

        with sqlite3.connect(db_path) as conn:
            df.to_sql(table_name, con=conn) # SQLiteにCSVをインポート

if __name__ == '__main__':
    insert_csv()

結果を確認するとこのように CSV がSQLiteにインポートされたことがわかります。

f:id:sanvarie:20201229185546p:plain

f:id:sanvarie:20201229183843p:plain

さいごに

CSV ファイルの SQLite へのインポートは思ったよりも簡単にできることがわかりました(Pandas は本当に便利なライブラリですね)。これでようやく「交通事故統計情報のオープンデータ」の可視化にとりかかれそうです。ちょっとまだすべてのデータを確認したわけではないので何とも言えないのですが、結構大変そうですね。ただ、今年の年末年始はコロナでたっぷり時間がありそうですし楽しみながら取り掛かってみようと思います。次回のエントリーもぜひ読んでみてください。