宙畑 Sorabatake

ニュース

【ゼロからのTellusの使い方】Jupyter Labで地域気象観測システム(アメダス)のデータを取得する

Tellusでは地域気象観測システム(アメダス)の雨量計、風速風向計、気温計、日照計、積雪計の1分毎の観測値を利用することができます。 本記事ではTellusのJupyter Labを使ってアメダスのデータを取得する方法を紹介します。

本記事ではTellusのJupyter Labを使って地域気象観測システム(アメダス)のデータを取得する方法を紹介します。

TellusでJupyter Lab(Jupyter Notebook)を使う方法は、「Tellusの開発環境でAPI引っ張ってみた」をご覧ください。

1. APIでアメダスのデータを取得する

TellusではAPIによりアメダス1分値(1分毎の観測データ)を取得することができます。

APIの使い方はリファレンスにまとめられています。
APIリファレンス
※データを利用する際は、データ詳細のデータポリシーを確認して、規約を違反しないように注意してください。

https://gisapi.tellusxdp.com/api/v1/amedas/1min/{yyyy}/{MM}/{dd}/{hh}/{mm}/{?min_lat,min_lon,max_lat,max_lon}

引数

yyyy 必須
str
観測時刻の西暦年 (4桁の数字で指定) 
※世界標準時
例:’2016’
MM 必須
str
観測時刻の月 01~12 (2桁の数字で指定)
※世界標準時
例:’01’
dd 必須
str
観測時刻の日 01~31 (2桁の数字で指定)
※世界標準時
例:’01’
hh 必須
str
観測時刻の時 00~23 (2桁の数字で指定)
※世界標準時
例:’12’
mm 必須
str
観測時刻の分 00~59 (2桁の数字で指定) 
※世界標準時
例:’00’
min_lat 任意
number
観測所の最小緯度(-90.0〜90.0)
min_lon 任意
number
観測所の最小経度(-180.0〜180.0)
max_lat 任意
number
観測所の最大緯度(-90.0〜90.0)
max_lon 任意
number
観測所の最大経度(-180.0〜180.0)

※現在Tellusで利用できるのは、世界標準時で2015年12月31日15:01から2017年12月31日15:00のデータとなります。

 

このAPIをJupyter LabのNotebookから呼び出してみましょう。

Tellus IDEを開き、「work」ディレクトリに移動してから「Python3」を選択し、以下のコードを貼り付けます。

「work」ディレクトリに移動(Jupyter Labの画面例)
「Python3」のパネルを選択

import requests, json

TOKEN = 'ここには自分のアカウントのトークンを貼り付ける'

def get_amedas_1min(year, month, day, hour, minute, payload={}):
    """
    /api/v1/amedas/1minを叩く

    Parameters
    ----------
    year : string
        西暦年 (4桁の数字で指定) UTC
    month : string
        月 01~12 (2桁の数字で指定) UTC
    day : string
        日 01~31 (2桁の数字で指定) UTC
    hour : string
        時 00~23 (2桁の数字で指定) UTC
    minute : string
        分 00~59 (2桁の数字で指定) UTC
    payload : dict
        パラメータ(min_lat,min_lon,max_lat,max_lon)
    
    Returns
    -------
    content : list
        結果
    """

    url = 'https://gisapi.tellusxdp.com/api/v1/amedas/1min/{}/{}/{}/{}/{}/'.format(year, month, day, hour, minute)
    
    headers = {
        'Authorization': 'Bearer ' + TOKEN
    }
    r = requests.get(url, headers=headers, params=payload)

    if r.status_code is not 200:
        print(r.content)
        raise ValueError('status error({}).'.format(r.status_code))
    return json.loads(r.content)  

amedas1 = get_amedas_1min('2015', '12', '31', '15', '01')
print(len(amedas1))
print(amedas1[0])

トークンはマイページのAPIアクセス設定(要ログイン)から取得してください。

上部の三角形アイコンをクリックしてコードを実行すると、指定した時刻のデータが返ってきます。

サンプルコードでは、日本時間で2016年1月1日0時1分のデータを取得しています。
一回のリクエストでその時刻における全国の観測所のデータが返ってきてます。
2016年1月の時点で観測所は1284件あるようです。

