気象データを Python でスクレイピングする方法

さて、本日は気象庁の HP から気象データをスクレイピングしようと思います。

スクレイピングとは

ウェブスクレイピングとは、ウェブサイトから情報を抽出するコンピュータソフトウェア技術のこと。ウェブ・クローラーあるいはウェブ・スパイダーとも呼ばれる。 通常このようなソフトウェアプログラムは低レベルのHTTPを実装することで、もしくはウェブブラウザを埋め込むことによって、WWWのコンテンツを取得する。（出典：Wikipedia）

関連エントリー

過去のエントリーで他サイトのスクレイピングについて書いていますので興味のある方はぜひ参考にしてみてください。

www.gis-py.com

www.gis-py.com

事の発端

以下のサイトから東京の一時間ごとの気象データをダウンロードしようとしたのですが、一回のダウンロードの容量制限があり数年間分のデータをダウンロードするのは少し面倒だなと思い今回のエントリーに至りました。

www.data.jma.go.jp

スクレイピングするサイト

上記のダウンロードサイトではなく過去の気象データ検索のサイトを利用してスクレイピングを行おうと思います。

www.data.jma.go.jp

使用するライブラリ

Beautiful Soup というライブラリを使用します。Beautiful Soupは、HTMLおよびXMLドキュメントを解析するためのPythonパッケージです。 HTMLからデータを抽出するために使用できる解析済みページの解析ツリーを作成します。これはWebスクレイピングに役立ちます。（出典：Wikipedia）

ダウンロードするデータ

今回は以下データをダウンロードします。場所、期間を変えたいという方は掲載するサンプルコードを少し変えれば問題ないかと思います。ただ、データの単位を変えたい場合は少し工夫が必要になりますが今回はそれについては割愛します。

• 場所・・・東京
• 期間・・・2016年~2020年
• データの単位・・・時間ごと

環境

Windows10 64bit
Python3.8.5

サンプルコード

気象データをスクレイピングするサンプルコードです。

# -*- coding: utf-8 -*-
import os
import datetime
import csv
import urllib.request
from bs4 import BeautifulSoup

def str2float(weather_data):
    try:
        return float(weather_data)
    except:
        return 0

def scraping(url, date):

    # 気象データのページを取得
    html = urllib.request.urlopen(url).read()
    soup = BeautifulSoup(html)
    trs = soup.find("table", { "class" : "data2_s" })

    data_list = []
    data_list_per_hour = []

    # table の中身を取得
    for tr in trs.findAll('tr')[2:]:
        tds = tr.findAll('td')

        if tds[1].string == None:
            break;

        data_list.append(date)
        data_list.append(tds[0].string)
        data_list.append(str2float(tds[1].string))
        data_list.append(str2float(tds[2].string))
        data_list.append(str2float(tds[3].string))
        data_list.append(str2float(tds[4].string))
        data_list.append(str2float(tds[5].string))
        data_list.append(str2float(tds[6].string))
        data_list.append(str2float(tds[7].string))
        data_list.append(str2float(tds[8].string))
        data_list.append(str2float(tds[9].string))
        data_list.append(str2float(tds[10].string))
        data_list.append(str2float(tds[11].string))
        data_list.append(str2float(tds[12].string))
        data_list.append(str2float(tds[13].string))

        data_list_per_hour.append(data_list)

        data_list = []

    return data_list_per_hour

def create_csv():
    # CSV 出力先ディレクトリ
    output_dir = r"D:\blog\data\weather"

    # 出力ファイル名
    output_file = "weather.csv"

    # データ取得開始・終了日
    start_date = datetime.date(2016, 1, 1)
    end_date   = datetime.date(2020, 12, 31)

    # CSV の列
    fields = ["年月日", "時間", "気圧（現地）", "気圧（海面）",
              "降水量", "気温", "露点湿度", "蒸気圧", "湿度",
              "風速", "風向", "日照時間", "全天日射量", "降雪", "積雪"] # 天気、雲量、視程は今回は対象外とする

    with open(os.path.join(output_dir, output_file), 'w') as f:
        writer = csv.writer(f, lineterminator='\n')
        writer.writerow(fields)

        date = start_date
        while date != end_date + datetime.timedelta(1):

            # 対象url（今回は東京）
            url = "http://www.data.jma.go.jp/obd/stats/etrn/view/hourly_s1.php?" \
                  "prec_no=44&block_no=47662&year=%d&month=%d&day=%d&view="%(date.year, date.month, date.day)

            data_per_day = scraping(url, date)

            for dpd in data_per_day:
                writer.writerow(dpd)

            date += datetime.timedelta(1)

