全観測地点・数十年分の気象データをPythonでスクレイピング

日本中の気温の推移を数十年単位で可視化してみたく、気象庁が公開している全ての観測地点の過去の気象データをスクレイピングしました。

以下のサイトを参考にしています。

https://www.gis-py.com/entry/scraping-weather-data

# 気象庁の過去データについて
気象庁の過去データ提供サービスは２種類あります。

本記事では２．過去の気象データ検索サービスの1日ごとのデータをスクレイピングします。

## 過去の気象データ・ダウンロードサービス
https://www.data.jma.go.jp/gmd/risk/obsdl/

こちらのサービスを使うと、特定の場所、特定の期間、特定の気象情報（気温、降水量、風速など）をcsv形式でダウンロードすることができます。

しかし一度にダウンロードできるデータ容量が限られているため、全観測地点のデータを長期間にわたって取得しようとすると、大変な回数ダウンロードしなければならず現実的ではありません。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2605106/3ee059f6-ceac-7b3d-fde2-02a3bfcce66a.png)

## 過去の気象データ検索サービス
https://www.data.jma.go.jp/obd/stats/etrn/index.php

こちらは、観測所、年月日を選択すると文字情報で気象情報を表示してくれるサイトです。10分毎、1時間ごと、1日ごとなど粒度別で取得が可能です。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2605106/b235dc56-fdad-281a-148e-6ef5b6f3b1f1.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2605106/70deff24-f153-c1be-0408-00bbf30fe504.png)
日ごとの値から、降水量の合計、平均気温・最高気温・最低気温、平均風速、降雪量の合計と最深積雪をスクレイピングで取得してみます。

## 地点コード一覧の取得
全観測地点のデータを取得するので、そのコード一覧が必要になります。

こちらのページに、観測地点のコードをまとめたobs_stations.xlsxというエクセルファイルがあるので、これをダウンロードしてください。

https://washitake.com/weather/amedas/obs_stations.md
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2605106/e4365b39-ed72-23dc-db78-fdbb531d306f.png)

# ライブラリ
使用したライブラリは以下の通りです。必要に応じてインストールしてください。
・pandas
・datetime
・dateutil.relativedelta
・urllib.request
・BeautifulSoup

# 実装
環境は以下を使用しました。
・Windows11
・Python3.9.15
・Jupyter notebook6.5.2

Jupyter notebook の形式で git にも保管しています。

https://github.com/BB-engineer/sandbox/tree/main/Weather


※アクセス集中を防ぐため1月取得するごとに５秒待つ処理を追加しました。
     
　　

パッケージのインポート
```python
import pandas as pd
import datetime
from dateutil.relativedelta import relativedelta
import urllib.request
from bs4 import BeautifulSoup
import time
```
　　
　　
観測地点一覧ファイルの読み込み。

一覧には観測を停止している地点も入っているので、観測を継続している地点のみを抽出します。また気温の観測を行っている地点に絞ります。
```python
obs_stations = pd.read_excel("obs_stations.xlsx")
obs_stations = obs_stations.query('ed_y == 9999')
obs_stations = obs_stations.query('気温 == "Y"')
```
　　

データ取得開始と終了年月日の設定。

好きな年月日を設定して下さい。

