Pythonのみで業務スーパーの各店舗と気象状況を可視化する

みなさん、業務スーパー使ってますか？私は、朝ごはんに業務スーパーのインスタントフォー（チキン味）を食べるようにしていることから、週1回は業務スーパーに通っています。最近、全都道府県に店舗が出店されたそうです（パチパチパチ）。

というわけで、業務スーパーの全店舗の位置をプロットしてみようかと思いました。あと可視化だけだとちょっと退屈なので、気象データを使っていくつかの店舗の気温の比較を行います。foliumで作った可視化は次のようになります。

開発環境:

Python 3.9.1
folium 0.12.1
pandas 1.2.3
requests 2.25.1
requests_html 0.10.0
xarray 0.16.2

ハンズオン開催します

はんなりPythonで、この記事であったrequests / requests-htmlを使ったデータ取得ハンズオンを行います。ご関心をお持ちの方はconnpassのサイトからお申し込みください。

利用データ

業務スーパーの位置データ

これは業務スーパーのサイトから取得しました。requests_htmlを使いました。データを集める際は便利なモジュールです。

サイトからデータを取得する際の注意としては次のサイトに詳しいので参照ください。

天候のデータ

天候のデータは気象庁の気象過去データで得たデータを利用します。様々なデータを取得できるデータの中から、2キロ格子間隔でデータが取得できる局地モデルのデータを利用します。これにより、全国の2km格子の気象データを使うことができ、アメダスよりも詳細な気象データを使うことができます。

位置データの取得

住所の取得

各店舗の住所のデータは業務スーパーのサイトから取得します。利用したパッケージはrequests_htmlです。requests_htmlはちょっとデータがとり難いサイトもレンダーさせてから取得出来るなど、便利なパッケージです。例えば、次のような感じで使えます。

from requests_html import HTMLSession

session = HTMLSession()
r = session.get('https://hogehoge.jp')
r.html.render() # サイトを描画させる

hoge_class = r.html.find('.hoge') # hoge クラスのデータを取得

上のような感じで簡単に使えるrequests_htmlを使って、908店舗分の住所と店名を集めました。

経度緯度の取得

住所を集めただけでは、位置情報を可視化することはできません。というか、コンピュータ上では住所はあまり役に立ちません。

というわけで、住所を経度緯度にしたデータを作成します。経度緯度の作成にはYahoo!のジオコーダAPIを使いました。このAPIはアプリケーションIDを取得すると使うことができます。

Yahoo!のサイトには次のようにリクエストパラメータ一覧があります。これのappidにアプリケーションIDを、queryに住所文字列を渡すと位置情報を返してくれます。

これで簡単という感じに思いますが、実際は位置情報が返されない住所があります。今回の場合90件くらいがそうなりました。そのような際にどうするか・・・というのはちょっと良いアイデアが浮かばなかったので、1文字ずつ文字列を削って、位置情報が見つかるまでというものにしました。コード自体は次のような感じです。

import time
import xml.etree.ElementTree as ET
import pandas as pd
import requests

def req_latlng(base_url: str, shop_address: str) -> str:
    '''
    apiの返り値が無効であれば
    住所から一文字削除してもう一度試す
    再帰関数
    shop_address str: 住所の文字列
    '''
    req_url = base_url + f'?appid={appid}&query={shop_address}'
    r = requests.get(req_url)
    root = ET.fromstring(r.content)
    time.sleep(5)
    try:
        print(shop_address)
        latlng = root[1][3][1].text 
        return latlng
    except IndexError:
        return req_latlng(shop_address[:-1])

再帰関数を作って文字を一文字ずつ削って、経度緯度があるか試すようなものとなっています。

位置情報の可視化

位置情報の可視化にはfoliumを使いました。foliumは地名などが入った地図の可視化をidなどを取得せずに可視化できる優れものです。そして、使用方法は容易ながら凄いものが作れます。

まずはシンプルに全店舗をプロット

まずはシンプルにすべての位置データをプロットします。データはcsvファイルに入れており、pandasを使ってデータを読み込みfoliumのMarkerを用いて可視化します。

import folium
import pandas as pd

data_path = 'hoge.csv'
data = pd.read_csv(data_path)