if __name__ == '__main__':
    create_csv()

やはり少々データ量が多いだけあって処理完了まで少しだけ時間がかかりました。ただ、結果はばっちりですね！

さいごに

気象庁のデータのダウンロードページは便利なのですが、様々な場所のデータを数年間分欲しいといった場合に手動でダウンロードを行うというのはあまり現実ではないかと思います。その際に本エントリーで紹介したようなスクレイピングが役に立つかと思います。また、世の中にある様々な Web サイトのデータ（例：食べログ）がほしいといった場合もスクレイピングをすることで今回のようにデータを取得することができます。興味のある方はぜひチャレンジしてみてください！本日は以上です。

機械学習を使って不動産取引価格の予測値を可視化してみよう！

さて、本日は機械学習ですね。今までのエントリーと異なりデータを加工・可視化するのではなく、収集したデータを使って将来の予測値を可視化してみようと思います。

機械学習とは

機械学習（きかいがくしゅう、英: Machine Learning）とは、経験からの学習により自動で改善するコンピューターアルゴリズムもしくはその研究領域で[1][2]、人工知能の一種であるとみなされている。「訓練データ」もしくは「学習データ」と呼ばれるデータを使って学習し、学習結果を使って何らかのタスクをこなす。例えば過去のスパムメールを訓練データとして用いて学習し、スパムフィルタリングというタスクをこなす、といった事が可能となる。（出典:Wikipedia）

今回使用する機械学習ライブラリ

今回は scikit-learn を使おうと思います。かなり広く使われているライブラリのようで、おそらく私のような初学者はまずこれを使うのがいいのではないかと思います。

scikit-learn.org

今回使用するデータ

不動産取引価格ダウンロード 国土交通省のデータを使用します。

2010年～2020年の横須賀市のデータをダウンロードしてみます。

データはこんなイメージですね。

各フィールド

データに関しては以下を参照してください。全てのフィールドは使用せずに必要なもののみをピックアップして機械学習を行おうと思います。

フィールド名	説明
No	連番
種類	「宅地(土地)」「中古マンション等」など。今回は「宅地(土地と建物)」を利用する。
地域	「住宅地」「商業値」など。 今回は「住宅地」を利用する。
市区町村コード	市区町村コード
都道府県名	都道府県名
市区町村名	市区町村名
地区名	地区名
最寄駅：名称	物件の最寄駅名
最寄駅：距離（分）	物件から最寄駅までの距離（分）
取引価格（総額）	物件の取引価格。これが今回推論したい値。
坪単価	坪単価
間取り	間取り
面積（㎡）	面積（㎡）
取引価格（㎡単価）	取引価格（㎡単価）
土地の形状	「長方形」「不整形」など
間口	間口
延床面積（㎡）	延床面積（㎡）
建築年	建築年。元号表記。
建物の構造	「ＳＲＣ」「木造」など
用途	「住宅」「共同住宅」など
今後の利用目的	「住宅」「事務所」など
前面道路：方位	前面道路の方位
前面道路：種類	前面道路の種類
前面道路：幅員（ｍ）	前面道路の幅員（ｍ）
都市計画	「商業地域」「第１種中高層住居専用地域」など
建ぺい率（％）	建ぺい率（％）
容積率（％）	容積率（％）
取引時点	「2016年第３四半期」など
改装	改装済みか否か
取引の事情等	特筆事項。「瑕疵有りの可能性」など

環境

Windows10 64bit
Python3.8.5
Tableau Desktop Public Edition 2021.1.0

手順