長くしすぎるとデータの取得に数日単位でかかってしまうので注意してください。
```python
start_date = datetime.date(2020, 1, 1)
end_date   = datetime.date(2020, 12, 31)
```
　　
　　
メイン部分
```python
Temp_data = pd.DataFrame(index=[], columns=['降水量','降雪量','積雪量','平均気温','最高気温','最低気温','平均風速'])

date = start_date
while date < end_date:
    Temp_data.loc[date] = -100.0
    date += relativedelta(days=1)

def str2float(weather_data):
    try:
        return float(weather_data)
    except:
        return -100

for i in obs_stations.index:
    Temp_data_s = Temp_data.copy()

    url_y = "https://www.data.jma.go.jp/obd/stats/etrn/view/annually_%s.php?" \
            "prec_no=%d&block_no=%04d&year=&month=&day=&view=" \
                %(str.lower(obs_stations['区分'].loc[i]), obs_stations['府県番号'].loc[i], obs_stations['地点コード'].loc[i])

    html = urllib.request.urlopen(url_y).read()
    soup = BeautifulSoup(html)
    trs = soup.find("table", { "class" : "data2_s" })

    if trs is None:
        continue
    
    tr = trs.findAll('tr')[3]
    tds = tr.findAll('td')[0].findAll('div')[0].findAll('a')

    date = start_date
    # date = datetime.date(max(int(tds[0].string),1872), 1, 1)
    while date < end_date:
        url_m = "https://www.data.jma.go.jp/obd/stats/etrn/view/daily_%s1.php?" \
                "prec_no=%s&block_no=%04d&year=%d&month=%d&day=&view=" \
                    %(str.lower(obs_stations['区分'].loc[i]), obs_stations['府県番号'].loc[i], obs_stations['地点コード'].loc[i], date.year, date.month)

        html = urllib.request.urlopen(url_m).read()
        soup = BeautifulSoup(html)
        trs = soup.find("table", { "class" : "data2_s" })

        if trs is None:
            print(f"\r{obs_stations['地点'].loc[i]}", end="")
            date += relativedelta(months=1)
            continue

        # table の中身を取得
        date_day = date
        print(f"\r{i,obs_stations['地点'].loc[i], obs_stations['地点コード'].loc[i], date_day.year, date_day.month}", end="")
        
        if obs_stations['区分'].loc[i] == 'S':
            for tr in trs.findAll('tr')[4:]:
                tds = tr.findAll('td')
                
                if tds[6].string == None or tds[7].string == None or tds[8].string == None:
                    break
                Temp_data_s['降水量'].loc[date_day] = str2float((tds[3].string).split(' ')[0])
                Temp_data_s['平均気温'].loc[date_day] = str2float((tds[6].string).split(' ')[0])
                Temp_data_s['最高気温'].loc[date_day] = str2float((tds[7].string).split(' ')[0])
                Temp_data_s['最低気温'].loc[date_day] = str2float((tds[8].string).split(' ')[0])
                Temp_data_s['降雪量'].loc[date_day] = str2float((tds[17].string).split(' ')[0])
                Temp_data_s['積雪量'].loc[date_day] = str2float((tds[18].string).split(' ')[0])
                Temp_data_s['平均風速'].loc[date_day] = str2float((tds[11].string).split(' ')[0])
                date_day += relativedelta(days=1)
        else:
            for tr in trs.findAll('tr')[3:]:
                tds = tr.findAll('td')
 
                if tds[4].string == None or tds[5].string == None or tds[6].string == None:
                    break
                Temp_data_s['降水量'].loc[date_day] = str2float((tds[1].string).split(' ')[0])
                Temp_data_s['平均気温'].loc[date_day] = str2float((tds[4].string).split(' ')[0])
                Temp_data_s['最高気温'].loc[date_day] = str2float((tds[5].string).split(' ')[0])
                Temp_data_s['最低気温'].loc[date_day] = str2float((tds[6].string).split(' ')[0])
                Temp_data_s['降雪量'].loc[date_day] = str2float((tds[16].string).split(' ')[0])
                Temp_data_s['積雪量'].loc[date_day] = str2float((tds[17].string).split(' ')[0])
                Temp_data_s['平均風速'].loc[date_day] = str2float((tds[9].string).split(' ')[0])
                date_day += relativedelta(days=1)

        date += relativedelta(months=1)

        time.sleep(5)

    Temp_data_s.to_csv("%d_%s.csv"%(obs_stations['地点コード'].loc[i],obs_stations['地点'].loc[i]))
```
　　
　　
なお上記のコードでは、欠測や降水がない場合などで数値が入っていないところを一律 -100 にしています。

実際には以下の記号が入っているので区別したい場合はカスタマイズしてください。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2605106/e68e0324-70ce-f026-8697-533f57457485.png)
# 結果
各観測地点の気象データがそれぞれのファイルで保存されました。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2605106/ae223967-412c-207c-3fb1-0a07016748b9.png)
中身はこんな感じです。ちゃんと値を取得できていますね。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2605106/6308c1a1-5697-a653-f181-679fcde0e6d7.png)
試しに 2020/1/1 の最高気温を地図上にプロットしてみました。

各観測地点の緯度経度は準備でダウンロードした obs_stations.xlsx に記載されています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2605106/2e38a3c5-e0d0-5b2b-3eb3-c0671768cfd3.png)

# まとめ
今回は、Python を使って、気象庁の過去の気象データ検索サービスから、特定期間の全ての観測地点の気象データをスクレイピングする方法について記載しました。