自分の家から一番近い気象台を調べよう！ - geopandasで位置情報分析 -
はじめに
天気予報がどのように発表されているかご存知ですか？
何気なく目にする天気予報は、さまざまなプロセスを経て我々の元に届きます。
[天気予報が発表されるまでのプロセス]
はじめに、気象台より観測データの収集が行われます。その後、数値予報モデルを実行し、大気や海洋・陸地の状態の変化をシミュレーションします。シミュレーションされた情報を元に予報官により天気予報の作成・発表が行われます。
今回は、気象台がどのように配置されているのか。また、自分の家の近くの気象台はどこなのかを調べていきます。
気象台について詳しく
気象台は天気予報や災害情報を収集・発信する施設です。日本には5つの管区気象台（札幌、仙台、東京、大阪、福岡）があり、その下に多数の地方気象台や測候所が配置されています。
[気象台の組織図]
高層気象台は大気の上層部における気象現象を観測しています。これは、雨や風、雪や雷などの天気現象は、はるか上空の大気の層の中で起きているためです。南極の昭和基地も高層気象台に分類されます。
geopandasを使って気象台をプロットしよう！
主なパッケージ
- pandas
- geopandas
- matplotlib
- shapely
地上観測所の情報をGeoDataFrameに格納
気象庁のページから気象観測所の一覧（地上気象観測 地点情報履歴ファイル）をダウンロードします。
地上気象観測フォーマットを参考に、ダウンロードしたファイルから地点番号や地点名などの要素を抽出しています。
ちなみに、観測所区分（OBSERVATORY_CODES）に「旧地台」や「新地台」という用語があります。これらはおそらく「旧・地方気象台」「新・地方気象台」を表していると考えられます。
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
KEEP_OBSERVING_CODE = '99999999'
OBSERVATORY_CODES = {
"1": "管区",
"2": "旧地台",
"3": "海洋気象台(2013 年 9 月 30 日以前)",
"4": "新地台",
"5": "測候所",
"6": "地球環境業務課",
"7": "施設等機関",
"8": "特別地域気象観測所"
}
def dms_to_decimal(dms, is_longitude):
"""度分秒を緯度経度に変換"""
if is_longitude:
# 経度の場合は3桁の度数を取得
degrees = int(dms[:3])
minutes = int(dms[3:5])
seconds = int(dms[5:7])
else:
# 緯度の場合は2桁の度数を取得
degrees = int(dms[:2])
minutes = int(dms[2:4])
seconds = int(dms[4:6])
return degrees + minutes / 60 + seconds / 3600
def observatory_code_to_site(code):
try:
return OBSERVATORY_CODES[code]
except KeyError:
raise ValueError(f"Invalid observatory code: {code}")
file_path = "./data/smaster.index"
with open(file_path, encoding="shift_jis") as f:
content = f.readlines()
# indexファイルから観測所情報を抽出
site_data = []
for line in content:
site_name = line[80:93].split(" ")[0]
bureau_name = line[80:93].split(" ")[-1]
observatory = line[14:15]
observation_period_end =line[72:80]
# 現在も使われている観測所のみ格納
if observation_period_end == KEEP_OBSERVING_CODE:
record = {
"地点番号": line[0:4],
"地点名": site_name,
"振興局": bureau_name,
"観測所区分": observatory_code_to_site(line[14:15]),
"カナ地点名": line[16:24].split(" ")[0],# 半角空白削除
"緯度": dms_to_decimal(line[36:42],is_longitude=False),
"経度": dms_to_decimal(line[42:50],is_longitude=True)
}
site_data.append(record)
df = pd.DataFrame(site_data)
geometry = gpd.points_from_xy(df["経度"],df["緯度"])
gdf = gpd.GeoDataFrame(df, geometry=geometry,crs="EPSG:4326")
gdf = gdf.loc[gdf["地点名"] != "昭和"] # 昭和基地の観測所を省く
gdfに観測所の位置情報を含んだデータが格納されています。
南極にある昭和基地観測所は、日本地図が小さく表示されてしまうため、取り除いてあります。
gdfの中身は以下のようになります。geometryカラムに観測所の位置が格納されています。
日本地図と重ねてプロットしてみます。
import geopandas as gpd
prefectures_file_path = "./weather/data/prefectures.geojson"
japan = gpd.read_file(prefectures_file_path)
ax = japan.plot(facecolor='lightgray')
gdf.plot(ax=ax, color='forestgreen',markersize=5)
日本全国に観測所が配置されているのがわかりますね。
蛇足：南鳥島
ちなみに右下にある、日本の最東端の島、南鳥島でも人が住んで365日24時間観測しているそうです。約1.5㎢の隔離された島に住むのはなかなかハードな気がしますが、一面海しかない環境だと邪念がなく作業できそうですね。
三角の島で可愛いですね。
家から最も近い観測所を表示
サンプルとして東京スカイツリーから最も近い観測所を抽出しました。
from shapely.ops import nearest_points
from shapely.geometry import Point
# 東京スカイツリーの緯度経度
target_point = Point(139.8140455 ,35.7067694)
nearest =nearest_points(target_point,gdf.geometry.unary_union)
nearest_x, nearest_y = nearest[1].x, nearest[1].y
# プロット
ax = japan.plot(facecolor='lightgray',edgecolor='white')
ax.scatter(target_point.x, target_point.y, color="red", s=50,label="Target Point")
ax.scatter(nearest[1].x, nearest[1].y, color="green", s=50,label="Nearest Point")
# 拡大表示
zoom_margin = 0.5
ax.set_xlim(nearest_x - zoom_margin, nearest_x + zoom_margin)
ax.set_ylim(nearest_y - zoom_margin, nearest_y + zoom_margin)
ax.legend()
データフレーム情報は以下で表示できます。
nearest_record = gdf[gdf.geometry == nearest[1]]
nearest_record
おわりに
geopandasを使って気象台の位置情報を可視化しました。ちなみに数値予想モデルにはスーパーコンピュータ「富岳」が使われているみたいですよ。
参考
Discussion