1. 横須賀市各駅の位置情報取得
2. データ加工
   1. 「種類」フィールドの値が「宅地(土地と建物)」のデータのみ使用
   2. 「地域」フィールドの値が「住宅地」のデータのみ使用
   3. 「用途」フィールドの値が「住宅」のデータのみ使用
   4. 今回推論に使用するフィールドにNULLがあるレコードの除去
   5. 「築年数」（元号表記）から築年数を算出
   6. 「築年数」（元号表記）から不正な値（「戦前」）を昭和15年に置換する
   7. 「最寄駅：距離（分）」の不正な値（「1H30?2H」など）を分に変換
   8. 「最寄駅：名称」の「逗子・葉山」を「逗子」に置換
   9. 外れ値の除去（1億5000万円より高い物件は今回は除去）
   10. 外れ値の除去（延床面積（㎡）が250㎡より大きいレコードは削除）
3. モデルの作成
4. モデルの評価

1.横須賀市各駅の位置情報取得

こちらのサイトのデータを利用し CSV を作成しました。（鎌倉や逗子が入っているというツッコミはなしでお願いします・・・）

2.データ加工

以下のコードで「手順」に記載した加工をすべて行っています。

import os
import pandas as pd

# 使用するファイル
input_dir = r"D:\blog\data\property"
property_data = "14201_20101_20204.csv"

# 出力するファイル
output_data = "output_data.csv"

# CSV読込
df = pd.read_csv(os.path.join(input_dir, property_data), encoding="cp932")

# 「種類」フィールドの値が「宅地(土地と建物)」のデータのみ使用
df = df[df['種類'] == "宅地(土地と建物)"]

# 「地域」フィールドの値が「住宅地」データのみ使用
df = df[df['地域'] == "住宅地"]

# 「用途」フィールドの値が「住宅」を含むデータのみ使用
df = df[df['用途'].str.contains("住宅") == True]

# 欠損値除去
df = df.dropna(subset=["最寄駅：距離（分）","延床面積（㎡）","建築年","建物の構造",'地域','用途'])

# 連番付与
serial_num = pd.RangeIndex(start=1, stop=len(df.index) + 1, step=1)
df["ID"] = serial_num

# 不要なフィールドを削除
df = df[["ID", "最寄駅：名称","最寄駅：距離（分）","延床面積（㎡）","建築年","建物の構造","取引価格（総額）"]]

# 最寄駅：距離（分）の不正な値を置換
df = df.replace({'最寄駅：距離（分）': {"1H?1H30": 60}})
df = df.replace({'最寄駅：距離（分）': {"1H30?2H": 90}})
df = df.replace({'最寄駅：距離（分）': {"2H?": 120}})
df = df.replace({'最寄駅：距離（分）': {"30分?60分": 30}})

# 建築年の不正な値を置換
df = df.replace({'建築年': {"戦前": "昭和15年"}})

# 建築年の元号表記を築年数に変換
df['年号'] = df['建築年'].str[:2]
df['和暦_年'] = df['建築年'].str[2:].str.replace('年','').astype(int)
df.loc[df['年号']=='昭和','築年数'] = 2021 - (df['和暦_年'] + 1925)
df.loc[df['年号']=='平成','築年数'] = 2021 - (df['和暦_年'] + 1988)
df.loc[df['年号']=='令和','築年数'] = 2021 - (df['和暦_年'] + 2019)

#「最寄駅：名称」の「逗子・葉山」を「逗子」に置換
df = df.replace({'最寄駅：名称': {"逗子・葉山": "逗子"}})

# 外れ値の除去（取引価格（総額）が1億50000万円より大きいレコードは削除）
df = df[df['取引価格（総額）'] <= 150000000]

# 外れ値の除去（延床面積（㎡）が250㎡より大きいレコードは削除）
df = df[df['延床面積（㎡）'].astype(int) <= 250]

# データの加工結果を出力
df.to_csv(os.path.join(input_dir, output_data))

出力した CSV を Tableau で確認してみます。各駅の取引価格（総額）を箱ひげ図で表示してみました。駅の順番はあやしいですが、「三崎口」や「三浦海岸」など三浦半島の先端に行くにしたがって金額が低くなることがわかります。また、サンプル数が少ないのですが、「逗子」や「鎌倉」は金額が高いですね。高級住宅街などがあるせいでしょうか。

価格と広さの関係性を確認してみようと思います。

R2乗値が約0.19ということは価格と広さにはあまり相関関係がないと言えそうです。

次に価格と築年数の関係性を確認してみようと思います。

R2乗値を見る限り、築年数が大きければ価格が下がる傾向にあると言えそうです。

