アメダスのデータをAPIを利用して取得する

Tag

はじめに

気象庁の保有する全国各地の地域気象観測所(アメダス)で得られた1分値データです。

降水量、気温、風向、風速、日照時間、積雪の深さなどが含まれています。

本チュートリアルではTellusのJupyterLabを使って地域気象観測システム(アメダス)の気象データを取得する方法を紹介します。

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

TellusではAPIによりアメダス1分値(1分毎の観測データ)を取得することができます。APIの使い方はリファレンスにまとめられています。 ※データを利用する際は、利用ポリシーを確認して、規約を違反しないように注意してください。

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

Name Input Description
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を利用してアメダスのデータを取得する関数を定義しましょう。

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

サンプルコードでは、日本時間で2016年1月1日0時1分のデータを取得しています。一回のリクエストでその時刻における全国の観測所のデータが返ってきます。2016年1月の時点で観測所は1284件と示されています。 APIからは観測時刻(timestamp)、座標(経度緯度,coordinates)、観測所情報(station)、雨量計データ(rain)、風速風向計データ(wind)、気温計データ(temperature)、日照計データ(sunshine)、積雪計データ(snow)が返されます。各観測値は以下の構造となっています。

{
    "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箇所の観測所データが返ってきました。

雨量計データ

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

amedas2[:][0]['rain']

devAmedas_api1_20200220_2.png

intensityに降水強度、intensity_flagにその品質を示す値が示されています。降雨強度を値の大きい順に並び替え、上から10ヶ所分表示してみます。

取得する観測時刻を日本時間で2017年10月28日8時30分とし、値が初期値(2147483647)の場合や、利用フラグが12より大きい(観測品質が悪い)場合は除外値とします。

上記の時間帯で観測品質が悪いデータがどれ程含まれているのかを最初に確認してみましょう。

amedas2 = get_amedas_1min('2017', '10', '28', '8', '30')
len([x['rain']['intensity_flag']['value'] < 13 for x in amedas2])

1282が返されました。観測時刻のデータは1282箇所あったので、一つとして観測品質が悪かったデータがなかったことになります。また、観測された値に最大のバイト数が入っていた場合にも除外をします。除外対象となる値をそのバイトをキーとした辞書型のデータ(byte_size)とし、関数内での条件分け基準として使うことにしましょう。

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']]

上で定義した関数を実行します。観測されたデータは降雨強度ソートされ、条件にあった値のみが取得されたデータがから抽出されます。その中から10個のみのデータを表示します。 *取得したrainのデータはキャリブレーションされています。そのため、値を10で割り、元の単位に戻す必要があります。

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]')

devAmedas_api1_20200220_3.png

結果から、2017年10月28日8時30分における降雨強度の全国1位は観測所番号87412、宮崎県の赤江観測所でした。50mm/hを越えると豪雨と言われるので、112.5 mm/hはかなり強い雨が降っていることになります。2位以下も宮崎や熊本、鹿児島と南側に偏っています。これは九州地方に台風22号が接近していた影響です。

devAmedas_api1_20200220_1.png

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

Credit : 気象庁

風速風向計データ

風速風向計データには、最大瞬間風速(3秒移動平均)、最大瞬間風速(3秒移動平均)、平均風速(10分移動平均)、平均風向(前10分間のベクトル平均)などの観測値と、データの質を示すフラグが含まれています。ここでは、台風の接近に伴い平均風速(velocity_ave)と平均風向16方位(direction_ave_16direction)が変化していく様子を観察します。

データは2017年に台風22号が西日本に接近した期間(10月28日21時から24時間)を1時間毎に取得し、観測所は四国の南西に位置する清水特別地域気象観測所(観測所番号:74516)に絞り込みます。このとき値が初期値(2147483647)の場合や、利用フラグが12より大きい(観測品質が悪い)の場合は観測対象から除きます。
*取得したwindのデータはキャリブレーションされています。そのため、値を10で割り、元の単位に戻す必要があります。

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)

devAmedas_api1_20200220_5.png

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

気温計データ

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

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

気温を日本時間の2016年8月1日0時から24時にかけて1時間毎に取得し、値が初期値(2147483647)の場合や、利用フラグが12より大きい(観測品質が悪い)の場合はデータから除外します。
*取得したtemperatureのデータはキャリブレーションされています。そのため、値を10で割り、元の単位に戻す必要があります。

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

devAmedas_api1_20200220_4.png

Credit : 気象庁

グラフにしたことにより、1日の中でどのように気温が変化したのかを示すことができました。多治見観測所では僅か数時間で気温が急激に下がっています。多治見では15時以降に雨雲が発達し、日照量は低下、その後に強い雨が振り(約20mm程度)急激な温度変化をもたらしたと考えられます。雨は次第に弱くなりましたが、その後20時くらいまで振り続けたため、気温は23度程度安定したと読み取ることができます。

TellusのJupyterLabを使って地域気象観測システム(アメダス)のデータを取得する方法をご紹介しました。アメダス1分値は、1年分のデータだけでもレコード数が6億を越える膨大なデータです。高い時間分解能を活かし、衛星データの補完にお役立てください。

今回のスクリプトです。

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])
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))
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]')
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)
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
get_ipython().run_line_magic('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()