APIからは観測時刻、座標(経度緯度)、観測所情報、雨量計データ、風速風向計データ、気温計データ、日照計データ、積雪計データが返されます。
各観測値は以下の構造となっており、

{
    "value": int,
    "size": int
}

センサからデータを取得できなかった観測においては、値(value)にビット長(size)の最大値が与えられています。(例:{“size”: 4, “value”: 2147483647})

また、緯度経度で取得する範囲を絞り込むこともできます。

payload = {
    'min_lat': 39.0,
    'min_lon': 139.0,
    'max_lat': 40.0,
    'max_lon': 140.0
}
amedas2 = get_amedas_1min('2015', '12', '31', '15', '01', payload)
print(len(amedas2))

サンプルでは、経度が139 から 140°、緯度が39 から 40°の範囲で絞ってリクエストを投げた結果、該当する5件の観測所データが返ってきました。

2. 雨量計データ

雨量計データには、前1分降水量、降水強度、最大降水強度などの観測値と、データの質を示すフラグが含まれています。

ここでは、降雨強度(intensity)を値の大きい順に並び替え、上から10ヶ所分表示してみます。

取得する観測時刻を日本時間で2017年10月28日17時30分とし、値が初期値(2147483647)の場合や、利用フラグが12以上(観測品質が悪い)の場合は、データをエラー値として今回の例では対象から除きます。

def is_valid(value_obj, flag_obj):
    byte_size = {
        1: 127,
        2: 32767,
        4: 2147483647
    }
    return flag_obj['value'] <= 12 and value_obj['value'] != byte_size[value_obj['size']]

amedas2 = get_amedas_1min('2017', '10', '28', '8', '30')
amedas2_sorted = sorted(amedas2, key=lambda d: d['rain']['intensity']['value'], reverse=True)
amedas2_filtered = list(filter(lambda d: is_valid(d['rain']['intensity'], d['rain']['intensity_flag']), amedas2_sorted))

for i in range(10):
    print('station_id : ' + str(amedas2_filtered[i]['station']['prefecture_no']['value']).zfill(2) + str(amedas2_filtered[i]['station']['observatory_no']['value']).zfill(3)
          + ', intensity : ' + str(amedas2_filtered[i]['rain']['intensity']['value']/10) + '[mm/h]')

降雨強度の全国1位は観測所番号87412、宮崎県の赤江観測所でした。
50mm/hを越えると豪雨と言われるので、112.5 mm/hはかなり強い雨が降っていることになります。

2位以下も宮崎や熊本、鹿児島の観測所ばかりですが、これは九州地方に台風22号が接近していた影響です。

Tellus OS上で2017年10月28日のひまわり8号の「赤外1(B13)」を選択肢、左下に台風の渦が確認できる Source : 気象庁

3. 風速風向計データ

風速風向計データには、最大瞬間風速(3秒移動平均)、最大瞬間風速(3秒移動平均)、平均風速(10分移動平均)、平均風向(前10分間のベクトル平均)などの観測値と、データの質を示すフラグが含まれています。

ここでは、台風の接近に伴い平均風速(velocity_ave)と平均風向16方位(direction_ave_16direction)が変化していく様子を観察します。

データは2017年に台風22号が西日本に接近した期間(10月28日21時から24時間)を1時間毎に取得し、観測所は四国の南西に位置する清水特別地域気象観測所(観測所番号:74516)に絞り込みます。
またこのとき値が初期値(2147483647)の場合や、利用フラグが12以上(観測品質が悪い)の場合は、データをエラー値として今回の例では対象から除きます。

from datetime import datetime
from datetime import timedelta

direction_label = ['無風', 
                   '北北東', '北東', '東北東', '東', 
                   '東南東', '南東', '南南東', '南',
                   '南南西', '南西', '西南西', '西', 
                   '西北西', '北西', '北北西', '北',]