また、「1.横須賀市各駅の位置情報取得」で作成した CSV を Tableau で読み込んでみました。

「取引価格（総額）」が高いほどポイントが赤色になるように設定しました。「逗子」や「鎌倉」エリアが高いことがわかります。

3.モデルの作成

まずは学習に必要なフィールドのみを抽出します。また、「1.横須賀市各駅の位置情報取得」で作成した CSV の緯度経度情報を output_data.csv にマージします。

import pandas as pd
import numpy as np

# 使用するファイル
input_dir = r"D:\blog\data\property"
output_data = "output_data.csv"
station_master = "yokosuka_station_master.csv"

# 前の手順で出力した CSV の読込
df_output = pd.read_csv(os.path.join(input_dir, output_data), encoding="utf-8")

# 駅マスタ読込
df_station = pd.read_csv(os.path.join(input_dir, station_master), encoding="utf-8")

# マージ
df = pd.merge(df_output, df_station, left_on='最寄駅：名称', right_on='station')
df = df.drop(["id", "station"], axis=1)

df.reset_index(drop=True, inplace=True)

# 出力変数
t = df["取引価格（総額）"]

# 学習に必要なカラムを抽出
x = df.iloc[:,[1,2,3,5,9]]

x.head()

ちなみに カラム抽出前のデータフレームは以下のように駅マスタの緯度経度が追加されています。

df.head()

object 型のカラムをダミー変数化します。

x = pd.get_dummies(x)
x.head()

データを訓練データと検証データに分割してモデルの学習を行います。

from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression

# 訓練データと検証データに分ける
x_train, x_val, t_train, t_val = train_test_split(x, t, test_size = 0.3, random_state = 0)

# 線形回帰モデル
model = LinearRegression()

model.fit(x_train, t_train)

モデルの評価を行います。訓練データで決定係数を確認しましたが、あまりいい値ではないですね。。

model.score(x_train, t_train)

検証データでも決定係数を確認してみます。

model.score(x_train, t_train)

予測値の推論

とりあえず今回は精度向上の試みはせずに先に進めようと思います。生成されたモデルを使って全てのデータの入力変数を使用して価格を推論します。

# 価格推論
pred = model.predict(x)

# 列名付与
y = pd.DataFrame(pred, columns=["predict"])

# 読み込んだ CSV の右端に推論結果を追加
results = pd.concat([df, y], axis=1)

# CSV 出力
results.to_csv(r"D:\blog\data\property\predict.csv", index=False)

4.モデルの評価

推論した価格と実際の価格を比較して誤差がどの程度かを Tableau で確認してみようと思います。

推論価格をX軸に、実際の価格をY軸にして散布図を作ってみました。本来は斜め45度のに分布するのが良いのですが、結果を見るとうーん、といった感じですね。

このように推論値がマイナスになっているものもあり、マイナスはいかんだろう感じです（笑）

他のを見ても推論値と実際の値とのずれが大きいことがわかります。

中にはかなり惜しいものもありますね。

「誤差」フィールドを作成して各駅ごとの推論値と実際値の誤差を表示してみました。「三浦海岸」「久里浜」「鎌倉」あたりが誤差が少ないことがわかります。

外れ値を確認してみると誤差が9千万以上ありますね。これはあまりいい結果とは言えないですね・・・

地図上で確認してみます。グレーになっているポイント（「三浦海岸」「久里浜」「鎌倉」など）が誤差が少ない駅ですね。

さいごに

はじめて本ブログで機械学習について触れてみましたがいかがでしたでしょうか。結果はいまいちでしたが、もう少しデータを詳細に分析し適切な学習モデルを作成することでより正確な推論値を得ることができるかと思います。これからの時代はただデータを GIS 上で可視化するだけでは差別化を行うことは難しく、機械学習などを用いて将来の予測値などを可視化する必要があるかと思っています。自分自身学ばなければならないことが多いのですが、時間をかけて学習する価値はあるかと思っています。これからも機械学習について継続的に本ブログで紹介しようかと思いますのでよろしくお願いします。本日は以上です。

Windows の Anaconda 環境に PyStan をインストールする方法

