# -*- coding: utf-8 -*-
"""
補助エンジン idx_engine.py (札幌向け 特殊指数エンジン)
役割
1. 本体エンジンが出力した GPV の npz (MSM / ANAL) を読み、特殊指数を計算して
IDX_{モデル}_{初期時刻}_FT{ft}.npz として出力する(ビューワは IDX_* だけ読めばよい)
2. 気象庁 MGDSST(日別海面水温, 0.25度, 数値データ)を自動取得して SST_latest.npz に保存する
3. 札幌周辺の領域統計を idx_stats.csv に追記する(閾値設計・検証用)
4. IDX_META.json(要素の表示定義)と idx_thresholds.json(危険度の目安・総合スコアの設定)を出力する
使い方
python idx_engine.py # GUIが開く(NPZフォルダ・出力先・範囲をGUIで指定)
python idx_engine.py --autostart # GUIを開いて、前回の設定ですぐ開始
python idx_engine.py --cli --npz-dir D:\\npz --out-dir D:\\idx # GUIなしで実行
python idx_engine.py --cli --npz-dir ... --out-dir ... --all --once # 全データを1回処理して終了
オプション: --hours 12(既定) / --all / --force(最初の1回だけ作り直し) / --sst-file / --no-sst-download / --sst-insecure
注意
* 閾値(idx_thresholds.json)の多くは仮置きです。札幌の過去事例で検証して調整してください。
* 各指標は「気象庁GPVの降水予報とは別系統の環境場診断」であり、降水量そのものの予測ではありません。
"""
import os, re, sys, json, gzip, time, csv, argparse, logging, threading, ssl, urllib.request, urllib.parse, urllib.error
from datetime import datetime, timedelta
import numpy as np
from scipy.ndimage import map_coordinates, gaussian_filter, distance_transform_edt
from scipy.interpolate import RegularGridInterpolator
# ---------------------------------------------------------------- 定数
G = 9.80665; RD = 287.04; CP = 1004.0; LV = 2.5e6; EPS = 0.622; KAPPA = 0.2854
SOURCES_ALL = ("MSM", "ANAL", "GSM_JP")
FILE_RE = re.compile(r'^(MSM|ANAL|GSM_JP)_(\d{14})_FT(\d{2,3})\.npz$')
UPPER_RE = re.compile(r'^(t|u|v|r|w|gh|tddep|ep)(\d{3,4})$|^vort500$')
SAPPORO = (141.328, 43.060)
AREAS = {
"札幌近郊": (141.05, 141.65, 42.85, 43.27),
"石狩湾周辺": (140.6, 141.5, 43.0, 43.8),
"北海道": (139.0, 146.0, 41.5, 45.6),
}
SST_BASES = ["https://www.data.jma.go.jp/goos/data/pub/JMA-product/mgd_sst_glb_D",
"https://www.data.jma.go.jp/gmd/goos/data/pub/JMA-product/mgd_sst_glb_D",
"http://www.data.jma.go.jp/goos/data/pub/JMA-product/mgd_sst_glb_D",
"http://www.data.jma.go.jp/gmd/goos/data/pub/JMA-product/mgd_sst_glb_D"]
SST_CROP = (120.0, 160.0, 20.0, 55.0) # lon_min, lon_max, lat_min, lat_max
BASIS_TEXT = {"A": "一般に使われる目安", "B": "経験則(札幌では未検証)", "C": "仮置き(未検証)", "U": "ご指定の値"}
def _ind(op, caution, warning, basis, note="", **extra):
d = {"op": op, "caution": caution, "warning": warning, "basis": basis, "note": note}
d.update(extra)
return d
_SSI_COLORS = {"caution": "red", "warning": "orange"} # ご指定: 0以下=赤、-3以下=橙(他の指標は 注意=橙/警戒=赤)
_WARM = [6, 7, 8, 9] # 暖候期だけ危険度を評価する指標用
DEFAULT_CFG = {
"version": 3,
"_note": ("危険度の目安と総合スコアの設定。数値はすべて調整可能です。変更後、地図とスコアに反映するには "
"idx_engine.py を --force 付きで再実行してください(図の危険度ラインと注釈は再起動なしで更新されます)。"
"basis は根拠の強さ: A=一般に使われる目安 / B=経験則(札幌未検証) / C=仮置き(未検証) / U=ご指定の値。"
"op が ge は『以上』、le は『以下』で危険。"
"colors で線と表の色を変えられます(caution/warning に red か orange)。"
"active_months を付けた指標は、その月だけ危険度を評価します(月以外は参考表示)。"
"実測データから閾値を提案するには --suggest-thresholds を使います。"),
"indicators": {
"MFC": _ind("ge", 2.0, 4.0, "C", "下層の水蒸気収束による降水の供給量(mm/h)"),
"DIFF_SUP": _ind("ge", 2.0, 4.0, "C", "供給量が予報降水を上回る分。大きいほどモデルが降水を出し切れていない可能性"),
"IPW": _ind("ge", 35.0, 45.0, "C", "北日本の夏としては多めの水蒸気量という感覚的な目安"),
"IVT": _ind("ge", 250.0, 500.0, "B", "大気の川の基準(250 kg/m/s)が目安。本指標は975-500hPa積分のためやや小さめに出る"),
"QFLX850": _ind("ge", 150.0, 250.0, "C", "850hPaの水蒸気流入の強さ"),
"QFLX925": _ind("ge", 150.0, 250.0, "C", "925hPaの水蒸気流入の強さ"),
"DIV925": _ind("le", -3.0, -6.0, "C", "下層収束(負)が強いほど持ち上げが強い"),
"K": _ind("ge", 30.0, 35.0, "A", "雷雨の目安としてよく使われる値(米国由来)。札幌では未検証"),
"SSI": _ind("le", 0.0, -3.0, "A", "850→500hPa。0以下で不安定、-3以下でかなり不安定という古典的な目安", colors=_SSI_COLORS),
"SSI_850_700": _ind("le", 0.0, -3.0, "U", "850→700hPa。SSIと同じ値をご指定により設定", colors=_SSI_COLORS),
"SSI_925_700": _ind("le", 0.0, -3.0, "U", "925→700hPa。SSIと同じ値をご指定により設定", colors=_SSI_COLORS),
"TT": _ind("ge", 44.0, 50.0, "A", "44以上で雷雨の可能性、50以上で強い、という古典的な目安"),
"DTHE": _ind("ge", 10.0, 20.0, "C", "θe(850)-θe(500)。大きいほど対流不安定"),
"WCD": _ind("ge", 3000.0, 4000.0, "C", "暖かい雲層が厚いほど降水効率が高い傾向"),
"SHR925_700": _ind("ge", 12.0, 18.0, "C", "暖候期のみ評価(冬は偏西風で常に大きいため対象外)。見直し: 季節別+実測分布で調整", active_months=_WARM),
"SHR925_500": _ind("ge", 18.0, 25.0, "C", "暖候期のみ評価(冬は偏西風で常に大きいため対象外)。見直し: 季節別+実測分布で調整", active_months=_WARM),
"DT850": _ind("ge", 15.0, 20.0, "B", "本州側日本海で筋状雲が発達する経験則。石狩湾は海面水温が低く要校正"),
"DT500": _ind("ge", 36.0, 40.0, "B", "同上。大きいほど不安定"),
"FETCH": _ind("ge", 150.0, 300.0, "C", "風上の連続した海上距離(km)"),
"HEATPATH": _ind("ge", 20.0, 40.0, "C", "風上の海上での気団変質の大きさ"),
"DT850_UP": _ind("ge", 15.0, 20.0, "B", "風上の海上経路での平均。DT850と同じ経験則"),
"DT500_UP": _ind("ge", 36.0, 40.0, "B", "風上の海上経路での平均。DT500と同じ経験則"),
"SENS": _ind("ge", 60.0, 120.0, "C", "風速×(SST-気温)。海面から大気への熱供給の目安"),
"LATENT": _ind("ge", 60.0, 120.0, "C", "風速×(海面と大気の比湿差)。水蒸気供給の目安"),
"DIV10": _ind("le", -4.0, -8.0, "C", "地上収束帯(石狩湾周辺など)の検出用"),
"MOISTTHICK": _ind("ge", 2000.0, 3000.0, "C", "雲層の厚さの代用"),
"DGZTHICK": _ind("ge", 600.0, 1000.0, "C", "-12〜-18℃かつ湿潤な層の厚さ。降雪効率の目安"),
},
"scores": {
"summer": {
"label": "夏型(大雨)",
"items": [{"key": k, "weight": 1} for k in ("MFC", "QFLX850", "IPW", "K", "SSI", "DIV925", "WCD")],
"judge": {"caution": 5, "warning": 9},
"gate": None,
},
"winter": {
"label": "冬型(大雪)",
"items": [{"key": k, "weight": 1} for k in ("DT850_UP", "DT500_UP", "FETCH", "DIV10", "DGZTHICK", "MOISTTHICK")],
"judge": {"caution": 4, "warning": 7},
"gate": {"key": "Z0", "op": "le", "value": 1500.0, "label": "0℃高度が1500m以下(寒気が入っている)"},
},
},
}
# ---------------------------------------------------------------- 表示定義(ビューワが読む)
CATEGORIES = ["総合", "量", "収束", "安定度", "雷", "気温・雲層", "風・雲の動き", "海上・冬型"]
def _el(key, label, unit, cat, grid, levels, cmap, extend, desc, short=None, overlay=True, ovl=None, ovnorm=False):
if not overlay: ovl = None
elif ovl is None:
ovl = dict(u='OV_U925', v='OV_V925', grid='p') if grid == 'p' else dict(u='OV_U10', v='OV_V10', grid='s')
else:
ovl = dict(ovl)
if ovl is not None: ovl['norm'] = bool(ovnorm)
return dict(key=key, label=label, short=short or label, unit=unit, season=cat, cat=cat, grid=grid, levels=levels,
cmap=cmap, extend=extend, desc=desc, overlay=overlay, ovl=ovl)
_OV850 = dict(u='OV_U850', v='OV_V850', grid='p')
META_ELEMENTS = [
# --- 総合
_el("SCORE_SUMMER", "夏の危険度スコア(大雨)", "点", "総合", "p", [0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5], "YlOrRd", "max",
"条件成立数に応じた点数。条件は idx_thresholds.json で調整。", short="夏スコア"),
_el("SCORE_WINTER", "冬の危険度スコア(大雪)", "点", "総合", "p", [0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5], "PuBu", "max",
"条件成立数に応じた点数。条件は idx_thresholds.json で調整。", short="冬スコア"),
# --- 量
_el("IPW", "可降水量(概算)", "mm", "量", "p", list(range(10, 66, 5)), "GnBu", "max",
"975〜300hPaの比湿を積分した概算値(地表〜975hPaは含まない)。", short="可降水量"),
_el("IVT", "IVT(水蒸気輸送量)", "kg/m/s", "量", "p", [100, 150, 200, 250, 300, 400, 500, 700], "PuRd", "max",
"975〜500hPaの q・V を積分した水蒸気輸送量の大きさ。", short="IVT"),
_el("QFLX850", "850hPa水蒸気フラックス(q×V)", "g/kg·m/s", "量", "p", [50, 100, 150, 200, 250, 300, 400], "YlOrRd", "max",
"比湿[g/kg]×風速[m/s]。下層への水蒸気流入の強さ。", short="水蒸気流入850"),
_el("QFLX925", "925hPa水蒸気フラックス(q×V)", "g/kg·m/s", "量", "p", [50, 100, 150, 200, 250, 300, 400], "YlOrRd", "max",
"比湿[g/kg]×風速[m/s]。", short="水蒸気流入925"),
_el("DIFF_SUP", "水蒸気供給量 − MSM予報降水量", "mm/h", "量", "p", [-8, -4, -2, 0, 2, 4, 6, 8, 12], "RdBu_r", "both",
"プラスが大きい所は、環境場の供給に対してモデルの降水が小さい(出し切れていない可能性)。要検証の目安。", short="供給−予報降水"),
# --- 収束
_el("MFC", "水蒸気収束供給量(975-700hPa)", "mm/h", "収束", "p", [0.5, 1, 2, 3, 4, 6, 8, 12, 16], "YlGnBu", "max",
"下層水蒸気の収束を鉛直積分して mm/h に換算した値。降水の『供給側』の目安で、降水量そのものではない。", short="水蒸気収束"),
_el("DIV925", "925hPa風の発散(負=収束)", "1e-5/s", "収束", "p", [-12, -8, -6, -4, -2, 0, 2, 4, 6], "RdBu", "both",
"負の領域が下層収束。降水の持ち上げ要因の目安。", short="発散925"),
_el("DIV10", "地上風の発散(負=収束)", "1e-5/s", "収束", "s", [-20, -12, -8, -4, -2, 0, 2, 4, 8], "RdBu", "both",
"石狩湾周辺の収束帯の検出用。格子が細かいので平滑化済み。", short="発散10m"),
# --- 安定度
_el("SSI", "ショワルター安定指数 SSI(850→500hPa)", "℃", "安定度", "p", [-8, -6, -4, -3, -2, 0, 2, 4, 6], "RdYlBu", "both",
"850hPa気塊を500hPaまで断熱上昇させた時の環境との差。負が不安定(寒色=安定)。", short="SSI 850-500"),
_el("SSI_850_700", "SSI(850→700hPa)", "℃", "安定度", "p", [-8, -6, -4, -3, -2, 0, 2, 4, 6], "RdYlBu", "both",
"850hPa気塊を700hPaまで持ち上げた時の環境との差。下層の不安定の目安。負が不安定。", short="SSI 850-700"),
_el("SSI_925_700", "SSI(925→700hPa)", "℃", "安定度", "p", [-8, -6, -4, -3, -2, 0, 2, 4, 6], "RdYlBu", "both",
"925hPa気塊を700hPaまで持ち上げた時の環境との差。最下層の不安定の目安。負が不安定。", short="SSI 925-700"),
_el("DTHE", "θe(850)−θe(500)", "K", "安定度", "p", [0, 5, 10, 15, 20, 25, 30], "YlOrRd", "max", "大きいほど対流不安定。", short="θe差"),
# --- 雷
_el("K", "Kインデックス", "℃", "雷", "p", [15, 20, 25, 30, 33, 36, 40], "YlOrRd", "max", "(T850−T500)+Td850−(T700−Td700)", short="K指数"),
_el("TT", "Total Totals", "℃", "雷", "p", [38, 42, 44, 46, 48, 50, 52, 56], "YlOrRd", "max", "T850+Td850−2×T500", short="TT"),
# --- 気温・雲層
_el("Z0", "0℃高度", "m", "気温・雲層", "p", [0, 1000, 2000, 3000, 3500, 4000, 4500, 5000], "coolwarm", "max", "気圧面データからの内挿。", short="0℃高度"),
_el("Z3", "3℃高度(気温が3℃になる高さ)", "m", "気温・雲層", "p", [0, 250, 500, 1000, 1500, 2000, 3000, 4000], "coolwarm", "max",
"気温が3℃になる高さ。最下層(975hPa)がすでに3℃未満なら0m。低いほど雪に近い環境(気圧面データからの内挿)。", short="3℃高度"),
_el("WCD", "暖かい雲層の厚さ(0℃高度−LCL高度)", "m", "気温・雲層", "p", [0, 1000, 2000, 3000, 3500, 4000, 4500, 5000], "YlGnBu", "max",
"厚いほど降水効率が高い傾向(概算)。", short="暖雲厚"),
_el("MOISTTHICK", "湿潤層厚(RH≥80%)", "m", "気温・雲層", "p", [0, 500, 1000, 1500, 2000, 2500, 3000, 4000], "GnBu", "max",
"雲層の厚さの代用(気圧面の層から概算)。", short="湿潤層厚"),
_el("CLDTHICK", "下層から連続する雲層厚(RH≥80%)", "m", "気温・雲層", "p", [0, 500, 1000, 1500, 2000, 3000, 4000, 6000], "GnBu", "max",
"最下層(975hPa)から上へ、RH≥80%が途切れず続く層の厚さ。最下層が乾いていれば0。", short="下層連続雲厚"),
_el("TOP80", "湿潤層の最高高度(雲頂の代用)", "m", "気温・雲層", "p", [0, 1000, 2000, 3000, 4000, 5000, 6000], "GnBu", "max",
"RH≥80%を満たす最高の気圧面の高度。", short="雲頂高度"),
_el("COLDCLD", "0℃以上の高さの雲層厚(0℃高度より上のRH≥80%)", "m", "気温・雲層", "p", [0, 500, 1000, 2000, 3000, 4000, 6000], "PuBu", "max",
"気温が0℃未満の側(0℃高度より上)にある湿潤層(RH≥80%)の厚さ。氷晶過程が働く雲の厚さの目安。", short="0℃上の雲厚"),
_el("DGZTHICK", "樹枝状結晶成長層の湿潤厚(−12〜−18℃)", "m", "気温・雲層", "p", [0, 200, 400, 600, 800, 1000, 1500, 2000], "PuBu", "max",
"−12〜−18℃かつRH≥80%の層の厚さ。降雪効率の目安。", short="DGZ厚"),
# --- 風・雲の動き
_el("WS925", "925hPa風速(下層の流入)", "m/s", "風・雲の動き", "p", [5, 8, 10, 15, 20, 25, 30], "PuBu", "max", "矢印=風向。水蒸気の流入の強さの目安。", short="風速925"),
_el("WS850", "850hPa風速(下層ジェット)", "m/s", "風・雲の動き", "p", [5, 8, 10, 15, 20, 25, 30], "PuBu", "max", "矢印=風向。下層ジェットの目安。",
short="風速850", ovl=_OV850),
_el("CLDWIND", "雲層平均風(850〜300hPa)", "m/s", "風・雲の動き", "p", [2, 5, 8, 10, 15, 20, 30], "YlGnBu", "max",
"850・700・500・300hPaの風の圧力加重平均。発達した雲が流される向きと速さの目安。矢印=向き。", short="雲層平均風",
ovl=dict(u='OV_UCLD', v='OV_VCLD', grid='p'), ovnorm=True),
_el("MCSMOVE", "積乱雲群(MCS)の移動速度(Corfidi法)", "m/s", "風・雲の動き", "p", [1, 3, 5, 8, 10, 15, 20], "YlOrBr", "max",
"2×雲層平均風 − 850hPa風。小さいほど同じ場所に雨が続きやすい(水蒸気量・不安定と合わせて見る)。矢印=移動の向き。", short="MCS移動速",
ovl=dict(u='OV_UMCS', v='OV_VMCS', grid='p'), ovnorm=True),
_el("PROPAG", "新規セル(子雲)の伝播ベクトル(Corfidi法)", "m/s", "風・雲の動き", "p", [1, 3, 5, 8, 10, 15, 20], "OrRd", "max",
"雲層平均風 − 850hPa風。新しい積乱雲が既存の雲に対してできやすい向きと速さの目安。矢印=新しいセルができやすい側。", short="新規セル伝播",
ovl=dict(u='OV_UPRP', v='OV_VPRP', grid='p'), ovnorm=True),
_el("SHR925_700", "鉛直シア(925-700hPa)", "m/s", "風・雲の動き", "p", [2, 4, 6, 8, 10, 15, 20], "PuBu", "max", "ベクトル差の大きさ。", short="シア925-700"),
_el("SHR925_500", "鉛直シア(925-500hPa)", "m/s", "風・雲の動き", "p", [5, 10, 15, 20, 25, 30, 40], "PuBu", "max",
"下層風と中層風の差(組織化・移動の目安)。", short="シア925-500"),
# --- 海上・冬型
_el("SST", "海面水温(MGDSST)", "℃", "海上・冬型", "p", [0, 4, 8, 12, 16, 20, 24, 28], "turbo", "max",
"気象庁MGDSSTの日別海面水温(解析日はSST状態表示を参照)。", short="SST"),
_el("DT850", "SST−T850(海上のみ)", "℃", "海上・冬型", "p", [10, 12, 15, 18, 20, 22, 25], "YlOrRd", "max",
"海面水温と850hPa気温の差。15℃以上で筋状雲が発達する目安(本州側日本海の経験則。札幌は要校正)。", short="SST−T850"),
_el("DT500", "SST−T500(海上のみ)", "℃", "海上・冬型", "p", [30, 33, 36, 38, 40, 42, 45], "YlOrRd", "max",
"海面水温と500hPa気温の差。大きいほど不安定(経験則。要校正)。", short="SST−T500"),
_el("DT850_UP", "風上平均 SST−T850", "℃", "海上・冬型", "p", [8, 10, 12, 15, 18, 20, 25], "YlOrRd", "max",
"風上の海上経路での平均値。陸上(札幌など)でも評価できる。", short="風上SST−T850"),
_el("DT500_UP", "風上平均 SST−T500", "℃", "海上・冬型", "p", [30, 33, 36, 38, 40, 42, 45], "YlOrRd", "max", "風上の海上経路での平均値。", short="風上SST−T500"),
_el("FETCH", "海上吹走距離(925hPa後方流跡)", "km", "海上・冬型", "p", [0, 50, 100, 150, 200, 300, 400, 600], "YlGnBu", "max",
"各地点から925hPa風の風上へ遡り、連続して海上だった距離。陸上の点も風上の海を評価する。", short="吹走距離"),
_el("HEATPATH", "吹走熱量指標(風上のSST−T925積算)", "℃·100km", "海上・冬型", "p", [0, 5, 10, 20, 30, 40, 60, 80], "YlOrRd", "max",
"風上の海上経路でΣ max(SST−T925,0)×距離/100km。気団変質の大きさの目安。", short="吹走熱量"),
_el("SENS", "海面顕熱指標 |V10|×(SST−T2m)", "℃·m/s", "海上・冬型", "s", [0, 20, 40, 60, 80, 100, 150, 200], "YlOrRd", "max",
"風を含めた海面顕熱供給の目安(海上のみ)。", short="海面顕熱"),
_el("LATENT", "海面潜熱指標 |V10|×(qs(SST)−q2m)", "g/kg·m/s", "海上・冬型", "s", [0, 20, 40, 60, 80, 100, 150, 200], "YlGnBu", "max",
"風を含めた海面水蒸気供給の目安(海上のみ)。", short="海面潜熱"),
]
# ---------------------------------------------------------------- 物理ユーティリティ
def es_hpa(tc):
return 6.112 * np.exp(17.67 * tc / (tc + 243.5))
def q_from_t_rh(tc, rh, p_hpa):
e = es_hpa(tc) * np.clip(rh, 0.1, 100.0) / 100.0
return EPS * e / (p_hpa - (1.0 - EPS) * e)
def td_from_t_rh(tc, rh):
e = es_hpa(tc) * np.clip(rh, 1.0, 100.0) / 100.0
ln = np.log(e / 6.112)
return 243.5 * ln / (17.67 - ln)
def theta_e(tc, rh, p_hpa):
tk = tc + 273.15
q = q_from_t_rh(tc, rh, p_hpa)
w = q / (1.0 - q)
return tk * (1000.0 / p_hpa) ** KAPPA * np.exp(LV * w / (CP * tk))
def wet_bulb(tc, rh): # Stull (2011)
rh = np.clip(rh, 5.0, 99.0)
return (tc * np.arctan(0.151977 * np.sqrt(rh + 8.313659)) + np.arctan(tc + rh)
- np.arctan(rh - 1.676331) + 0.00391838 * rh ** 1.5 * np.arctan(0.023101 * rh) - 4.686035)
def smooth(a, sigma):
m = np.isfinite(a)
if not m.any():
return a
f = gaussian_filter(np.where(m, a, 0.0), sigma)
w = gaussian_filter(m.astype(float), sigma)
return np.where(m, f / np.maximum(w, 1e-6), np.nan)
def div_field(u, v, lon, lat):
"""球面上の発散 (1/s)。lon, lat は昇順の1次元。"""
R = 6371000.0
dlon = np.deg2rad(np.gradient(lon)); dlat = np.deg2rad(np.gradient(lat))
cl = np.cos(np.deg2rad(lat))[:, None]
t1 = np.gradient(u, axis=1) / (R * cl * dlon[None, :])
t2 = np.gradient(v * cl, axis=0) / (R * cl * dlat[:, None])
return t1 + t2
def regrid(arr, src_lon, src_lat, dst_lon, dst_lat):
f = RegularGridInterpolator((src_lat, src_lon), arr, bounds_error=False, fill_value=np.nan)
LAT, LON = np.meshgrid(dst_lat, dst_lon, indexing='ij')
return f(np.stack([LAT.ravel(), LON.ravel()], -1)).reshape(LAT.shape)
def sst_to_grid(sst, slon, slat, dlon, dlat):
"""SST(陸=NaN)を目的格子へ。陸の判定は有効値マスクの内挿(>=0.5)で行う。"""
valid = np.isfinite(sst)
if not valid.any():
return np.full((len(dlat), len(dlon)), np.nan)
idx = distance_transform_edt(~valid, return_distances=False, return_indices=True)
filled = sst[tuple(idx)]
out = regrid(filled, slon, slat, dlon, dlat)
m = regrid(valid.astype(float), slon, slat, dlon, dlat) >= 0.5
return np.where(m, out, np.nan)
def level_heights(A, lv):
"""各気圧面の高度(m)を返す。ghがある面はそれを使い、無い面は下の面から静力学式で補う。"""
Z = {}
for i, L in enumerate(lv):
g_ = A.get(f'gh{L}')
if g_ is not None and np.isfinite(g_).any():
Z[L] = g_; continue
if i == 0:
z_std = 44330.0 * (1.0 - (L / 1013.25) ** 0.1903)
Z[L] = np.full(A[f't{L}'].shape, z_std)
else:
Lp = lv[i - 1]
tm = 0.5 * (A[f't{L}'] + A[f't{Lp}']) + 273.15
Z[L] = Z[Lp] + RD * tm / G * np.log(Lp / L)
return Z
def lift_parcel(t0c, td0c, p0, p_end=500.0, n=30):
"""p0(hPa)の気塊を p_end まで持ち上げた時の気温(K)を返す。
乾燥断熱でLCLまで(解析解)→ LCL から p_end まで湿潤断熱(2次精度のルンゲクッタ、n分割)。
旧方式(10hPa刻みの前進オイラー+LCLのステップ誤差)は気塊温度が最大約1K高めに出ていたため置き換えた。"""
t0 = t0c + 273.15; td0 = np.minimum(td0c, t0c) + 273.15
tlcl = 1.0 / (1.0 / (td0 - 56.0) + np.log(t0 / td0) / 800.0) + 56.0
plcl = p0 * (tlcl / t0) ** (1.0 / KAPPA)
t_dry_end = t0 * (p_end / p0) ** KAPPA
def f(T, p):
es = es_hpa(T - 273.15); rs = EPS * es / np.maximum(p - es, 1.0)
return (RD * T + LV * rs) / (p * (CP + (LV ** 2 * rs * EPS) / (RD * T ** 2)))
p = np.maximum(plcl, p_end); h = (p - p_end) / n # 点ごとに刻み幅を変える(LCLを厳密に扱うため)
T = np.where(plcl > p_end, tlcl, t_dry_end)
for _ in range(n):
Tm = T - f(T, p) * h / 2.0
T = T - f(Tm, p - h / 2.0) * h
p = p - h
return np.where(plcl > p_end, T, t_dry_end)
# ---------------------------------------------------------------- npz 読み込み
def _pack(lon, lat, a):
lon = lon.astype(float); lat = lat.astype(float)
if lat[0] > lat[-1]:
lat = lat[::-1]; a = {k: v[::-1, :] for k, v in a.items()}
if lon[0] > lon[-1]:
lon = lon[::-1]; a = {k: v[:, ::-1] for k, v in a.items()}
clean = {}
for k, v in a.items():
v = np.asarray(v, dtype=np.float64)
v = np.where(np.abs(v) > 1e15, np.nan, v)
if k.startswith('gh') and np.isfinite(v).any() and np.nanmax(v) > 30000.0:
v = v / G # ジオポテンシャル(m2/s2)で格納されていた場合はmへ換算
clean[k] = v
return {'lon': lon, 'lat': lat, 'a': clean}
def load_fields(path):
z = np.load(path, allow_pickle=False)
arrs = {k: np.asarray(z[k]) for k in z.files}
def pick_grid(prefs, ref):
shp = arrs[ref].shape if ref in arrs else None
for lk, tk in prefs:
if lk in arrs and tk in arrs and arrs[lk].ndim == 1 and arrs[tk].ndim == 1:
if shp is None or (len(arrs[tk]), len(arrs[lk])) == shp:
return arrs[lk], arrs[tk]
return None
up = [k for k in arrs if UPPER_RE.match(k) and arrs[k].ndim == 2]
sf = [k for k in arrs if arrs[k].ndim == 2 and not UPPER_RE.match(k) and not k.startswith(('lon', 'lat'))]
res = {'P': None, 'S': None}
if up:
g = pick_grid([('lon_pall', 'lat_pall'), ('lon_p', 'lat_p'), ('lon', 'lat')], up[0])
if g is not None:
shp = (len(g[1]), len(g[0]))
res['P'] = _pack(g[0], g[1], {k: arrs[k] for k in up if arrs[k].shape == shp})
if sf:
ref = 't2m' if 't2m' in arrs else sf[0]
g = pick_grid([('lon_surf', 'lat_surf'), ('lon', 'lat')], ref)
if g is not None:
shp = (len(g[1]), len(g[0]))
res['S'] = _pack(g[0], g[1], {k: arrs[k] for k in sf if arrs[k].shape == shp})
return res
# ---------------------------------------------------------------- 指標計算
def _levels(A, names, cands):
return [L for L in cands if all(f'{n}{L}' in A for n in names)]
def upstream_fetch(u, v, sst, t925, t850, t500, lon, lat, ds_km=10.0, n_steps=60):
"""各格子点から風上へ遡り、連続した海上区間の距離・熱量・SST-T平均を求める。"""
ny, nx = u.shape
LON, LAT = np.meshgrid(lon, lat)
dlon = (lon[-1] - lon[0]) / (nx - 1); dlat = (lat[-1] - lat[0]) / (ny - 1)
sea = np.isfinite(sst)
u0 = np.nan_to_num(u); v0 = np.nan_to_num(v); sst0 = np.where(sea, sst, 0.0); seaf = sea.astype(float)
t925_0 = np.nan_to_num(t925); t850_0 = np.nan_to_num(t850); t500_0 = np.nan_to_num(t500)
def samp(A, x, y):
return map_coordinates(A, [(y - lat[0]) / dlat, (x - lon[0]) / dlon], order=1, mode='nearest')
x = LON.copy(); y = LAT.copy()
alive = np.ones((ny, nx), bool); seen = np.zeros((ny, nx), bool)
dist = np.zeros((ny, nx)); heat = np.zeros((ny, nx)); d850 = np.zeros((ny, nx)); d500 = np.zeros((ny, nx)); cnt = np.zeros((ny, nx))
for _ in range(n_steps):
uu = samp(u0, x, y); vv = samp(v0, x, y)
sp = np.maximum(np.hypot(uu, vv), 0.5)
x = x - uu / sp * ds_km / (111.0 * np.cos(np.deg2rad(y)))
y = y - vv / sp * ds_km / 111.0
inside = (x >= lon[0]) & (x <= lon[-1]) & (y >= lat[0]) & (y <= lat[-1])
is_sea = samp(seaf, x, y) > 0.5
alive &= inside
alive &= ~(seen & ~is_sea)
ok = alive & is_sea
seen |= ok
s = samp(sst0, x, y)
dist += ok * ds_km
heat += ok * np.maximum(s - samp(t925_0, x, y), 0.0) * ds_km / 100.0
d850 += ok * (s - samp(t850_0, x, y)); d500 += ok * (s - samp(t500_0, x, y)); cnt += ok
with np.errstate(all='ignore'):
d850u = np.where(cnt > 0, d850 / cnt, np.nan); d500u = np.where(cnt > 0, d500 / cnt, np.nan)
return dist, heat, d850u, d500u
def pts_arr(a, ind):
"""0点(基準未満) / 1点(注意) / 2点(警戒)"""
if ind['op'] == 'ge':
p = np.where(a >= ind['warning'], 2.0, np.where(a >= ind['caution'], 1.0, 0.0))
else:
p = np.where(a <= ind['warning'], 2.0, np.where(a <= ind['caution'], 1.0, 0.0))
return np.where(np.isfinite(a), p, 0.0)
def pts_scalar(v, ind):
if v is None: return None
if ind['op'] == 'ge': return 2 if v >= ind['warning'] else (1 if v >= ind['caution'] else 0)
return 2 if v <= ind['warning'] else (1 if v <= ind['caution'] else 0)
def _fl(x):
x = float(x)
return x if np.isfinite(x) else None
def level_of(total, judge):
if total is None: return None
return '警戒' if total >= judge['warning'] else ('注意' if total >= judge['caution'] else '平常')
def build_table(cfg, out, P, get_p, raws):
"""総合判定表(札幌点と札幌近郊の、各条件の値・閾値・点数)をJSON化できる形で返す。"""
lon, lat = P['lon'], P['lat']
j = int(np.abs(lat - SAPPORO[1]).argmin()); i = int(np.abs(lon - SAPPORO[0]).argmin())
x0, x1, y0, y1 = AREAS['札幌近郊']
mx = (lon >= x0) & (lon <= x1); my = (lat >= y0) & (lat <= y1)
has_box = bool(mx.any() and my.any())
def at_pt(a): return _fl(a[j, i])
def at_box(a, op):
if not has_box: return None
sub = a[np.ix_(my, mx)]
if not np.isfinite(sub).any(): return None
return _fl(np.nanmax(sub) if op == 'ge' else np.nanmin(sub))
elabel = {e['key']: (e['label'], e['unit']) for e in META_ELEMENTS}
tbl = {'version': 1, 'sapporo': list(SAPPORO), 'grid_pt': [_fl(lon[i]), _fl(lat[j])], 'area': '札幌近郊'}
for sk, okey in (('summer', 'SCORE_SUMMER'), ('winter', 'SCORE_WINTER')):
sc = cfg.get('scores', {}).get(sk)
if not sc or sk not in raws: continue
rows = []
for it in sc['items']:
key = it['key']; ind = cfg['indicators'].get(key); lab, unit = elabel.get(key, (key, ''))
arr = get_p(key)
row = {'key': key, 'label': lab, 'unit': unit, 'weight': it.get('weight', 1)}
if ind is None or arr is None:
row['missing'] = True; rows.append(row); continue
vp = at_pt(arr); vb = at_box(arr, ind['op'])
row.update(op=ind['op'], caution=ind['caution'], warning=ind['warning'], basis=ind.get('basis', ''),
v_pt=vp, v_box=vb, p_pt=pts_scalar(vp, ind), p_box=pts_scalar(vb, ind))
rows.append(row)
shown = out[okey][1]; judge = sc['judge']
rp = at_pt(raws[sk]); rb = at_box(raws[sk], 'ge'); sp = at_pt(shown); sb = at_box(shown, 'ge')
gate = None; g = sc.get('gate')
if g:
ga = get_p(g['key'])
if ga is not None:
gop = g.get('op', 'le'); gv = at_pt(ga); gb = at_box(ga, 'ge' if gop == 'ge' else 'le')
okf = lambda v: None if v is None else bool((v <= g['value']) if gop == 'le' else (v >= g['value']))
gate = {'label': g.get('label', g['key']), 'key': g['key'], 'op': gop, 'value': g['value'],
'unit': elabel.get(g['key'], ('', ''))[1], 'v_pt': gv, 'ok_pt': okf(gv), 'v_box': gb, 'ok_box': okf(gb)}
tbl[sk] = {'label': sc.get('label', sk), 'max_points': 2 * sum(float(it.get('weight', 1)) for it in sc['items']),
'judge': judge, 'gate': gate, 'rows': rows,
'pt': {'raw': rp, 'shown': sp, 'level': level_of(sp, judge), 'raw_level': level_of(rp, judge)},
'box': {'raw': rb, 'shown': sb, 'level': level_of(sb, judge), 'raw_level': level_of(rb, judge)}}
return tbl
def compute_indices(F, sst, thr, log):
"""F: load_fields の結果, sst: dict(lon,lat,sst) or None。戻り値: {key: (grid, array)}"""
out = {}
P = F.get('P'); S = F.get('S')
def guard(label, fn):
try:
fn()
except Exception as e: # 1要素の失敗で全体を止めない
log.warning(f" 指標 {label} をスキップ: {e}")
sstP = sstS = None
if sst is not None:
if P is not None: sstP = sst_to_grid(sst['sst'], sst['lon'], sst['lat'], P['lon'], P['lat'])
if S is not None: sstS = sst_to_grid(sst['sst'], sst['lon'], sst['lat'], S['lon'], S['lat'])
with np.errstate(all='ignore'):
if P is not None:
A = P['a']; lon = P['lon']; lat = P['lat']
def _mfc():
lv = _levels(A, 'truv', [975, 950, 925, 850, 700])
if len(lv) < 3: return
divq = []
for L in lv:
q = q_from_t_rh(A[f't{L}'], A[f'r{L}'], L)
divq.append(div_field(q * A[f'u{L}'], q * A[f'v{L}'], lon, lat))
m = np.zeros_like(divq[0])
for i in range(len(lv) - 1):
m += -0.5 * (divq[i] + divq[i + 1]) * (lv[i] - lv[i + 1]) * 100.0 / G
out['MFC'] = ('p', smooth(m * 3600.0, 1.0))
guard('MFC', _mfc)
def _ipw():
lv = _levels(A, 'tr', [975, 950, 925, 850, 700, 600, 500, 300])
if len(lv) < 3: return
qs = [q_from_t_rh(A[f't{L}'], A[f'r{L}'], L) for L in lv]
s = sum(0.5 * (qs[i] + qs[i + 1]) * (lv[i] - lv[i + 1]) * 100.0 / G for i in range(len(lv) - 1))
out['IPW'] = ('p', s)
guard('IPW', _ipw)
def _ivt():
lv = _levels(A, 'truv', [975, 950, 925, 850, 700, 600, 500])
if len(lv) < 3: return
qu = []; qv = []
for L in lv:
q = q_from_t_rh(A[f't{L}'], A[f'r{L}'], L); qu.append(q * A[f'u{L}']); qv.append(q * A[f'v{L}'])
iu = sum(0.5 * (qu[i] + qu[i + 1]) * (lv[i] - lv[i + 1]) * 100.0 / G for i in range(len(lv) - 1))
iv = sum(0.5 * (qv[i] + qv[i + 1]) * (lv[i] - lv[i + 1]) * 100.0 / G for i in range(len(lv) - 1))
out['IVT'] = ('p', np.hypot(iu, iv))
guard('IVT', _ivt)
def _qflx():
for L in (850, 925):
if all(f'{n}{L}' in A for n in 'truv'):
q = q_from_t_rh(A[f't{L}'], A[f'r{L}'], L) * 1000.0
out[f'QFLX{L}'] = ('p', q * np.hypot(A[f'u{L}'], A[f'v{L}']))
guard('QFLX', _qflx)
def _div925():
for L in (925, 950, 975):
if f'u{L}' in A and f'v{L}' in A:
out['DIV925'] = ('p', smooth(div_field(A[f'u{L}'], A[f'v{L}'], lon, lat) * 1e5, 1.0)); return
guard('DIV925', _div925)
def _stab():
if all(k in A for k in ('t850', 'r850', 't500')):
td850 = td_from_t_rh(A['t850'], A['r850'])
tp = lift_parcel(A['t850'], td850, 850.0)
out['SSI'] = ('p', A['t500'] - (tp - 273.15))
out['TT'] = ('p', A['t850'] + td850 - 2.0 * A['t500'])
if all(k in A for k in ('t700', 'r700')):
td700 = td_from_t_rh(A['t700'], A['r700'])
out['K'] = ('p', (A['t850'] - A['t500']) + td850 - (A['t700'] - td700))
out['DTHE'] = ('p', theta_e(A['t850'], A['r850'], 850.0) - theta_e(A['t500'], A['r500'], 500.0)) if 'r500' in A else None
if out.get('DTHE') is None: out.pop('DTHE', None)
guard('安定度', _stab)
def _ssi_extra():
for b, t_, key in ((850, 700, 'SSI_850_700'), (925, 700, 'SSI_925_700')):
if all(k in A for k in (f't{b}', f'r{b}', f't{t_}')):
tdb = td_from_t_rh(A[f't{b}'], A[f'r{b}'])
tp = lift_parcel(A[f't{b}'], tdb, float(b), p_end=float(t_))
out[key] = ('p', A[f't{t_}'] - (tp - 273.15))
guard('SSI(下層)', _ssi_extra)
def _zero():
lv = _levels(A, ('t',), [975, 950, 925, 850, 700, 600, 500])
if len(lv) < 3:
log.warning(f' Z0/WCD: 必要な気圧面が不足 {lv}'); return
Zd = level_heights(A, lv)
T = [A[f't{L}'] for L in lv]; Z = [Zd[L] for L in lv]
z0 = np.full(T[0].shape, np.nan)
for i in range(len(lv) - 1):
cross = np.isnan(z0) & (T[i] >= 0) & (T[i + 1] < 0)
zz = Z[i] + T[i] / np.where(T[i] - T[i + 1] == 0, np.nan, T[i] - T[i + 1]) * (Z[i + 1] - Z[i])
z0 = np.where(cross, zz, z0)
z0 = np.where(np.isnan(z0) & (T[0] < 0), 0.0, z0)
z0 = np.where(np.isnan(z0) & (T[-1] >= 0), Z[-1], z0)
out['Z0'] = ('p', z0)
def _tlevel(tc_):
z_ = np.full(T[0].shape, np.nan)
for i in range(len(lv) - 1):
den = np.where(T[i] - T[i + 1] == 0, np.nan, T[i] - T[i + 1])
cross = np.isnan(z_) & (T[i] >= tc_) & (T[i + 1] < tc_)
z_ = np.where(cross, Z[i] + (T[i] - tc_) / den * (Z[i + 1] - Z[i]), z_)
z_ = np.where(np.isnan(z_) & (T[0] < tc_), 0.0, z_)
return np.where(np.isnan(z_) & (T[-1] >= tc_), Z[-1], z_)
out['Z3'] = ('p', _tlevel(3.0))
for L in (975, 950, 925):
if L in Zd and f'r{L}' in A:
zl = Zd[L] + 125.0 * (A[f't{L}'] - td_from_t_rh(A[f't{L}'], A[f'r{L}']))
out['WCD'] = ('p', np.maximum(z0 - zl, 0.0)); break
guard('Z0/WCD', _zero)
def _shear():
if all(k in A for k in ('u925', 'v925', 'u700', 'v700')):
out['SHR925_700'] = ('p', np.hypot(A['u925'] - A['u700'], A['v925'] - A['v700']))
if all(k in A for k in ('u925', 'v925', 'u500', 'v500')):
out['SHR925_500'] = ('p', np.hypot(A['u925'] - A['u500'], A['v925'] - A['v500']))
guard('シア', _shear)
def _winds():
for L in (925, 850):
if f'u{L}' in A and f'v{L}' in A:
out[f'WS{L}'] = ('p', np.hypot(A[f'u{L}'], A[f'v{L}']))
if 'u850' in A and 'v850' in A:
out['OV_U850'] = ('p', A['u850']); out['OV_V850'] = ('p', A['v850'])
lv = [L for L in (850, 700, 500, 300) if f'u{L}' in A and f'v{L}' in A]
if len(lv) < 3: return
ps = np.array(lv, float); w = np.zeros(len(ps))
for i in range(len(ps) - 1):
d_ = ps[i] - ps[i + 1]; w[i] += d_ / 2.0; w[i + 1] += d_ / 2.0
w = w / w.sum() # 圧力で重みづけした雲層平均風
uc = sum(wi * A[f'u{L}'] for wi, L in zip(w, lv)); vc = sum(wi * A[f'v{L}'] for wi, L in zip(w, lv))
out['CLDWIND'] = ('p', np.hypot(uc, vc)); out['OV_UCLD'] = ('p', uc); out['OV_VCLD'] = ('p', vc)
for L in (850, 925): # 下層ジェットの代表面
if f'u{L}' in A and f'v{L}' in A:
ul, vl = A[f'u{L}'], A[f'v{L}']; break
else:
return
up_, vp_ = uc - ul, vc - vl # Corfidi: 伝播ベクトル = 雲層平均風 − 下層風
um, vm = uc + up_, vc + vp_ # Corfidi: MCS移動 = 雲層平均風 + 伝播
out['PROPAG'] = ('p', np.hypot(up_, vp_)); out['OV_UPRP'] = ('p', up_); out['OV_VPRP'] = ('p', vp_)
out['MCSMOVE'] = ('p', np.hypot(um, vm)); out['OV_UMCS'] = ('p', um); out['OV_VMCS'] = ('p', vm)
guard('風・雲の動き', _winds)
def _moist():
lv = _levels(A, ('t', 'r'), [975, 950, 925, 850, 700, 600, 500])
if len(lv) < 3:
log.warning(f' 湿潤層/DGZ: 必要な気圧面が不足 {lv}'); return
Zd = level_heights(A, lv)
T = [A[f't{L}'] for L in lv]; R = [A[f'r{L}'] for L in lv]; Z = [Zd[L] for L in lv]
shp = T[0].shape; moist = np.zeros(shp); dgz = np.zeros(shp); top = np.zeros(shp); NS = 8
cont = np.zeros(shp); cold = np.zeros(shp); alive = np.ones(shp, bool)
for i in range(len(lv) - 1):
dz = np.maximum(Z[i + 1] - Z[i], 0.0) / NS
for k in range(NS):
f = (k + 0.5) / NS
Ti = T[i] + (T[i + 1] - T[i]) * f; Ri = R[i] + (R[i + 1] - R[i]) * f; Zi = Z[i] + (Z[i + 1] - Z[i]) * f
m = Ri >= 80.0
moist += np.where(m, dz, 0.0)
alive &= m
cont += np.where(alive, dz, 0.0)
cold += np.where(m & (Ti < 0.0), dz, 0.0)
dgz += np.where(m & (Ti <= -12.0) & (Ti >= -18.0), dz, 0.0)
top = np.where(m, np.maximum(top, Zi), top)
out['MOISTTHICK'] = ('p', moist); out['DGZTHICK'] = ('p', dgz); out['TOP80'] = ('p', top)
out['CLDTHICK'] = ('p', cont); out['COLDCLD'] = ('p', cold)
guard('湿潤層/DGZ', _moist)
# --- SST を使う指標(P格子)
if sstP is not None:
out['SST'] = ('p', sstP)
def _dt():
if 't850' in A: out['DT850'] = ('p', sstP - A['t850'])
if 't500' in A: out['DT500'] = ('p', sstP - A['t500'])
guard('SST-T', _dt)
def _fetch():
for L in (925, 950, 975, 850):
if f'u{L}' in A and f'v{L}' in A: break
else:
return
if not all(k in A for k in ('t925', 't850', 't500')): return
dist, heat, d850u, d500u = upstream_fetch(A[f'u{L}'], A[f'v{L}'], sstP, A['t925'], A['t850'], A['t500'], lon, lat)
out['FETCH'] = ('p', dist); out['HEATPATH'] = ('p', heat)
out['DT850_UP'] = ('p', d850u); out['DT500_UP'] = ('p', d500u)
guard('吹走距離', _fetch)
# --- 供給量−予報降水
def _diff():
if S is not None and 'precip' in S['a'] and 'MFC' in out:
pr = regrid(S['a']['precip'], S['lon'], S['lat'], lon, lat)
out['DIFF_SUP'] = ('p', out['MFC'][1] - pr)
guard('供給-降水', _diff)
if 'u925' in A and 'v925' in A:
out['OV_U925'] = ('p', A['u925']); out['OV_V925'] = ('p', A['v925'])
if S is not None:
B = S['a']; slon = S['lon']; slat = S['lat']
def _div10():
if 'u10' in B and 'v10' in B: out['DIV10'] = ('s', smooth(div_field(B['u10'], B['v10'], slon, slat) * 1e5, 1.5))
guard('DIV10', _div10)
if sstS is not None:
def _flux():
if all(k in B for k in ('u10', 'v10', 't2m')):
sp = np.hypot(B['u10'], B['v10'])
out['SENS'] = ('s', sp * (sstS - B['t2m']))
if 'rh2m' in B:
qs = EPS * es_hpa(sstS) / (1010.0 - es_hpa(sstS)) * 1000.0
q2 = q_from_t_rh(B['t2m'], B['rh2m'], 1010.0) * 1000.0
out['LATENT'] = ('s', sp * (qs - q2))
guard('海面フラックス', _flux)
if 'u10' in B and 'v10' in B:
out['OV_U10'] = ('s', B['u10']); out['OV_V10'] = ('s', B['v10'])
# --- スコア
if P is not None:
def get_p(key):
if key not in out: return None
g, a = out[key]
if g == 'p': return a
return regrid(a, S['lon'], S['lat'], P['lon'], P['lat'])
raws = {}
def _score():
for sk, okey in (('summer', 'SCORE_SUMMER'), ('winter', 'SCORE_WINTER')):
sc = thr.get('scores', {}).get(sk)
if not sc: continue
tot = None
for it in sc['items']:
ind = thr['indicators'].get(it['key']); arr = get_p(it['key'])
if ind is None or arr is None: continue
p = pts_arr(arr, ind) * float(it.get('weight', 1))
tot = p if tot is None else tot + p
if tot is None: continue
shown = tot; g = sc.get('gate')
if g:
ga = get_p(g['key'])
if ga is not None:
ok = (ga <= g['value']) if g.get('op', 'le') == 'le' else (ga >= g['value'])
shown = np.where(np.isfinite(ga) & ok, tot, np.nan)
raws[sk] = tot; out[okey] = ('p', shown)
if raws:
out['_TABLE'] = ('x', json.dumps(build_table(thr, out, P, get_p, raws), ensure_ascii=False))
guard('スコア', _score)
def _ptvals():
grids = {}
if P is not None: grids['p'] = P
if S is not None: grids['s'] = S
x0, x1, y0, y1 = AREAS['札幌近郊']; idx = {}
for g_, G_ in grids.items():
j = int(np.abs(G_['lat'] - SAPPORO[1]).argmin()); i = int(np.abs(G_['lon'] - SAPPORO[0]).argmin())
mx = (G_['lon'] >= x0) & (G_['lon'] <= x1); my = (G_['lat'] >= y0) & (G_['lat'] <= y1)
idx[g_] = (j, i, np.ix_(my, mx) if (mx.any() and my.any()) else None)
pt = {}; bx = {}
for key, (g_, arr) in out.items():
if key.startswith(('OV_', '_')) or g_ not in idx: continue
arr = np.asarray(arr, float); j, i, box = idx[g_]
pt[key] = _fl(arr[j, i])
if box is None: bx[key] = None; continue
sub = arr[box]
if not np.isfinite(sub).any(): bx[key] = None; continue
ind = thr['indicators'].get(key)
if key in ('SCORE_SUMMER', 'SCORE_WINTER'): bx[key] = _fl(np.nanmax(sub))
elif ind: bx[key] = _fl(np.nanmax(sub) if ind['op'] == 'ge' else np.nanmin(sub)) # 危険側
else: bx[key] = _fl(np.nanmean(sub))
out['_PTVALS'] = ('x', json.dumps({'has_p': P is not None, 'pt': pt, 'box': bx}, ensure_ascii=False))
guard('点値', _ptvals)
return out
# ---------------------------------------------------------------- SST (MGDSST) 取得・解析
def is_ssl_error(e):
r = getattr(e, 'reason', None)
return isinstance(e, ssl.SSLError) or isinstance(r, ssl.SSLError) or 'ASN1' in str(e) or 'CERTIFICATE' in str(e).upper()
SSL_HINT = ('→ Windowsの証明書ストアの読み込みに失敗している可能性があります。'
'「SSL証明書の検証を行わない」にチェックして再取得するか、ブラウザで取得したSSTファイルを「SSTファイルを指定…」で読み込んでください')
def describe_err(e):
"""例外を、画面に出せる短い日本語にする"""
if is_ssl_error(e):
return f'SSL/証明書エラー({getattr(e, "reason", None) or e}) {SSL_HINT}'
if isinstance(e, urllib.error.HTTPError): return f'HTTP {e.code} {e.reason}'
if isinstance(e, urllib.error.URLError): return f'接続エラー({e.reason})'
if isinstance(e, TimeoutError) or 'timed out' in str(e).lower(): return f'タイムアウト({e})'
return f'{type(e).__name__}: {e}'
def is_network_error(e):
"""HTTPで応答が返った(404など)ものと、SSLの問題を除く、通信そのものの失敗"""
if is_ssl_error(e): return False
return not isinstance(e, urllib.error.HTTPError) and isinstance(e, (urllib.error.URLError, TimeoutError, OSError))
_CTX_CACHE = {}
def make_ssl_context(insecure=False):
"""SSLコンテキストを作る。ssl.create_default_context() が証明書ストアの不具合で失敗する環境(Windows)に備え、
順に代替手段を試す。結果(失敗も)はキャッシュする。
insecure=True : 検証しないコンテキスト(システムの証明書を読まない)
通常 : 1) 標準 → 2) Windowsの証明書を1件ずつ読み、壊れた証明書だけ飛ばす → 3) certifi の証明書"""
key = bool(insecure)
if key in _CTX_CACHE:
c = _CTX_CACHE[key]
if isinstance(c, Exception): raise c
return c
try:
if insecure:
ctx = ssl._create_unverified_context()
else:
ctx = None; first = None
try:
ctx = ssl.create_default_context()
except Exception as e:
first = e
if ctx is None and hasattr(ssl, 'enum_certificates'): # Windows
try:
c2 = ssl.SSLContext(ssl.PROTOCOL_TLS_CLIENT); c2.check_hostname = True; c2.verify_mode = ssl.CERT_REQUIRED
n_ok = 0
for store in ('ROOT', 'CA'):
try: certs = ssl.enum_certificates(store)
except Exception: continue
for cert, enc, trust in certs:
if enc != 'x509_asn': continue
try: c2.load_verify_locations(cadata=cert); n_ok += 1
except Exception: pass # 壊れた証明書は飛ばす
if n_ok > 0: ctx = c2
except Exception:
ctx = None
if ctx is None:
try:
import certifi
c3 = ssl.SSLContext(ssl.PROTOCOL_TLS_CLIENT); c3.check_hostname = True; c3.verify_mode = ssl.CERT_REQUIRED
c3.load_verify_locations(cafile=certifi.where()); ctx = c3
except Exception:
ctx = None
if ctx is None: raise first if first is not None else ssl.SSLError('SSLコンテキストを作れません')
except Exception as e:
_CTX_CACHE[key] = e; raise
_CTX_CACHE[key] = ctx
return ctx
def http_get(url, timeout=60, insecure=False, retries=2):
last = None
for k in range(retries + 1):
try:
try: ctx = make_ssl_context(insecure)
except Exception:
if url.lower().startswith('https'): raise # httpsなのに安全な接続を作れない
ctx = None # httpなら証明書は不要
opener = urllib.request.build_opener(urllib.request.HTTPSHandler(context=ctx)) if ctx is not None \
else urllib.request.build_opener()
req = urllib.request.Request(url, headers={'User-Agent': 'Mozilla/5.0 (idx_engine)'})
with opener.open(req, timeout=timeout) as r:
return r.read()
except urllib.error.HTTPError as e:
if e.code in (400, 401, 403, 404, 410): raise # 再試行しても無駄
last = e
except Exception as e:
last = e
if is_ssl_error(e): raise # SSLの問題は再試行しても直らない
if k < retries: time.sleep(1.5 * (k + 1))
raise last
def parse_mgdsst(raw):
"""MGDSST(テキスト, 0.25度, 720行×1440列, 各3桁, 0.1℃単位, 999=陸, 888=海氷)を解析。"""
if raw[:2] == b'\x1f\x8b':
raw = gzip.decompress(raw)
elif raw[:2] == b'\x1f\x9d':
raise RuntimeError('.Z(compress)形式は未対応です。gzip形式か展開済みテキストを --sst-file で指定してください')
lines = [l for l in raw.decode('ascii', 'ignore').splitlines() if l.strip()]
if len(lines) < 721:
raise RuntimeError(f'MGDSSTの行数が想定外です: {len(lines)}')
rows = lines[1:721]
arr = np.empty((720, 1440), dtype=np.int32)
for j, l in enumerate(rows):
if len(l) < 4320: raise RuntimeError(f'MGDSSTの行長が想定外です: {len(l)}')
arr[j] = [int(l[i:i + 3]) for i in range(0, 4320, 3)]
sst = arr / 10.0
sst[arr == 999] = np.nan; sst[arr == 888] = -1.8
nvalid = int(np.isfinite(sst).sum())
if nvalid < 100000 or np.nanmax(sst) > 45 or np.nanmin(sst) < -3:
raise RuntimeError(f'MGDSSTの値が想定外です(有効格子数 {nvalid}, 範囲 {np.nanmin(sst):.1f}〜{np.nanmax(sst):.1f})')
lat = 89.875 - 0.25 * np.arange(720); lon = 0.125 + 0.25 * np.arange(1440)
lat = lat[::-1]; sst = sst[::-1, :] # 昇順に
m_lon = (lon >= SST_CROP[0]) & (lon <= SST_CROP[1]); m_lat = (lat >= SST_CROP[2]) & (lat <= SST_CROP[3])
return lon[m_lon], lat[m_lat], sst[np.ix_(m_lat, m_lon)]
def save_sst(out_dir, lon, lat, sst, date_str):
for name in (f'SST_{date_str}.npz', 'SST_latest.npz'):
tmp = os.path.join(out_dir, '~tmp_' + name)
np.savez_compressed(tmp, lon=lon, lat=lat, sst=sst, date=np.array(date_str))
os.replace(tmp if tmp.endswith('.npz') else tmp + '.npz', os.path.join(out_dir, name))
def load_sst(out_dir):
p = os.path.join(out_dir, 'SST_latest.npz')
if not os.path.exists(p): return None
z = np.load(p)
return {'lon': z['lon'], 'lat': z['lat'], 'sst': z['sst'], 'date': str(z['date'])}
LIST_RE = re.compile(r'href\s*=\s*["\']([^"\']*mgd_sst_glb_D(\d{8})[^"\']*)["\']', re.I)
def sst_list_candidates(now, insecure, errs):
"""年別フォルダの一覧(HTML)から {日付: [URL,...]} を作る。戻り値: (候補, 通信不可か, httpsが使えるか)"""
cands = {}; net_err = 0; https_ok = True
for year in (now.year, now.year - 1):
for base in SST_BASES:
if not https_ok and base.startswith('https'): continue
url_dir = f'{base}/{year}/'
try:
html = http_get(url_dir, 20, insecure, retries=0).decode('utf-8', 'ignore')
except Exception as e:
errs.append(f'{url_dir} : {describe_err(e)}')
if is_ssl_error(e): https_ok = False
elif is_network_error(e):
net_err += 1
if net_err >= 2: return cands, True, https_ok
continue
for href, d in LIST_RE.findall(html):
if 'norm' in href.lower(): continue
u = urllib.parse.urljoin(url_dir, href)
lst = cands.setdefault(d, [])
if u not in lst: lst.append(u)
if cands: break
if cands: break
return cands, False, https_ok
def sst_guess_candidates(now, https_ok=True):
"""一覧が読めない場合の保険。ファイル名を直接指定して試す"""
cands = {}
bases = [b for b in SST_BASES if https_ok or not b.startswith('https')][:2]
for k in range(0, 10):
d = now - timedelta(days=k); ds = d.strftime('%Y%m%d')
for base in bases:
for pre in ('', 're_'):
cands.setdefault(ds, []).append(f'{base}/{d.year}/{pre}mgd_sst_glb_D{ds}.txt.gz')
return cands
def fetch_latest_sst(out_dir, log, insecure=False, max_try=8):
"""最新の日別MGDSSTを取得して保存する。戻り値: (成功したか, 画面表示用メッセージ)
一覧→新しい日付から順に試す→失敗したら古い日付へ。全部失敗しても手元のSSTはそのまま使う。"""
now = datetime.utcnow() + timedelta(hours=9)
cur = load_sst(out_dir)
errs = []
cands, net_down, https_ok = sst_list_candidates(now, insecure, errs)
if not cands and not net_down:
for e in errs[:3]: log.info(f' SST一覧: {e}')
log.info('SST一覧を読めないため、ファイル名を直接指定して試します')
cands = sst_guess_candidates(now, https_ok)
if not https_ok and not insecure:
# httpsが使えず、httpでも一覧が取れなかった場合は、原因をそのまま伝える
ssl_msg = next((m for m in errs if 'SSL' in m), None)
if ssl_msg and not cands:
msg = 'SSTを取得できませんでした: ' + ssl_msg
if cur: msg += f' / 手元のSST(解析日 {cur["date"]})を引き続き使用します'
log.warning(msg); return False, msg
if net_down:
msg = 'SSTを取得できませんでした: 気象庁サーバーに接続できません(' + (errs[-1] if errs else '') + ')'
if cur: msg += f' / 手元のSST(解析日 {cur["date"]})を引き続き使用します'
log.warning(msg); return False, msg
dates = sorted(cands, reverse=True)
if cur and dates and cur['date'] >= dates[0]:
msg = f'SSTは最新です(解析日 {cur["date"]})'; log.info(msg); return True, msg
fails = []; tried = 0; net_err = 0
for d in dates:
if cur and cur['date'] >= d: break
tried += 1
if tried > max_try: break
for url in cands[d]:
try:
log.info(f'SSTを取得: {url}')
lon, lat, sst = parse_mgdsst(http_get(url, 120, insecure))
save_sst(out_dir, lon, lat, sst, d)
msg = f'SSTを取得しました(解析日 {d})'; log.info(msg); return True, msg
except Exception as e:
m = f'{os.path.basename(url)} : {describe_err(e)}'; fails.append(m); log.warning(f' 失敗 {m}')
if is_ssl_error(e): net_err = 99; break
if is_network_error(e):
net_err += 1
if net_err >= 2: break
if net_err >= 2: break
msg = 'SSTを取得できませんでした: ' + (fails[-1] if fails else '取得できる新しいファイルがありません')
if cur: msg += f' / 手元のSST(解析日 {cur["date"]})を引き続き使用します'
log.warning(msg); return False, msg
def sst_status_text(out_dir):
cur = load_sst(out_dir)
if not cur: return 'SST: なし(SSTを使う指標は計算されません)'
try:
days = ((datetime.utcnow() + timedelta(hours=9)).date() - datetime.strptime(cur['date'], '%Y%m%d').date()).days
except Exception:
days = None
t = f'SST: 解析日 {cur["date"]}' + (f'({days}日前)' if days is not None else '')
if days is not None and days >= 7: t += ' ※古いデータです'
return t
def load_sst_file(path, out_dir, log):
with open(path, 'rb') as f: raw = f.read()
m = re.search(r'(\d{8})', os.path.basename(path))
lon, lat, sst = parse_mgdsst(raw)
save_sst(out_dir, lon, lat, sst, m.group(1) if m else datetime.utcnow().strftime('%Y%m%d'))
log.info(f'SSTファイルを読み込みました: {path}')
# ---------------------------------------------------------------- 出力・統計
def write_meta(out_dir, cfg):
els = []
for e in META_ELEMENTS:
e = dict(e)
ind = cfg.get('indicators', {}).get(e['key'])
if ind:
e['danger'] = dict(ind, basis_text=BASIS_TEXT.get(ind.get('basis', ''), ''))
sk = {'SCORE_SUMMER': 'summer', 'SCORE_WINTER': 'winter'}.get(e['key'])
if sk and sk in cfg.get('scores', {}):
sc = cfg['scores'][sk]; mp = int(round(2 * sum(float(it.get('weight', 1)) for it in sc['items'])))
e['levels'] = [i + 0.5 for i in range(mp)]
e['danger'] = {'op': 'ge', 'caution': sc['judge']['caution'], 'warning': sc['judge']['warning'], 'basis': 'C',
'basis_text': BASIS_TEXT['C'], 'note': f'総合点({mp}点満点)の判定ライン'}
g = sc.get('gate')
e['desc'] = (f"{sc.get('label', sk)}の総合点。各条件を 0点(基準未満)/1点(注意)/2点(警戒) で採点して合計(最大{mp}点)。"
"内訳は『スコア内訳』タブで確認できます。"
+ (f" ゲート: {g.get('label', g['key'])}を満たす場所だけ表示します。" if g else ""))
els.append(e)
text = json.dumps({'version': 2, 'sapporo': SAPPORO, 'areas': AREAS, 'basis_text': BASIS_TEXT, 'categories': CATEGORIES, 'elements': els},
ensure_ascii=False, indent=1)
p = os.path.join(out_dir, 'IDX_META.json')
try:
with open(p, 'r', encoding='utf-8') as f:
if f.read() == text: return
except Exception:
pass
with open(p, 'w', encoding='utf-8') as f: f.write(text)
def load_thresholds(out_dir, log=None):
p = os.path.join(out_dir, 'idx_thresholds.json')
def write_default():
with open(p, 'w', encoding='utf-8') as f: json.dump(DEFAULT_CFG, f, ensure_ascii=False, indent=1)
return json.loads(json.dumps(DEFAULT_CFG))
if not os.path.exists(p): return write_default()
try:
with open(p, 'r', encoding='utf-8') as f: cfg = json.load(f)
except Exception as e:
if log: log.warning(f'idx_thresholds.json を読めません({e})。既定値を使います(ファイルは変更しません)')
return json.loads(json.dumps(DEFAULT_CFG))
if not isinstance(cfg, dict) or 'indicators' not in cfg: # 旧形式(v1)
try: os.replace(p, os.path.join(out_dir, 'idx_thresholds_old_v1.json'))
except Exception: pass
if log: log.info('旧形式の idx_thresholds.json を idx_thresholds_old_v1.json に退避し、新形式で作り直しました')
return write_default()
ver = cfg.get('version', 1)
if ver < DEFAULT_CFG['version']: # 旧版の設定は退避して、新しい既定値で作り直す(指標の構成が変わったため)
try: os.replace(p, os.path.join(out_dir, f'idx_thresholds_old_v{ver}.json'))
except Exception: pass
if log: log.info(f'旧版(v{ver})の idx_thresholds.json を idx_thresholds_old_v{ver}.json に退避し、新しい既定値で作り直しました')
return write_default()
base = json.loads(json.dumps(DEFAULT_CFG))
base['indicators'].update(cfg.get('indicators', {}))
base['scores'].update(cfg.get('scores', {}))
return base
def save_idx(out_dir, name, out, F):
data = {}
for key, (grid, arr) in out.items():
data[key] = np.array(arr) if key in ('_TABLE', '_PTVALS') else np.asarray(arr, dtype=np.float32)
if F.get('P') is not None:
data['lon'] = F['P']['lon'].astype(np.float32); data['lat'] = F['P']['lat'].astype(np.float32)
if F.get('S') is not None:
data['lon_s'] = F['S']['lon'].astype(np.float32); data['lat_s'] = F['S']['lat'].astype(np.float32)
tmp = os.path.join(out_dir, '~tmp_' + name)
np.savez_compressed(tmp, **data)
os.replace(tmp if os.path.exists(tmp) else tmp + '.npz', os.path.join(out_dir, name))
def write_stats(out_dir, source, init, ft, out, F):
path = os.path.join(out_dir, 'idx_stats.csv')
new = not os.path.exists(path)
valid = datetime.strptime(init, '%Y%m%d%H%M%S') + timedelta(hours=9 + ft)
with open(path, 'a', newline='', encoding='utf-8-sig') as f:
w = csv.writer(f)
if new: w.writerow(['init_utc', 'source', 'ft', 'valid_jst', 'area', 'key', 'max', 'mean'])
for key, (grid, arr) in out.items():
if key.startswith(('OV_', '_')): continue
G_ = F['P'] if grid == 'p' else F['S']
if G_ is None: continue
for an, (x0, x1, y0, y1) in AREAS.items():
mx = (G_['lon'] >= x0) & (G_['lon'] <= x1); my = (G_['lat'] >= y0) & (G_['lat'] <= y1)
if not mx.any() or not my.any(): continue
sub = np.asarray(arr)[np.ix_(my, mx)]
if not np.isfinite(sub).any(): continue
w.writerow([init, source, ft, valid.strftime('%Y-%m-%d %H:%M'), an, key,
f'{np.nanmax(sub):.3f}', f'{np.nanmean(sub):.3f}'])
def suggest_thresholds(out_dir, log, area='札幌近郊', months=None, q_caution=90.0, q_warning=97.5):
"""idx_stats.csv に溜まった実際の分布から、閾値の提案(上位10%/2.5%など)と、現在の閾値の超過率を出す。
ge の指標は領域最大値、le の指標は領域平均値(CSVに最小値が無いため)で評価する。"""
path = os.path.join(out_dir, 'idx_stats.csv')
if not os.path.exists(path):
log.warning('idx_stats.csv がありません。エンジンを動かしてデータを溜めてから実行してください'); return None
cfg = load_thresholds(out_dir, log)
vals = {}
with open(path, newline='', encoding='utf-8-sig') as f:
for r in csv.DictReader(f):
if r.get('area') != area: continue
try:
if months and int(r['valid_jst'][5:7]) not in months: continue
ind = cfg['indicators'].get(r['key'])
if ind is None: continue
vals.setdefault(r['key'], []).append(float(r['max'] if ind['op'] == 'ge' else r['mean']))
except Exception:
continue
mtxt = f'{",".join(str(m) for m in months)}月' if months else '全期間'
log.info(f'[閾値の提案] 領域={area} / 期間={mtxt} / 注意=上位{100 - q_caution:g}% 警戒=上位{100 - q_warning:g}% の値')
sug = {}; lacking = []
for key, ind in cfg['indicators'].items():
v = vals.get(key)
if not v or len(v) < 30:
lacking.append(key); continue
a = np.asarray(v)
if ind['op'] == 'ge':
c = float(np.percentile(a, q_caution)); w = float(np.percentile(a, q_warning))
ec = float((a >= ind['caution']).mean()); ew = float((a >= ind['warning']).mean())
else:
c = float(np.percentile(a, 100 - q_caution)); w = float(np.percentile(a, 100 - q_warning))
ec = float((a <= ind['caution']).mean()); ew = float((a <= ind['warning']).mean())
sug[key] = {'op': ind['op'], 'caution': round(c, 2), 'warning': round(w, 2), 'n': len(v),
'now_caution': ind['caution'], 'now_warning': ind['warning'],
'now_exceed_caution_pct': round(ec * 100, 1), 'now_exceed_warning_pct': round(ew * 100, 1)}
flag = ' ←常時近く超過' if ec > 0.5 else ''
log.info(f' {key:12s} n={len(v):5d} 現在 注意{ind["caution"]:g}/警戒{ind["warning"]:g}(超過率 {ec*100:.0f}%/{ew*100:.0f}%)'
f' → 提案 {c:.2f}/{w:.2f}{flag}')
if lacking: log.info(' データ不足(30件未満)で提案なし: ' + ', '.join(lacking))
outp = os.path.join(out_dir, 'idx_thresholds_suggest.json')
with open(outp, 'w', encoding='utf-8') as f:
json.dump({'area': area, 'months': months, 'q_caution': q_caution, 'q_warning': q_warning, 'suggest': sug}, f, ensure_ascii=False, indent=1)
log.info(f'提案を {outp} に保存しました(idx_thresholds.json は変更していません)')
return sug
# ---------------------------------------------------------------- メイン処理
def scan(npz_dir, models, hours, all_inits, stat=None):
now = datetime.utcnow(); items = []
st = stat if stat is not None else {}
st.update(total=0, name_ok=0, too_old=0, newest=None)
for fn in os.listdir(npz_dir):
if fn.lower().endswith('.npz'): st['total'] += 1
m = FILE_RE.match(fn)
if not m or m.group(1) not in models: continue
st['name_ok'] += 1
src, init, ft = m.group(1), m.group(2), int(m.group(3))
if st['newest'] is None or init > st['newest']: st['newest'] = init
if not all_inits:
try:
if (now - datetime.strptime(init, '%Y%m%d%H%M%S')).total_seconds() > hours * 3600:
st['too_old'] += 1; continue
except Exception: pass
items.append((src, init, ft, os.path.join(npz_dir, fn)))
items.sort(key=lambda t: (t[1], t[2]), reverse=True) # 新しい初期時刻から
return items
_last_diag = [None]
def diag(log, msg):
if msg != _last_diag[0]:
_last_diag[0] = msg; log.info(msg)
class EngineConfig:
"""エンジンの設定(GUI・コマンドライン共通)"""
def __init__(self, npz_dir, out_dir, models=('MSM', 'ANAL'), hours=12.0, all=False, force=False,
interval=60.0, sst_auto=True, sst_insecure=False, sst_file=None):
self.npz_dir = npz_dir; self.out_dir = out_dir; self.models = tuple(models); self.hours = float(hours)
self.all = bool(all); self.force = bool(force); self.interval = float(interval)
self.sst_auto = bool(sst_auto); self.sst_insecure = bool(sst_insecure); self.sst_file = sst_file
def run_cycle(cfg, sst, thr, log, force=False, stop=None, status=None):
n_done = 0; n_exist = 0; n_empty = 0; n_fail = 0
stat = {}
if not os.path.isdir(cfg.npz_dir):
diag(log, f'[診断] 入力フォルダが存在しません: {cfg.npz_dir}'); return 0
items = scan(cfg.npz_dir, cfg.models, cfg.hours, cfg.all, stat)
for src, init, ft, path in items:
if stop is not None and stop.is_set(): break
name = f'IDX_{src}_{init}_FT{ft:02d}.npz'
dst = os.path.join(cfg.out_dir, name)
if os.path.exists(dst) and not force: n_exist += 1; continue
if time.time() - os.path.getmtime(path) < 2: continue
try:
F = load_fields(path)
if F.get('P') is None and F.get('S') is None:
keys = sorted(np.load(path).files)[:20]
log.warning(f' 格子を認識できません {os.path.basename(path)} keys={keys}'); n_empty += 1; continue
out = compute_indices(F, sst, thr, log)
if not any(not k.startswith(('OV_', '_')) for k in out):
log.warning(f' 計算できた指標が0件です {os.path.basename(path)} '
f'(気圧面={"あり" if F.get("P") else "なし"}, 地上={"あり" if F.get("S") else "なし"})')
n_empty += 1; continue
save_idx(cfg.out_dir, name, out, F)
write_stats(cfg.out_dir, src, init, ft, out, F)
n_done += 1
log.info(f'出力: {name} ({len([k for k in out if not k.startswith(("OV_", "_"))])}要素)')
except Exception as e:
n_fail += 1
log.warning(f'処理失敗 {os.path.basename(path)}: {e}')
new_txt = (f'{stat["newest"][:8]} {stat["newest"][8:12]}UTC' if stat.get('newest') else 'なし')
rng = '全期間' if cfg.all else f'{cfg.hours:g}時間以内'
msg = (f'[診断] npz総数={stat.get("total")} / 対象モデル名一致={stat.get("name_ok")} / '
f'期間外({rng}の外)={stat.get("too_old")} / 処理対象={len(items)} / 出力済み={n_exist} / '
f'今回出力={n_done} / 失敗={n_fail} / 空={n_empty} / 最新の初期時刻={new_txt}'
+ ('' if stat.get('name_ok') else ' → 名前が MSM_YYYYMMDDHHMMSS_FTnn.npz 形式のファイルが見つかりません')
+ (' → すべて期間外です。「すべてのデータ」を選んでください' if stat.get('name_ok') and stat.get('too_old') == stat.get('name_ok') else ''))
diag(log, msg)
if status is not None:
status['summary'] = msg; status['last_cycle'] = time.time()
if n_done: status['last_output'] = time.time()
return n_done
def engine_loop(cfg, log, stop, status=None, once=False, sst_now=None):
"""メインループ。stop(Event)がセットされると終了する。sst_now(Event)がセットされるとSSTを今すぐ取得する。"""
status = status if status is not None else {}
_last_diag[0] = None
os.makedirs(cfg.out_dir, exist_ok=True)
thr = load_thresholds(cfg.out_dir, log); write_meta(cfg.out_dir, thr)
log.info(f'入力: {cfg.npz_dir} / 出力: {cfg.out_dir} / 対象: {",".join(cfg.models)} / '
f'範囲: {"すべてのデータ" if cfg.all else f"{cfg.hours:g}時間以内"}'
+ (' / 既存IDXも最初に1回だけ再計算' if cfg.force else ''))
if cfg.sst_file:
try: load_sst_file(cfg.sst_file, cfg.out_dir, log); status.update(sst_msg=f'SSTファイルを読み込みました: {cfg.sst_file}', sst_ok=True)
except Exception as e: log.warning(f'SSTファイル読み込み失敗: {e}'); status.update(sst_msg=f'SSTファイル読み込み失敗: {e}', sst_ok=False)
status['running'] = True
force = cfg.force
next_sst = 0.0
while not stop.is_set():
want = sst_now is not None and sst_now.is_set()
if want or (cfg.sst_auto and time.time() >= next_sst):
if sst_now is not None: sst_now.clear()
try: ok, msg = fetch_latest_sst(cfg.out_dir, log, cfg.sst_insecure)
except Exception as e: ok, msg = False, f'SST取得で予期しないエラー: {describe_err(e)}'; log.warning(msg)
status.update(sst_msg=msg, sst_ok=ok, sst_time=time.time())
next_sst = time.time() + (3 * 3600 if ok else 15 * 60) # 失敗時は15分後に再試行
sst = load_sst(cfg.out_dir)
if sst is None: diag(log, 'SSTなし: SSTを使う指標(DT850,DT500,吹走距離,海面フラックス等)はスキップします')
thr = load_thresholds(cfg.out_dir, log); write_meta(cfg.out_dir, thr)
try:
n = run_cycle(cfg, sst, thr, log, force, stop, status)
except Exception as e:
log.warning(f'サイクル失敗: {e}'); n = 0
force = False # 強制再計算は最初の1周だけ(毎周やり直さない)
if once: break
stop.wait(cfg.interval if n == 0 else 1.0)
status['running'] = False
# ---------------------------------------------------------------- GUI
def run_gui(args):
from PyQt6.QtWidgets import (QApplication, QWidget, QVBoxLayout, QHBoxLayout, QGridLayout, QLabel, QLineEdit,
QPushButton, QRadioButton, QCheckBox, QPlainTextEdit, QFileDialog, QGroupBox, QMessageBox)
from PyQt6.QtCore import QObject, pyqtSignal, QTimer, QSettings
class Bridge(QObject):
log = pyqtSignal(str)
class QtHandler(logging.Handler):
def __init__(self, bridge):
super().__init__(); self.bridge = bridge
def emit(self, record):
try: self.bridge.log.emit(self.format(record))
except Exception: pass
class Win(QWidget):
def __init__(self):
super().__init__()
self.setWindowTitle('特殊指数エンジン (idx_engine)'); self.resize(900, 720)
self.settings = QSettings('GPVTools', 'IdxEngine')
self.thread = None; self.stop = threading.Event(); self.sst_now = threading.Event(); self.status = {}
self.file_handler = None; self._sst_cache = (None, '')
self.bridge = Bridge(); self.bridge.log.connect(self.append_log)
self.log = logging.getLogger('idx'); self.log.setLevel(logging.INFO)
self.qh = QtHandler(self.bridge); self.qh.setFormatter(logging.Formatter('%(asctime)s %(message)s', '%H:%M:%S'))
self.log.addHandler(self.qh)
def pick(name, cli, default):
return cli if cli else str(self.settings.value(name, default))
def pbool(name, default):
v = self.settings.value(name, default)
return v if isinstance(v, bool) else str(v).lower() in ('true', '1')
lay = QVBoxLayout(self)
g = QGroupBox('フォルダ'); gl = QGridLayout(g)
self.ed_npz = QLineEdit(pick('npz_dir', args.npz_dir, '')); self.ed_out = QLineEdit(pick('out_dir', args.out_dir, ''))
b1 = QPushButton('参照…'); b2 = QPushButton('参照…')
gl.addWidget(QLabel('NPZフォルダ(本体エンジンが出力したGPVのnpz)'), 0, 0); gl.addWidget(self.ed_npz, 0, 1); gl.addWidget(b1, 0, 2)
gl.addWidget(QLabel('出力先フォルダ(IDX・SST・設定ファイル)'), 1, 0); gl.addWidget(self.ed_out, 1, 1); gl.addWidget(b2, 1, 2)
gl.setColumnStretch(1, 1)
b1.clicked.connect(lambda: self.browse(self.ed_npz)); b2.clicked.connect(lambda: self.browse(self.ed_out))
lay.addWidget(g)
g2 = QGroupBox('処理するデータの範囲(初期時刻が基準)'); hl = QHBoxLayout(g2)
self.rb_12 = QRadioButton('12時間以内のデータ'); self.rb_all = QRadioButton('すべてのデータ(件数が多いと時間がかかります)')
use_all = args.all or pbool('use_all', False)
self.rb_all.setChecked(use_all); self.rb_12.setChecked(not use_all)
hl.addWidget(self.rb_12); hl.addWidget(self.rb_all); lay.addWidget(g2)
g3 = QGroupBox('対象'); h3 = QHBoxLayout(g3)
models = args.models.split(',') if args.models else None
self.chk_msm = QCheckBox('MSM'); self.chk_msm.setChecked(('MSM' in models) if models else pbool('m_msm', True))
self.chk_anal = QCheckBox('ANAL'); self.chk_anal.setChecked(('ANAL' in models) if models else pbool('m_anal', True))
self.chk_force = QCheckBox('開始時に、作成済みのIDXも作り直す(最初の1回だけ。閾値を変えた時など)')
self.chk_force.setChecked(bool(args.force))
h3.addWidget(self.chk_msm); h3.addWidget(self.chk_anal); h3.addWidget(self.chk_force, 1); lay.addWidget(g3)
g4 = QGroupBox('SST(海面水温)'); v4 = QVBoxLayout(g4)
self.lab_sst = QLabel(''); self.lab_sst.setWordWrap(True)
row = QHBoxLayout()
self.chk_auto = QCheckBox('SSTを自動で取得する(気象庁MGDSST)'); self.chk_auto.setChecked(not args.no_sst_download and pbool('sst_auto', True))
self.chk_insecure = QCheckBox('SSL証明書の検証を行わない(社内ネットワーク等で取得に失敗する場合のみ)')
self.chk_insecure.setChecked(bool(args.sst_insecure) or pbool('sst_insecure', False))
row.addWidget(self.chk_auto); row.addWidget(self.chk_insecure, 1)
row2 = QHBoxLayout(); self.btn_sst = QPushButton('SSTを今すぐ取得'); self.btn_sstfile = QPushButton('SSTファイルを指定…')
row2.addWidget(self.btn_sst); row2.addWidget(self.btn_sstfile); row2.addStretch(1)
v4.addWidget(self.lab_sst); v4.addLayout(row); v4.addLayout(row2); lay.addWidget(g4)
self.btn_sst.clicked.connect(self.sst_now_clicked); self.btn_sstfile.clicked.connect(self.sst_file_clicked)
hb = QHBoxLayout()
self.btn_start = QPushButton('開始'); self.btn_stop = QPushButton('停止'); self.btn_stop.setEnabled(False)
self.lab_state = QLabel('停止中'); self.lab_state.setStyleSheet('font-weight:bold;')
hb.addWidget(self.btn_start); hb.addWidget(self.btn_stop); hb.addWidget(self.lab_state, 1); lay.addLayout(hb)
self.lab_sum = QLabel(''); self.lab_sum.setWordWrap(True); lay.addWidget(self.lab_sum)
self.txt = QPlainTextEdit(); self.txt.setReadOnly(True); self.txt.setMaximumBlockCount(3000); lay.addWidget(self.txt, 1)
self.btn_start.clicked.connect(self.start); self.btn_stop.clicked.connect(self.stop_clicked)
self.timer = QTimer(self); self.timer.timeout.connect(self.refresh); self.timer.start(1000)
self.refresh()
# ---- 操作
def append_log(self, s):
self.txt.appendPlainText(s)
def browse(self, edit):
d = QFileDialog.getExistingDirectory(self, 'フォルダを選択', edit.text() or os.getcwd())
if d: edit.setText(d)
def warn(self, msg):
QMessageBox.warning(self, '確認', msg)
def running(self):
return self.thread is not None and self.thread.is_alive()
def save_settings(self):
st = self.settings
st.setValue('npz_dir', self.ed_npz.text()); st.setValue('out_dir', self.ed_out.text())
st.setValue('use_all', self.rb_all.isChecked()); st.setValue('m_msm', self.chk_msm.isChecked())
st.setValue('m_anal', self.chk_anal.isChecked()); st.setValue('sst_auto', self.chk_auto.isChecked())
st.setValue('sst_insecure', self.chk_insecure.isChecked())
def start(self):
if self.running(): return
npz = self.ed_npz.text().strip(); out = self.ed_out.text().strip()
if not npz or not os.path.isdir(npz): self.warn('NPZフォルダが存在しません。「参照…」で選んでください。'); return
if not out: self.warn('出力先フォルダを指定してください。'); return
try: os.makedirs(out, exist_ok=True)
except Exception as e: self.warn(f'出力先フォルダを作成できません: {e}'); return
models = tuple(m for m, c in (('MSM', self.chk_msm), ('ANAL', self.chk_anal)) if c.isChecked())
if not models: self.warn('MSM か ANAL の少なくとも一方を選んでください。'); return
cfg = EngineConfig(npz, out, models, hours=12.0, all=self.rb_all.isChecked(), force=self.chk_force.isChecked(),
interval=60.0, sst_auto=self.chk_auto.isChecked(), sst_insecure=self.chk_insecure.isChecked())
self.save_settings()
self.attach_file_log(out)
self.stop.clear(); self.status = {}
self.thread = threading.Thread(target=self._run, args=(cfg,), daemon=True); self.thread.start()
self.refresh()
def _run(self, cfg):
try: engine_loop(cfg, self.log, self.stop, self.status, False, self.sst_now)
except Exception as e: self.log.warning(f'エンジンが異常終了しました: {e}')
finally: self.status['running'] = False
def stop_clicked(self):
self.stop.set(); self.btn_stop.setEnabled(False); self.lab_state.setText('停止しています…')
def attach_file_log(self, out):
try:
if self.file_handler is not None: self.log.removeHandler(self.file_handler); self.file_handler.close()
self.file_handler = logging.FileHandler(os.path.join(out, 'idx_engine.log'), encoding='utf-8')
self.file_handler.setFormatter(logging.Formatter('%(asctime)s %(message)s', '%H:%M:%S'))
self.log.addHandler(self.file_handler)
except Exception:
self.file_handler = None
def sst_now_clicked(self):
out = self.ed_out.text().strip()
if not out: self.warn('先に出力先フォルダを指定してください。'); return
if self.running():
self.sst_now.set(); self.log.info('SSTの取得を要求しました(現在の処理が終わり次第、実行します)'); return
try: os.makedirs(out, exist_ok=True)
except Exception as e: self.warn(f'出力先フォルダを作成できません: {e}'); return
insecure = self.chk_insecure.isChecked()
def work():
try: ok, msg = fetch_latest_sst(out, self.log, insecure)
except Exception as e: ok, msg = False, f'SST取得で予期しないエラー: {describe_err(e)}'; self.log.warning(msg)
self.status.update(sst_msg=msg, sst_ok=ok, sst_time=time.time())
threading.Thread(target=work, daemon=True).start()
def sst_file_clicked(self):
out = self.ed_out.text().strip()
if not out: self.warn('先に出力先フォルダを指定してください。'); return
path, _ = QFileDialog.getOpenFileName(self, 'MGDSSTファイルを選択(mgd_sst_glb_Dyyyymmdd.txt.gz など)', os.getcwd(),
'MGDSST (*.txt *.gz);;すべてのファイル (*)')
if not path: return
def work():
try:
os.makedirs(out, exist_ok=True); load_sst_file(path, out, self.log)
self.status.update(sst_msg=f'SSTファイルを読み込みました: {os.path.basename(path)}', sst_ok=True, sst_time=time.time())
except Exception as e:
self.log.warning(f'SSTファイル読み込み失敗: {e}')
self.status.update(sst_msg=f'SSTファイル読み込み失敗: {e}', sst_ok=False, sst_time=time.time())
threading.Thread(target=work, daemon=True).start()
# ---- 表示更新(1秒ごと)
def refresh(self):
run = self.running()
if not run and self.btn_stop.isEnabled() is False and self.lab_state.text() == '停止しています…':
self.lab_state.setText('停止中')
for w in (self.ed_npz, self.ed_out, self.rb_12, self.rb_all, self.chk_msm, self.chk_anal, self.chk_force):
w.setEnabled(not run)
self.btn_start.setEnabled(not run)
if run and self.lab_state.text() != '停止しています…':
self.btn_stop.setEnabled(True); self.lab_state.setText('稼働中(60秒ごとに新しいnpzを確認)')
elif not run:
self.btn_stop.setEnabled(False)
if self.lab_state.text() != '停止中': self.lab_state.setText('停止中')
out = self.ed_out.text().strip()
sst_line = ''
if out and os.path.isdir(out):
p = os.path.join(out, 'SST_latest.npz')
try: mt = os.path.getmtime(p)
except Exception: mt = None
if self._sst_cache[0] != (out, mt): self._sst_cache = ((out, mt), sst_status_text(out))
sst_line = self._sst_cache[1]
else:
sst_line = 'SST: 出力先フォルダを指定すると状態が表示されます'
msg = self.status.get('sst_msg')
if msg:
ts = time.strftime('%m/%d %H:%M', time.localtime(self.status.get('sst_time', time.time())))
sst_line += f'\n最後の取得結果({ts}): {msg}'
self.lab_sst.setText(sst_line)
self.lab_sum.setText(self.status.get('summary', ''))
def closeEvent(self, ev):
self.stop.set()
if self.thread is not None: self.thread.join(3)
super().closeEvent(ev)
app = QApplication(sys.argv)
w = Win(); w.show()
if args.autostart: QTimer.singleShot(300, w.start)
sys.exit(app.exec())
# ---------------------------------------------------------------- コマンドライン
def build_parser():
ap = argparse.ArgumentParser(description='札幌向け 特殊指数 補助エンジン(既定はGUI。--cli でコマンドライン実行)')
ap.add_argument('--npz-dir', default=None); ap.add_argument('--out-dir', default=None)
ap.add_argument('--models', default=None, help='例: MSM,ANAL')
ap.add_argument('--hours', type=float, default=12.0, help='この時間以内の初期時刻だけ処理(既定12)')
ap.add_argument('--all', action='store_true', help='時間制限なしで全初期時刻を処理')
ap.add_argument('--interval', type=float, default=60.0)
ap.add_argument('--once', action='store_true'); ap.add_argument('--force', action='store_true', help='作成済みIDXも最初の1回だけ作り直す')
ap.add_argument('--sst-file', default=None); ap.add_argument('--no-sst-download', action='store_true')
ap.add_argument('--sst-insecure', action='store_true', help='SSL証明書を検証しない(社内ネットワーク用)')
ap.add_argument('--cli', action='store_true', help='GUIを出さずに実行する')
ap.add_argument('--autostart', action='store_true', help='GUI起動と同時に処理を開始する')
ap.add_argument('--suggest-thresholds', action='store_true', help='idx_stats.csv から閾値の提案を作る(--out-dir 必須)')
ap.add_argument('--months', default=None, help='提案に使う月(例: 6,7,8,9)')
ap.add_argument('--area', default='札幌近郊', help='提案に使う領域名')
return ap
def main():
args = build_parser().parse_args()
if args.suggest_thresholds:
if not args.out_dir: print('--out-dir を指定してください'); return
log = logging.getLogger('idx'); log.setLevel(logging.INFO)
try: sys.stdout.reconfigure(encoding='utf-8', errors='replace')
except Exception: pass
log.addHandler(logging.StreamHandler(sys.stdout))
months = [int(m) for m in args.months.split(',')] if args.months else None
suggest_thresholds(args.out_dir, log, area=args.area, months=months); return
if not args.cli:
try:
import PyQt6 # noqa: F401
except ImportError:
print('PyQt6 が見つかりません。コマンドライン(--cli)で実行します。')
args.cli = True
if not args.cli:
run_gui(args); return
npz_dir = args.npz_dir or os.path.join(os.getcwd(), 'gpv_cache_npz')
out_dir = args.out_dir or os.path.join(os.getcwd(), 'idx_cache_npz')
models = tuple(m.strip() for m in (args.models or 'MSM,ANAL').split(',') if m.strip())
os.makedirs(out_dir, exist_ok=True)
log = logging.getLogger('idx'); log.setLevel(logging.INFO)
fmt = logging.Formatter('%(asctime)s %(message)s', '%H:%M:%S')
try: sys.stdout.reconfigure(encoding='utf-8', errors='replace')
except Exception: pass
sh = logging.StreamHandler(sys.stdout); sh.setFormatter(fmt); log.addHandler(sh)
fh = logging.FileHandler(os.path.join(out_dir, 'idx_engine.log'), encoding='utf-8'); fh.setFormatter(fmt); log.addHandler(fh)
cfg = EngineConfig(npz_dir, out_dir, models, hours=args.hours, all=args.all, force=args.force, interval=args.interval,
sst_auto=not args.no_sst_download, sst_insecure=args.sst_insecure, sst_file=args.sst_file)
stop = threading.Event()
try:
engine_loop(cfg, log, stop, {}, once=args.once)
except KeyboardInterrupt:
stop.set(); log.info('停止しました')
if __name__ == '__main__':
main()
未分類