time = datetime(2017,10,28,21,0,0) - timedelta(hours=9)
for i in range(24):
    amedas_3 = get_amedas_1min(str(time.year), str(time.month).zfill(2), str(time.day).zfill(2), str(time.hour).zfill(2), str(time.minute).zfill(2))
    amedas_91011 = next(d for d in amedas_3 if str(d['station']['prefecture_no']['value']).zfill(2) + str(d['station']['observatory_no']['value']).zfill(3) == '74516')
    
    if is_valid(amedas_91011['wind']['velocity_ave'], amedas_91011['wind']['velocity_ave_flag']):
        print(amedas_91011['timestamp'], '{}[m/s]'.format(amedas_91011['wind']['velocity_ave']['value']/10), direction_label[amedas_91011['wind']['direction_ave_16direction']['value']])
    
    time += timedelta(hours=1)

台風22号は四国の南の沖を西から東へ抜けていきましたが、観測結果からも台風の接近に合わせて風向きが東、北、西と、台風の中心方向へ変化していく様子が伺えます。

4. 気温計データ

気温計データには、気温、最高気温 (前1分間)、最低気温 (前1分間)などの観測値と、データの質を示すフラグが含まれています。

ここでは、最高気温記録の上位にランクインしている熊谷地方気象台(観測所番号:43056)、江川崎地域気象観測所(観測所番号:74381)、多治見地域気象観測所(観測所番号:52606)の気温をグラフにして重ねてみます。

データは日本時間の2016年8月1日0時から24時にかけて1時間毎に取得し、値が初期値(2147483647)の場合や、利用フラグが12以上(観測品質が悪い)の場合は、データをエラー値として今回は対象から除きます。

import matplotlib.pyplot as plt
import matplotlib.dates as mdates
%matplotlib inline

time = datetime(2016,8,1,0,0,0) - timedelta(hours=9)

kumagaya = []
egawa = []
tajimi = []
periods = []

for i in range(24):
    amedas_4 = get_amedas_1min(str(time.year), str(time.month).zfill(2), str(time.day).zfill(2), str(time.hour).zfill(2), str(time.minute).zfill(2))
    # 熊谷
    amedas_43056 = next(d for d in amedas_4 if str(d['station']['prefecture_no']['value']).zfill(2) + str(d['station']['observatory_no']['value']).zfill(3) == '43056')
    # 江川崎
    amedas_74381 = next(d for d in amedas_4 if str(d['station']['prefecture_no']['value']).zfill(2) + str(d['station']['observatory_no']['value']).zfill(3) == '74381')
    # 多治見
    amedas_52606 = next(d for d in amedas_4 if str(d['station']['prefecture_no']['value']).zfill(2) + str(d['station']['observatory_no']['value']).zfill(3) == '52606')
    
    if is_valid(amedas_43056['temperature']['deg'], amedas_43056['temperature']['deg_flag']):
        kumagaya.append(amedas_43056['temperature']['deg']['value']/10)
    else:
        kumagaya.append(0)
    if is_valid(amedas_74381['temperature']['deg'], amedas_74381['temperature']['deg_flag']):
        egawa.append(amedas_74381['temperature']['deg']['value']/10)
    else:
        egawa.append(0)
    if is_valid(amedas_52606['temperature']['deg'], amedas_52606['temperature']['deg_flag']):
        tajimi.append(amedas_52606['temperature']['deg']['value']/10)
    else:
        tajimi.append(0)
    periods.append(time + timedelta(hours=9))
    time += timedelta(hours=1)

fig, ax= plt.subplots(figsize=(12, 6))
ax.plot(periods, kumagaya, label = 'no marker')
ax.plot(periods, egawa, label = 'no marker')
ax.plot(periods, tajimi, label = 'no marker')
ax.legend(['kumagaya', 'egawa', 'tajimi'])
xfmt = mdates.DateFormatter('%H')
xloc = mdates.HourLocator()
ax.xaxis.set_major_locator(xloc)
ax.xaxis.set_major_formatter(xfmt)
ax.set_xlabel('hour(JST)')
ax.set_ylabel('temp')
ax.set_xlim(periods[0], periods[-1]) 
plt.show()

この日は3箇所のうち、多治見が最も暑かったようです。

以上が、TellusのJupyter Labを使って地域気象観測システム(アメダス)のデータを取得する方法をご紹介しました。

アメダス1分値は、1年分のデータだけでもレコード数が6億を越える膨大なデータです。高い時間分解能を活かし、衛星データの補完にお役立てください。