さて、本日は Windows の Anaconda 環境に PyStan をインストールする方法について書いてみようと思います。というのも先日この作業にはまったので備忘録として本エントリーを残そうと思います。本エントリーが同じ問題に直面した方の手助けになれば幸いです。

PyStan とは

Stan（スタン）はC++で書かれた統計的推論のための確率的プログラミング言語。 Stan言語では、対数確率密度関数を計算する命令型プログラムを使用して、（ベイジアン） 統計モデルを実装できる（出典: Wikipedia）。PyStan は Python で Stan を使用するためのライブラリです。

環境

Windows 10 64bit
Anaconda3-2020.11

Anaconda をインストール

以下のサイトで Windows 版の Anaconda をダウンロードしてインストールを行います。また、環境変数の設定も一緒に行います。

www.anaconda.com

私は以下に Anaconda 環境を入れました。
D:\anaconda3

環境変数 Path に以下を追加します。
D:\anaconda3
D:\anaconda3\Scripts

Anaconda 環境に PyStan をインストール

以下コマンドで Anaconda の仮想環境を構築することができます。今回は「1.仮想環境構築（ライブラリ無し）」で試してみます。

1.仮想環境構築（ライブラリ無し）

conda create --name XXX python=y.y

2.仮想環境構築（全ライブラリ含む）

conda create --name XXX python=y.y anaconda

仮想環境構築

Anaconda Prompt を開いて以下を実行します。Python のバージョンは3.7以上をインストールする必要があります。

conda create --name my-stan python=3.8.5

以下のコマンドで構築した仮想環境をアクティブにします。

conda activate my-stan

仮想環境をアクティブにしたので、必要なライブラリをインストールします。

conda install Cython Numpy Pandas matplotlib jupyter

pystanも入れてみます！

conda install pystan

C++のコンパイラーも必要になるので、以下コマンドでインストールします。

conda install libpython m2w64-toolchain -c msys2

インストール後、PYTHONPATH\envs\Lib\distutils に distutils.cfg というファイルが作成されます。もしなかったら手動で作成してみてください。ファイルの中身は以下のようになっています。

[build]           
compiler=mingw32 

Jupyter Notebook から作成した仮想環境に切り替えられるように以下コマンドを実行します。

ipython kernel install --user --name=my-stan

念のため再起動してこれで準備は完了です！

PyStan のインポートができました！使い方などはまたあらためて紹介できればと思います。本日は以上です。

Tableau の空間関数を使ってみよう！

さて、本日も前回のエントリーに続き Tableau について書いてみようと思います。以下のエントリーで Tableau を使って地図データを可視化してみましたが、今回はもう一歩踏み込んで空間関数を使ってみようと思います。

www.gis-py.com

使用するデータ

ESRIジャパンさんが主催した「ArcGIS 開発者のための最新アプリ開発塾 2020」で使用された 各店舗売上.csvを使用します。

community.esri.com

環境

Windows10 64bit
Tableau Desktop Public Edition 2021.1.0 64bit

使用する空間関数

• MakePoint・・・緯度経度からポイントを作成
• MakeLine・・・二つのポイントからラインを作成
• Distance・・・2 つのポイント間の距離を計測
• バッファー・・・ポイントのバッファーを作成

バッファーを作れるなんて素晴らしいですね！そして Tableau 恐るべしです。

MakePoint

緯度経度からポイントを作成します。

1.分析＞計算フィールドの作成

2.MAKEPOINT([緯度],[経度]) と入力＞OK

左側のデータ欄に「ポイント」が作成されたことがわかります。

3.「ポイント」をダブルクリック ポイントが可視化されました。

MakeLine

二つのポイントからラインを作成します。ここでは「蒲田駅」から「横須賀中央駅」の距離を計測してみようと思います。それぞれの緯度経度は以下です。

• 蒲田駅
  • 緯度 35.562479
  • 経度 139.716073
• 横須賀中央駅
  • 緯度 35.278699
  • 経度 139.670040

1.分析＞計算フィールドの作成

2.MAKELINE(MAKEPOINT(35.562479,139.716073),MAKEPOINT(35.278699,139.670040))＞OK

左側のデータ欄に「ライン」が作成されたことがわかります。

3.「ライン」をダブルクリック
蒲田駅～横須賀中央駅間でラインが可視化されました。

拡大してみましたが、駅が確認できません。