mean_y = data['y'].mean() # データの中心点を取る（x,y）
mean_x = data['x'].mean()
m = folium.Map([mean_x, mean_y], zoom_start=5) #Mapオブジェクトを作る
for i in range(len(data)): # データを渡しプロットする
    folium.Marker([data.loc[i, 'x'], data.loc[i, 'y']], popup=f'店舗名: {data.loc[i, "shop_name"]}').add_to(m) 

結果は次のようなものとなります。

拡大してもすごい・・・クリックして店舗名が出せるのですが、関東に土地勘のない私には、この店舗が何県のものか分かりません。これではせっかく詳細な地図が見えるのに役に立ってない感が・・・

mean_y = data['y'].mean()
mean_x = data['x'].mean()
m = folium.Map([mean_x, mean_y], zoom_start=5)
marker_cluster = MarkerCluster().add_to(m) # 1.付け加え

for i in range(len(data)):
    folium.Marker([data.loc[i, 'x'], data.loc[i, 'y']], popup=f'店舗名: {data.loc[i, "shop_name"]}').add_to(marker_cluster) # 2.add_to()に渡すのをm => marker_clusterに変更

拡大すると次のような感じになります。地名の情報が見えるようになったので、エリアを意識してみるとかが簡単にできるようになりました。すげー

気温の比較

可視化するだけだと面白くないので、気象庁の局地モデルの気温のデータを使って店舗ごとの気温を比較してみます。局地モデルには、気温以外にも湿度、気圧、風、雲量なんかのデータが含まれます。

データ自体はgrib2形式という、気象データによくあるデータ形式で配信されています。データはpygribというモジュールを使って読み込めるのですが、xarrayを使いたいので、netcdf4形式に変換してから読み込み、データの取得を行います。

grib2 => netcfd4

ファイルの変換にはNOAAが作成しているwgrib2というソフトを用いました。導入するファイルの変換は簡単です。

$ wgrib2 grib2-path -netcdf4 netcdf4-path

pythonから使う場合はsubprocessモジュールを使います。

import subprocess 

subprocess.run('wgrib2 grib2-path -netcdf4 netcdf4-path')

xarrayを使ってデータを読み込み指定した位置のデータを取得する

作成したファイルを読み込み、指定した位置のデータを取得していきます。今回は12月30日9時から31日11時までの1時間ごとのデータを使って、気温の変化を観察します。

xarrayを使ったデータの読み込みですが、次のようにファイルをまとめて読み込めます。

import xarray as xr

data_paths = 'lfm_nc/2020123*.nc'
ds = xr.open_mfdataset(data_paths)

読み込んだデータセットは次のように表示されます。Coordnatesをみると緯度、経度、時間の次元であること、Data Variablesをみると気圧、風、気温などのデータが含まれることが分かります。

xarrayではpandasのようにラベル名で、データを切り出すことができます。ilocに対応するのがisel、locに対応するのがselです。500番目から600番目までの経度緯度の気温データを切り抜きたい場合は次のような感じで使えます。

ds.isel(latitude=slice(500, 600), longitude=slice(500, 600))['TMP_1D5maboveground']

北海道と沖縄の店舗の温度比較

さて最後に北海道と沖縄の店舗の温度比較を行いたいと思います。xarrayを使うとpandasのデータフレームに簡単に変換することができ、可視化も簡単に行えます。可視化はplotly.expressを用いて行います。

北海道の店舗は「アツベツニシテン」、沖縄の店舗は「ハエバルテン」とします。それぞれのデータを取得して可視化したのが次のグラフです。

これを見ると、沖縄の最低気温は11度、一方で北海道は―18度とかなり気温差は大きく、それゆえ商品展開とか色々あるんだろうなぁ・・・。インスタントフォーは沖縄や北海道には売っているのだろうか・・・と思いをはせられます。

同じ地方の店舗比較

近い地域でも気温を見てみたら結構違いはあるのだろうかと思い、九州の南北の店舗で確認してみました。

利用したのは「ハカタセンショウテン」と「ウエアラタテン」です。位置関係的には次のような感じです。

気温の可視化は次の通りです。予想通り北の方が気温は低いですね。この程度の気温差だとどんな感じで売れ行きに差が出るのか？ちょっと知りたいですね。

まとめ

今回の記事では業務スーパーの店舗をプロットして、ちょっと気温を見てみました。

参考にさせていただきた記事は次の4つです。ありがとうございます。

気象データを探していた時に過去データというものがあり、それでかなりな詳細なデータが取れるということを知ったのは次の記事でした。

wgrib2のインストール

xarrayの解説記事