import os
import glob
import re
import numpy as np
import pandas as pd
from datetime import datetime, timedelta
==========================================
1. 設定・地点データ定義
==========================================
アプリAのNPZデータが保存されているディレクトリ
CACHE_DIR = “gpv_cache_npz”
OUTPUT_EXCEL = “MSM_Sapporo_Local_Forecast.xlsx”
抽出対象の地点リスト（Location, Latitude, Longitude）
LOCATIONS = [
(“中央区土木センター”, 43.0715, 141.3184), (“北区土木センター”, 43.1257, 141.3444),
(“東区土木センター”, 43.1065, 141.3732), (“白石区土木センター”, 43.0425, 141.4239),
(“厚別区土木センター”, 43.0298, 141.4682), (“豊平区土木センター”, 43.0229, 141.4012),
(“南区土木センター”, 43.0134, 141.3397), (“西区土木センター”, 43.0763, 141.2722),
(“手稲区土木センター”, 43.1118, 141.2559), (“南16条西4丁目道路用地”, 43.0401, 141.3535),
(“円山公園坂下野球場”, 43.0526, 141.3148), (“盤渓配水池”, 43.0315, 141.2861),
(“北26条通り西8丁目”, 43.0905, 141.3417), (“新琴似グリーン公園”, 43.1165, 141.3211),
(“拓北水再生プラザ”, 43.1585, 141.3852), (“茨戸東部中継ポンプ場”, 43.1561, 141.3725),
(“農業支援センター”, 43.1154, 141.4429), (“札苗公園”, 43.1054, 141.4137),
(“一条大橋橋台広場”, 43.0592, 141.3653), (“豊平川水再生プラザ”, 43.0560, 141.3831),
(“厚別水再生プラザ”, 43.0255, 141.4617), (“大谷地流通団地東緑地”, 43.0238, 141.4552),
(“平岡公園”, 43.0062, 141.4632), (“清田区土木センター”, 43.0084, 141.4398),
(“西岡公園”, 42.9961, 141.3934), (“里塚霊園”, 42.9818, 141.4568),
(“中の沢雨水調整池”, 42.9933, 141.3175), (“石山1条4丁目緑地”, 42.9649, 141.3409),
(“藤野公園”, 42.9550, 141.3228), (“定山渓水再生プラザ”, 42.9734, 141.1718),
(“滝野自然学園”, 42.9161, 141.3957), (“福井中央公園”, 43.0416, 141.2676),
(“平和霊園”, 43.0317, 141.2585), (“手稲水再生プラザ”, 43.1257, 141.2335),
(“花川汚水中継ポンプ場”, 43.1678, 141.3155), (“高岡小中学校跡地”, 43.2081, 141.4239),
(“石狩斎場”, 43.2185, 141.3204), (“茨戸水再生プラザ”, 43.1612, 141.3456),
(“太美町汚水処理センター”, 43.1764, 141.4468), (“当別町下水終末処理場”, 43.2182, 141.5178),
(“金沢会館”, 43.2505, 141.5645), (“榎本公園”, 43.1091, 141.5312),
(“上江別南町公園”, 43.0963, 141.5547), (“大麻東公園”, 43.0805, 141.5037),
(“野幌農村環境改善センター”, 43.0645, 141.5204), (“北広島市白樺プール”, 42.9732, 141.5658),
(“島松寿町会館”, 42.9026, 141.5724), (“西部小学校”, 42.8557, 141.5303),
(“西の里”, 43.0232, 141.5199), (“東8丁目アンダーパス”, 43.0682, 141.3595),
(“苗穂アンダーパス”, 43.0684, 141.3725), (“篠路アンダーパス”, 43.1361, 141.3533),
(“百合が原アンダーパス”, 43.1205, 141.3644), (“新川アンダーパス”, 43.0988, 141.3235),
(“菊水アンダーパス”, 43.0618, 141.3776), (“もみじ台通アンダーパス”, 43.0245, 141.4831),
(“上野幌アンダーパス”, 43.0185, 141.4842), (“創成トンネル”, 43.0601, 141.3551),
(“エルムトンネル”, 43.0768, 141.3392), (“平岡跨道橋”, 43.0165, 141.4589),
(“大谷地跨道橋”, 43.0242, 141.4591)
]
loc_names = [loc[0] for loc in LOCATIONS]
==========================================
2. 汎用抽出・計算ロジック
==========================================
def get_latest_init_time(cache_dir):
“””最新のMSM初期時刻を取得”””
files = glob.glob(os.path.join(cache_dir, “MSM_*.npz”))
if not files:
return None
times = sorted(list(set(re.search(r’MSM_(\d{14})‘, f).group(1) for f in files if re.search(r’MSM(\d{14})_’, f))), reverse=True)
return times[0] if times else None
def fuzzy_key_search(data_obj, keywords, exclude=None):
“””NPZファイルから対象変数のキーを柔軟に検索”””
if exclude is None: exclude = [‘lon’, ‘lat’, ‘time’]
for k in data_obj.files:
kl = k.lower()
if any(ex in kl for ex in exclude): continue
if any(kw in kl for kw in keywords): return k
return None
def get_coords_for_data(data_obj, data_array):
“””変数に対応する緯度経度配列を取得”””
ny = data_array.shape[0]; nx = data_array.shape[1] if data_array.ndim > 1 else len(data_array)
for k in data_obj.files:
if ‘lon’ in k.lower():
if (data_obj[k].ndim == 1 and data_obj[k].shape[0] == nx) or (data_obj[k].ndim == 2 and data_obj[k].shape == data_array.shape):
lat_k = k.replace(‘lon’, ‘lat’).replace(‘Lon’, ‘Lat’)
if lat_k in data_obj.files: return data_obj[k], data_obj[lat_k]
return None, None
def get_nearest_index(lon_arr, lat_arr, target_lon, target_lat):
“””指定座標に最も近いグリッドのインデックスを取得”””
if lon_arr.ndim == 1:
lon_2d, lat_2d = np.meshgrid(lon_arr, lat_arr)
else:
lon_2d, lat_2d = lon_arr, lat_arr
dist = (lon_2d – target_lon)2 + (lat_2d – target_lat)2
return np.unravel_index(np.argmin(dist), dist.shape)
def calc_wind_direction(u, v):
“””U,V成分から16方位文字列を計算 (気象学的風向)”””
if u == 0 and v == 0: return “静穏”
deg = (270 – np.degrees(np.arctan2(v, u))) % 360
dirs = [‘北’, ‘北北東’, ‘北東’, ‘東北東’, ‘東’, ‘東南東’, ‘南東’, ‘南南東’,
‘南’, ‘南南西’, ‘南西’, ‘西南西’, ‘西’, ‘西北西’, ‘北西’, ‘北北西’]
idx = int((deg + 11.25) / 22.5) % 16
return dirs[idx]
==========================================
3. メイン処理
==========================================
def main():
if not os.path.exists(CACHE_DIR):
print(f”エラー: ディレクトリ ‘{CACHE_DIR}’ が見つかりません。”)
return
init_time_str = get_latest_init_time(CACHE_DIR)
if not init_time_str:
print("エラー: MSMのNPZデータが見つかりません。")
return
init_dt = datetime.strptime(init_time_str, "%Y%m%d%H%M%S") + timedelta(hours=9)
print(f"対象データ初期時刻: {init_dt.strftime('%Y/%m/%d %H:00')} (JST)")
# 存在するFTファイルのリストを取得
files = glob.glob(os.path.join(CACHE_DIR, f"MSM_{init_time_str}_FT*.npz"))
fts = sorted([int(re.search(r'_FT(\d+)\.npz', f).group(1)) for f in files])
if not fts:
print("エラー: FTファイルが見つかりません。")
return
# 時間ヘッダー作成 (例: "1/20 15:00")
time_headers = [(init_dt + timedelta(hours=ft)).strftime("%m/%d %H:00") for ft in fts]
# 格納用辞書の初期化 (Index: 地点名, Columns: 予測時刻)
df_temp = pd.DataFrame(index=loc_names, columns=time_headers)
df_wdir = pd.DataFrame(index=loc_names, columns=time_headers)
df_wspd = pd.DataFrame(index=loc_names, columns=time_headers)
df_precip = pd.DataFrame(index=loc_names, columns=time_headers)
df_snow = pd.DataFrame(index=loc_names, columns=time_headers)
df_cloud = pd.DataFrame(index=loc_names, columns=time_headers)
print("データ抽出を開始します...")
# 各FTごとにデータを処理
for ft, t_head in zip(fts, time_headers):
file_path = os.path.join(CACHE_DIR, f"MSM_{init_time_str}_FT{ft:02d}.npz")
try:
data = np.load(file_path)
# 各要素のキーを特定
k_t2m = fuzzy_key_search(data, ['t2m', 'tmp2m', 'temp2m'])
k_u10 = fuzzy_key_search(data, ['u10', '10u', 'u_10m'])
k_v10 = fuzzy_key_search(data, ['v10', '10v', 'v_10m'])
k_pr = fuzzy_key_search(data, ['apcp', 'pr', 'precip', 'rain'])
k_sn = fuzzy_key_search(data, ['snow', 'snwe', 'csnow', 'pos']) # 降雪量のキー
k_cc = fuzzy_key_search(data, ['tcc', 'tcdc', 'lcdc', 'cloud'])
# 座標データの取得 (代表として気温の座標を使用)
if k_t2m:
lon_arr, lat_arr = get_coords_for_data(data, data[k_t2m])
else:
continue # 座標がない場合はスキップ
# 各地点ごとの値を抽出
for loc_name, lat, lon in LOCATIONS:
sy, sx = get_nearest_index(lon_arr, lat_arr, lon, lat)
# 1. 気温 (℃)
if k_t2m:
t_val = float(data[k_t2m][sy, sx])
if t_val > 100: t_val -= 273.15 # Kelvin -> Celsius
df_temp.at[loc_name, t_head] = round(t_val, 1)
# 2 & 3. 風向・風速 (16方位, m/s)
if k_u10 and k_v10:
u_val = float(data[k_u10][sy, sx])
v_val = float(data[k_v10][sy, sx])
wspd = np.hypot(u_val, v_val)
wdir = calc_wind_direction(u_val, v_val)
df_wspd.at[loc_name, t_head] = round(wspd, 1)
df_wdir.at[loc_name, t_head] = wdir
# 4. 降水量 (mm)
if k_pr:
pr_val = float(data[k_pr][sy, sx])
df_precip.at[loc_name, t_head] = round(max(pr_val, 0.0), 1)
# 5. 降雪量 (エンジン側の仕様依存。データがない場合は便宜上NaN)
if k_sn:
sn_val = float(data[k_sn][sy, sx])
df_snow.at[loc_name, t_head] = round(max(sn_val, 0.0), 1)
else:
df_snow.at[loc_name, t_head] = np.nan
# 6. 雲量 (%)
if k_cc:
cc_val = float(data[k_cc][sy, sx])
if cc_val <= 1.0 and cc_val > 0: cc_val *= 100 # 0-1スケールをパーセントに変換
df_cloud.at[loc_name, t_head] = round(cc_val, 0)
data.close()
except Exception as e:
print(f"FT={ft:02d} の読み込みでエラーが発生しました: {e}")
# ==========================================
# 4. Excelへの書き出し (C3基準レイアウト)
# ==========================================
print(f"Excelファイル '{OUTPUT_EXCEL}' に出力中...")
with pd.ExcelWriter(OUTPUT_EXCEL, engine='openpyxl') as writer:
# startrow=1, startcol=1 に設定することで、
# 行・列のインデックスが「0始まり」で計算され、データ開始位置が正確に「C3」セルになります。
df_temp.to_excel(writer, sheet_name="気温", startrow=1, startcol=1, index_label="地点名")
df_wdir.to_excel(writer, sheet_name="風向", startrow=1, startcol=1, index_label="地点名")
df_wspd.to_excel(writer, sheet_name="風速", startrow=1, startcol=1, index_label="地点名")
df_precip.to_excel(writer, sheet_name="降水量", startrow=1, startcol=1, index_label="地点名")
df_snow.to_excel(writer, sheet_name="降雪量", startrow=1, startcol=1, index_label="地点名")
df_cloud.to_excel(writer, sheet_name="雲量", startrow=1, startcol=1, index_label="地点名")
print("✅ すべての処理が完了しました！")if name == “main“:
main()
コメント