しかし、Tableau は以下のような手順で背景図を選ぶことができます。素晴らしいですね。

このように作成したラインが本当に蒲田駅～横須賀中央駅かどうか確認することができました。

Distance

2点間の距離を測ります。

1.分析＞計算フィールドの作成

2.DISTANCE(MAKEPOINT(35.562479,139.716073),MAKEPOINT(35.278699,139.670040),"km")と入力＞OK

左側のデータ欄に「距離」が作成されたことがわかります。

3.「距離」を「ラベル」にドラッグアンドドロップします。

4.ラベル＞フォント＞フォントサイズを24にします。

「距離」が表示されました。しかし、数値がおかしい気がします。蒲田駅～横須賀中央駅間は3049kmもありませんね・・・

5.合計(距離)＞メジャー(合計)＞平均を押下します。

以下のエントリーと同じく蒲田駅～横須賀中央駅間は31.76kmという結果になりました。

www.gis-py.com

バッファー

ポイントに対して300mのバッファーを作ります。

1.分析＞計算フィールドの作成

2.BUFFER([ポイント],300,"m")と入力＞OK

左側のデータ欄に「バッファー」が作成されたことがわかります。

3.「バッファー」をダブルクリック
バッファーが作成されたことがわかります。

さいごに

まさか Tableau でバッファーまで作れるとは思ってなかったので、これは意外な発見でした。まだまだ知らない機能がありそうなので、引き続き色々触って勉強してみようと思います。本日は以上です。

Tableau を使って地図データを可視化してみよう！

さて、本日は Tableau を使って地図データを可視化してみようと思います。Tableau と言えば BI ツールとして有名だと思いますが、なんと これで地図データも可視化できるようです。ちょっと地図データを可視化したいくらいだったら GIS を導入しなくても BI ツールで対応できてしまうということですね。おそらく地図関係の機能もどんどん強化されていくでしょうし、GIS 関係の会社さんは戦々恐々としているのではないでしょうか。

Tableau とは

Tableauは、世界的に人気な BI ツールです。誰でも簡単にデータ分析や可視化を行うことができます。

Tableau Desktop をインストールしてみる

Tableau は色々な製品があるのですが、今回はTableau Desktop Public Edition という無料で使える製品を使ってみます。以下のサイトからインストーラーをダウンロードします。

public.tableau.com

インストーラーを起動してインストールを進めます。

インストール完了後、Tableau を起動してみます。機能は少し制限されているそうなのですが、少しつかうくらいだったらこれで十分だと思います。

使用するデータ

• 全国市区町村界・・・ESRIジャパンさんが公開しているデータです。以下関連エントリーなので、興味がある方はぜひ読んでみてください。

www.gis-py.com

データを可視化してみる

空間ファイルを選択します。

全国市区町村界データを読み込むとこのような画面になります。「シート1」をクリックします。

このような表示になります。「ジオメトリ」をダブルクリックします。

おぉー全国市区町村界データが表示されました。背景には Mapbox の地図を使っていますね。

次に「KEN」を「ラベル」にドラッグアンドドロップします。

県名がラベリングされました。

同じように「KEN」を「色」にドラッグアンドドロップするとこのように県ごとの色分けが行われます。

こんな感じで属性テーブルも使うことができますね。

さいごに

やはり 地図機能に関しては GIS 製品と比べると見劣りするかなという感じでしたが、ちょっと地図データを可視化させてみたいという場合は Tableau で十分かなと思いました。位置情報に関する重要性が年々増してきているかと思うので、これから Tableau の地図機能はさらに充実してくるのではと思っています。私も仕事で Tableau を使うようになったのでこれを機に本ブログで Tableau について色々紹介できればと思います。本日は以上です。

Elasticsearch と Kibana を使って地図データを可視化してみよう！

さて、本日は Elasticsearch と Kibana を使って地図データを可視化してみようと思います。ただ、Elasticsearch と Kibana のインストールや使い方については本エントリーでは紹介はしませんので、興味のある方は公式ドキュメントを読んでみてください。

www.elastic.co

Elasticsearch とは

Elasticsearch は Apache Lucene をベースにしたオープンソースの検索エンジンです。近年では利用者が急増して、現在では最も人気のある検索エンジンとなっています。Elasticsearch は元々全文検索の用途で開発されましたが、それ以外にもログ収集・解析の用途としても多く利用されています。

Elasticsearch の特徴

• 分散配置による高速化と高可用性
• REST API によるアクセス
• JSON フォーマットに対応したドキュメント志向データベース
• ログ収集・可視化などの多様な関連ソフトウエアとの連携

Kibana とは

Kibana は Elasticsearch に格納されているデータをブラウザを用いて可視化、分析するためのツールです。グラフ、ダッシュボード、地図などを使ってデータを可視化することができます。

環境

CentOS 7
Elasticsearch7.6
Kibana7.6

手順

1. データ作成
2. データ可視化

データ作成

今回は2つデータを作成しました（駅名とそこの緯度経度を持つデータ）。curl コマンドか Kibana コンソールを使ってデータを作成できるので、両方のやり方を紹介します。

1.curl コマンド

curl -XPUT "http://localhost:9200/test_index" -H 'Content-Type: application/json' -d'{  "mappings": {    "properties": {      "location": {        "type": "geo_point"      }    }  }}'
curl -XPUT "http://localhost:9200/test_index/_doc/1" -H 'Content-Type: application/json' -d'{  "text": "蒲田駅",  "location": {     "lat": 35.562479,    "lon": 139.716073  }}'
curl -XPUT "http://localhost:9200/test_index/_doc/2" -H 'Content-Type: application/json' -d'{  "text": "横須賀中央駅",  "location": {     "lat": 35.278699,    "lon": 139.670040  }}'
curl -XGET "http://localhost:9200/test_index/_search"

2.Kibana コンソール

PUT test_index
{
  "mappings": {
    "properties": {
      "location": {
        "type": "geo_point"
      }
    }
  }
}

PUT test_index/_doc/1
{
  "text": "蒲田駅",
  "location": { 
    "lat": 35.562479,
    "lon": 139.716073
  }
}

PUT test_index/_doc/2
{
  "text": "横須賀中央駅",
  "location": { 
    "lat": 35.278699,
    "lon": 139.670040
  }
}

GET test_index/_search

このような形でデータが格納されていることがわかります。

{
  "took" : 1,
  "timed_out" : false,
  "_shards" : {
    "total" : 1,
    "successful" : 1,
    "skipped" : 0,
    "failed" : 0
  },
  "hits" : {
    "total" : {
      "value" : 2,
      "relation" : "eq"
    },
    "max_score" : 1.0,
    "hits" : [
      {
        "_index" : "test_index",
        "_type" : "_doc",
        "_id" : "1",
        "_score" : 1.0,
        "_source" : {
          "text" : "蒲田駅",
          "location" : {
            "lat" : 35.562479,
            "lon" : 139.716073
          }
        }
      },
      {
        "_index" : "test_index",
        "_type" : "_doc",
        "_id" : "2",
        "_score" : 1.0,
        "_source" : {
          "text" : "横須賀中央駅",
          "location" : {
            "lat" : 35.278699,
            "lon" : 139.67004
          }
        }
      }
    ]
  }
}

データ可視化

まずはKibana のマップを表示させてみました。続いて作成したデータを Kibana で可視化してみようと思います。

ポイントが二つプロットされていることがわかります。それぞれ拡大して確認してみようと思います。

蒲田駅と横須賀中央駅にそれぞれポイントがプロットされていることがわかります

さいごに

ものすごく簡単にでしたが Elasticsearch と Kibana を使って地図データの可視化に成功しました。データの表示だけではなくデータの分析もできるっぽいので、もう少し使い方を覚えたらそれについても紹介したいと思います。Elasticsearch に関しては今後ますます人気が高まるでしょうし、興味のある方はぜひ使ってみてください。本日は以上です。

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

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

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

距離計算方法

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

ヒュベニの公式

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

Ax：２点の経度の差
Ay：２点の緯度の差
P：２点の緯度の平均

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

e = SQRT((Rx² - Ry²)/Rx²)：第1離心率
Rx：長半径（赤道半径）
Ry：短半径（極半径）
E2 = e²：第2離心率

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

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

• 国土地理院方式
• 球面三角法
• 測地線航海算法

距離計測する区間

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

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

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

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

実行環境

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()

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

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

さいごに

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

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

参考サイト

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