# ==========================================
# 統合ワーカーコード
# ==========================================
WORKER_CODE = """
import sys, traceback
def main():
try:
import os, gc, json, subprocess, time, re
import numpy as np
import warnings; warnings.filterwarnings("ignore")
import cfgrib
os.environ["ECCODES_MAX_VALUES"] = "5000000"
if sys.platform == "win32":
conda_dir = os.path.dirname(sys.executable)
dll_paths = [os.path.join(conda_dir, "Library", "bin"), os.path.join(conda_dir, "bin")]
for p in dll_paths:
if os.path.exists(p):
os.environ["PATH"] = f"{p};{os.environ.get('PATH', '')}"
try: os.add_dll_directory(p)
except Exception: pass
def calculate_vorticity(u, v, lon, lat):
try:
R = 6371000.0; LON, LAT = np.meshgrid(lon, lat) if lon.ndim == 1 else (lon, lat)
rad_lat = np.deg2rad(LAT); rad_lon = np.deg2rad(LON)
dy = R * np.gradient(rad_lat, axis=0); dx = R * np.cos(rad_lat) * np.gradient(rad_lon, axis=1)
dx[dx == 0] = 1e-10; dy[dy == 0] = 1e-10
return ((np.gradient(v, axis=1) / dx) - (np.gradient(u, axis=0) / dy)) * 1e5
except Exception: return np.zeros_like(u)
model = sys.argv[1]; mode = sys.argv[2]; cache = sys.argv[3]
init = sys.argv[4]; target_fts_str = sys.argv[5]
f1 = sys.argv[6]; f2 = sys.argv[7]; wgrib2_path = sys.argv[8]
target_fts = [int(x) for x in target_fts_str.split(',')]
creationflags = 0x08000000 if sys.platform == "win32" else 0
# === 全モデル共通: wgrib2による高速スライス関数 ===
def slice_with_wgrib2(fin, ft_val, suffix):
if fin == "NONE" or not os.path.exists(fin): return "NONE"
fout = os.path.join(cache, f"temp_{model}_{init}_{ft_val}_{suffix}.bin")
if ft_val == 0:
match_str = r":(anl|0 hour [^:]+|[0-9]*[-]0 hour [^:]+):"
else:
patterns = [f"{ft_val} hour", f"[0-9]*[-]{ft_val} hour"]
if ft_val % 24 == 0:
days = ft_val // 24
patterns.extend([f"{days} day", f"[0-9]*[-]{days} day", f"{days} d", f"[0-9]*[-]{days} d"])
pat_joined = "|".join(patterns)
match_str = f":({pat_joined})[^:]*:"
cmd = [wgrib2_path, fin, "-match", match_str, "-grib", fout]
try:
subprocess.run(cmd, creationflags=creationflags, capture_output=True, text=True, timeout=120)
if os.path.exists(fout) and os.path.getsize(fout) > 0: return fout
except subprocess.TimeoutExpired:
print(f"WARN: slice timeout for FT={ft_val}")
except Exception: pass
return "NONE"
# ============================================================
# 【MEPS_GUID専用】 高速-match抽出 + P90/P10/超過確率計算
# ============================================================
if mode == "MEPS_GUID":
f_pall = f1; f_prrst = f2
d_p90 = {ft: {} for ft in target_fts}
d_p10 = {ft: {} for ft in target_fts}
d_prob1 = {ft: {} for ft in target_fts}
d_prob10 = {ft: {} for ft in target_fts}
d_prob20 = {ft: {} for ft in target_fts}
d_prob30 = {ft: {} for ft in target_fts}
for t_ft in target_fts:
# 巨大ファイルをFTごとに1回だけスライス（数十倍高速化）
f1_mini = slice_with_wgrib2(f_pall, t_ft, "pall")
f2_mini = slice_with_wgrib2(f_prrst, t_ft, "prrsf")
def extract_and_process(fin, src_name):
if fin == "NONE" or not os.path.exists(fin): return
# 軽量化されたminiファイルのインベントリを取得
res = subprocess.run([wgrib2_path, fin, "-s"], capture_output=True, text=True, creationflags=creationflags)
if not res.stdout.strip(): return
parsed_lines = []
for line in res.stdout.strip().split('\\n'):
if not line: continue
parts = line.split(':')
if len(parts) < 6: continue
param = parts[3]
time_str = parts[5]
m = re.search(r'([0-9]+)-([0-9]+)\\s+(hour|day)', time_str)
if not m:
m2 = re.search(r'([0-9]+)\\s+(hour|day)', time_str)
if m2:
end_val = int(m2.group(1)); unit = m2.group(2)
if unit == 'day': end_val *= 24
duration = 0
else: continue
else:
start_val = int(m.group(1)); end_val = int(m.group(2)); unit = m.group(3)
if unit == 'day': end_val *= 24; start_val *= 24
duration = end_val - start_val
mapped_ft = end_val
if mapped_ft != t_ft: continue
step_type = "accum"
if " max " in time_str: step_type = "max"
elif " min " in time_str: step_type = "min"
elif " ave " in time_str: step_type = "ave"
elif "inst" in time_str: step_type = "inst"
var_name = param.lower()
if param == "TPRATE": var_name = "tp"
elif param == "TSRATE": var_name = "asnow"
elif param == "TSTM": var_name = "tstm"
parsed_lines.append( (var_name, step_type, duration, line) )
grouped = {}
for meta in parsed_lines:
k = (meta[0], meta[1], meta[2])
if k not in grouped: grouped[k] = []
grouped[k].append(meta[3])
for k, group_lines in grouped.items():
var_name, step_type, duration = k
fixed_key = None
if var_name == "tstm": fixed_key = "thund"
elif var_name == "tp":
if src_name == "pall": fixed_key = "precip"
elif src_name == "prrsf":
if duration == 1: fixed_key = "precip1max"
elif duration == 3:
if step_type == "max": fixed_key = "precip1max"
else: fixed_key = "precip3max"
elif duration == 24: fixed_key = "precip24max"
elif var_name == "asnow":
if src_name == "prrsf" and duration in [3, 6, 12, 24]:
fixed_key = f"snow{duration}"
if not fixed_key:
expected_start = max(0, t_ft - (duration if duration else 3))
fixed_key = f"{var_name}_{src_name}_{step_type}_{expected_start}_{t_ft}"
inv_input = "\\n".join(group_lines) + "\\n"
temp_bin = os.path.join(cache, f"~tmp_ext_{model}_{init}_{t_ft}_{src_name}_{fixed_key}.bin")
try:
cmd_extract = [wgrib2_path, fin, "-i", "-grib", temp_bin]
process = subprocess.Popen(cmd_extract, stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE, creationflags=creationflags)
process.communicate(input=inv_input.encode('utf-8'), timeout=60)
except subprocess.TimeoutExpired:
process.kill()
continue
if os.path.exists(temp_bin) and os.path.getsize(temp_bin) > 0:
try:
dss = cfgrib.open_datasets(temp_bin, backend_kwargs={'indexpath': ''})
all_fields = []; lon_1d, lat_1d = None, None
for ds in dss:
if lon_1d is None:
lon = ds.coords['longitude'].values if 'longitude' in ds.coords else (ds.longitude.values if hasattr(ds, 'longitude') else None)
lat = ds.coords['latitude'].values if 'latitude' in ds.coords else (ds.latitude.values if hasattr(ds, 'latitude') else None)
if lon is not None and lat is not None:
lon_1d = lon[0, :] if lon.ndim == 2 else lon
lat_1d = lat[:, 0] if lat.ndim == 2 else lat
for v in ds.data_vars:
da = ds[v]
da_step = da.isel(step=0) if 'step' in da.dims else da
val = da_step.values.copy()
if 'number' in da_step.dims:
num_axis = da_step.dims.index('number')
val = np.moveaxis(val, num_axis, 0)
for i in range(val.shape[0]):
all_fields.append(val[i])
else:
all_fields.append(val)
for ds in dss: ds.close()
if all_fields and lon_1d is not None and lat_1d is not None:
stacked_val = np.stack(all_fields, axis=0)
grid_size = lon_1d.size * lat_1d.size
def set_coords(d_dict):
d_dict[f'lon_{fixed_key}'] = lon_1d.copy()
d_dict[f'lat_{fixed_key}'] = lat_1d.copy()
if grid_size > 100000 or 'lon' not in d_dict:
d_dict['lon'] = lon_1d.copy()
d_dict['lat'] = lat_1d.copy()
for d_dict in [d_p90[t_ft], d_p10[t_ft], d_prob1[t_ft], d_prob10[t_ft], d_prob20[t_ft], d_prob30[t_ft]]:
set_coords(d_dict)
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=RuntimeWarning)
d_p90[t_ft][fixed_key] = np.nanpercentile(stacked_val, 90, axis=0)
d_p10[t_ft][fixed_key] = np.nanpercentile(stacked_val, 10, axis=0)
valid_counts = np.sum(~np.isnan(stacked_val), axis=0)
with np.errstate(divide='ignore', invalid='ignore'):
d_prob1[t_ft][fixed_key] = np.where(valid_counts > 0, np.sum(stacked_val >= 1.0, axis=0) / valid_counts * 100.0, np.nan)
d_prob10[t_ft][fixed_key] = np.where(valid_counts > 0, np.sum(stacked_val >= 10.0, axis=0) / valid_counts * 100.0, np.nan)
d_prob20[t_ft][fixed_key] = np.where(valid_counts > 0, np.sum(stacked_val >= 20.0, axis=0) / valid_counts * 100.0, np.nan)
d_prob30[t_ft][fixed_key] = np.where(valid_counts > 0, np.sum(stacked_val >= 30.0, axis=0) / valid_counts * 100.0, np.nan)
except Exception as e: pass
finally:
gc.collect()
try: os.remove(temp_bin)
except: pass
extract_and_process(f1_mini, "pall")
extract_and_process(f2_mini, "prrsf")
try:
if f1_mini != "NONE" and os.path.exists(f1_mini): os.remove(f1_mini)
if f2_mini != "NONE" and os.path.exists(f2_mini): os.remove(f2_mini)
except: pass
if d_p90[t_ft]:
stat_dicts = {"MAX": d_p90[t_ft], "MIN": d_p10[t_ft], "PROB1": d_prob1[t_ft], "PROB10": d_prob10[t_ft], "PROB20": d_prob20[t_ft], "PROB30": d_prob30[t_ft]}
for stat_name, data_dict in stat_dicts.items():
if not data_dict: continue
final_filepath = os.path.join(cache, f"{model}_GUID_{stat_name}_{init}_FT{t_ft:02d}.npz")
temp_filepath = os.path.join(cache, f"~tmp_{model}_GUID_{stat_name}_{init}_FT{t_ft:02d}.npz")
try:
np.savez_compressed(temp_filepath, **data_dict)
success = False
for _ in range(10):
try:
os.replace(temp_filepath, final_filepath); success = True; break
except PermissionError: time.sleep(0.5)
if not success:
try: os.remove(final_filepath)
except: pass
np.savez_compressed(final_filepath, **data_dict)
try:
if os.path.exists(temp_filepath): os.remove(temp_filepath)
except: pass
except Exception: pass
print(f"SUCCESS:{t_ft}", flush=True)
return # MEPS_GUIDの処理はここで完結
# ============================================================
# 【GPV / ANAL / GUID (MSM,GSM,GSM_JP,ANAL)】 元のロジック(無変更)
# ============================================================
d_all = {ft: {} for ft in target_fts}
for t_ft in target_fts:
d = d_all[t_ft]
f1_mini = slice_with_wgrib2(f1, t_ft, "1")
f2_mini = slice_with_wgrib2(f2, t_ft, "2")
def process_file(filepath):
if filepath == "NONE": return
try:
dss = cfgrib.open_datasets(filepath, backend_kwargs={'indexpath': ''})
for ds in dss:
lon = ds.longitude.values if hasattr(ds, 'longitude') else None
lat = ds.latitude.values if hasattr(ds, 'latitude') else None
for v in ds.data_vars:
da = ds[v]
da_step = da.isel(step=0) if 'step' in da.dims else da
if 'number' in da_step.coords: da_step = da_step.isel(number=0)
sName = str(da.attrs.get('GRIB_shortName', v)).lower()
attrs_str = str(da.attrs).lower()
disc = da.attrs.get('GRIB_discipline', -1)
cat = da.attrs.get('GRIB_parameterCategory', -1)
num = da.attrs.get('GRIB_parameterNumber', -1)
is_upper = 'isobaricInhPa' in da_step.coords or 'level' in da_step.coords
def assign(k, val_to_assign, is_upper_flag):
if lon is None or lat is None:
d[k] = val_to_assign
return
lon_1d = lon[0, :] if lon.ndim == 2 else lon
lat_1d = lat[:, 0] if lat.ndim == 2 else lat
if k in ['thund', 'tstm'] and val_to_assign.shape != (len(lat_1d), len(lon_1d)):
target_lat_len, target_lon_len = val_to_assign.shape
lat_1d = np.linspace(lat_1d[0], lat_1d[-1], target_lat_len)
lon_1d = np.linspace(lon_1d[0], lon_1d[-1], target_lon_len)
lon_key = 'lon_pall' if is_upper_flag else 'lon_surf'
lat_key = 'lat_pall' if is_upper_flag else 'lat_surf'
if k in ['thund', 'tstm']:
lon_key = f'lon_{k}'; lat_key = f'lat_{k}'
if k not in d:
d[k] = val_to_assign.copy()
d[lon_key] = lon_1d.copy(); d[lat_key] = lat_1d.copy()
if not is_upper_flag and 'lon' not in d:
d['lon'] = lon_1d.copy(); d['lat'] = lat_1d.copy()
else:
old_lon = d[lon_key]; old_lat = d[lat_key]; old_val = d[k]
if old_val.shape == val_to_assign.shape and np.array_equal(old_lon, lon_1d):
d[k] = val_to_assign.copy()
return
new_lon = np.unique(np.concatenate([np.round(old_lon, 4), np.round(lon_1d, 4)]))
new_lat = np.unique(np.concatenate([np.round(old_lat, 4), np.round(lat_1d, 4)]))
new_lat = np.sort(new_lat)[::-1] if old_lat[0] > old_lat[-1] else np.sort(new_lat)
new_lon = np.sort(new_lon)
canvas = np.full((len(new_lat), len(new_lon)), np.nan)
lat_idx_old = np.where(np.isin(np.round(new_lat, 4), np.round(old_lat, 4)))[0]
lon_idx_old = np.where(np.isin(np.round(new_lon, 4), np.round(old_lon, 4)))[0]
if canvas[np.ix_(lat_idx_old, lon_idx_old)].shape == old_val.shape:
canvas[np.ix_(lat_idx_old, lon_idx_old)] = old_val
lat_idx_new = np.where(np.isin(np.round(new_lat, 4), np.round(lat_1d, 4)))[0]
lon_idx_new = np.where(np.isin(np.round(new_lon, 4), np.round(lon_1d, 4)))[0]
if canvas[np.ix_(lat_idx_new, lon_idx_new)].shape == val_to_assign.shape:
existing = canvas[np.ix_(lat_idx_new, lon_idx_new)]
mask = np.isnan(existing)
existing[mask] = val_to_assign[mask]
canvas[np.ix_(lat_idx_new, lon_idx_new)] = existing
d[k] = canvas; d[lon_key] = new_lon; d[lat_key] = new_lat
if k in ['thund', 'tstm'] or 'lon' not in d:
d['lon'] = new_lon; d['lat'] = new_lat
if mode in ["GPV", "ANAL"]:
if not is_upper:
val = da_step.values.copy()
while val.ndim > 2: val = val[0]
if sName in ['prmsl', 'msl', 'mslet']:
assign('slp', val / 100.0 if np.nanmax(val) > 2000 else val, is_upper)
elif sName in ['pres', 'sp'] or 'pressure' in attrs_str:
if 'slp' not in d and 'pres' not in d: assign('pres', val / 100.0 if np.nanmax(val) > 2000 else val, is_upper)
elif sName in ['10u', 'u', 'u10'] or ('u-component' in attrs_str and '10' in attrs_str):
assign('u10', val, is_upper)
elif sName in ['10v', 'v', 'v10'] or ('v-component' in attrs_str and '10' in attrs_str):
assign('v10', val, is_upper)
elif sName in ['2t', 't', 't2m', 'temp'] or ('temperature' in attrs_str and '2' in attrs_str):
assign('t2m', val - 273.15 if np.nanmax(val) > 150 else val, is_upper)
elif sName in ['2r', 'r', 'rh2m', 'rh'] or ('humidity' in attrs_str):
assign('rh2m', val, is_upper)
elif sName in ['tcc', 'hcc', 'mcc', 'lcc'] or (disc == 0 and cat == 6):
assign(sName if sName != 'unknown' else f"var_{disc}_{cat}_{num}", val, is_upper)
elif sName in ['tp', 'apcp', 'pr', 'precip'] or 'precip' in attrs_str or 'accum' in attrs_str:
if 'precip' not in d: assign('precip', np.nan_to_num(val, nan=0.0), is_upper)
else:
levels = []
if 'isobaricInhPa' in da_step.coords: levels = np.atleast_1d(da_step.isobaricInhPa.values)
elif 'level' in da_step.coords: levels = np.atleast_1d(da_step.level.values)
for l_idx, lvl in enumerate(levels):
lvl = int(lvl)
if lvl in [300, 500, 600, 700, 850, 925, 950, 975]:
if len(levels) > 1:
dim_n = 'isobaricInhPa' if 'isobaricInhPa' in da_step.coords else 'level'
val_l = da_step.isel(**{dim_n: l_idx}).values.copy()
else:
val_l = da_step.values.copy()
while val_l.ndim > 2: val_l = val_l[0]
if sName in ['t', 'temp'] or 'temperature' in attrs_str: assign(f't{lvl}', val_l - 273.15 if np.nanmax(val_l) > 150 else val_l, is_upper)
elif sName in ['u', 'u-component']: assign(f'u{lvl}', val_l, is_upper)
elif sName in ['v', 'v-component']: assign(f'v{lvl}', val_l, is_upper)
elif sName in ['r', 'rh', 'humidity']: assign(f'r{lvl}', val_l, is_upper)
elif sName in ['w', 'v-velocity', 'dz']: assign(f'w{lvl}', val_l, is_upper)
elif sName in ['gh', 'z', 'geopotential']: assign(f'gh{lvl}', val_l, is_upper)
elif mode == "GUID":
val = da_step.values.copy()
while val.ndim > 2: val = val[0]
if sName != 'unknown': assign(sName, val, is_upper)
else: assign(f"var_{disc}_{cat}_{num}", val, is_upper)
if sName in ['2t', 't', 't2m', 'tmp', 'temp'] or (disc == 0 and cat == 0 and num == 0) or 'temperature' in attrs_str:
assign('t2m', val - 273.15 if np.nanmax(val) > 150 else val, is_upper)
elif sName in ['2r', 'r', 'rh2m', 'rh'] or (disc == 0 and cat == 1 and num == 1) or 'humidity' in attrs_str:
assign('rh2m', val, is_upper)
elif sName in ['10u', 'u', 'u10', 'ugrd'] or (disc == 0 and cat == 2 and num == 2):
assign('u10', val, is_upper)
elif sName in ['10v', 'v', 'v10', 'vgrd'] or (disc == 0 and cat == 2 and num == 3):
assign('v10', val, is_upper)
elif sName in ['tp', 'apcp', 'pr', 'precip'] or (disc == 0 and cat == 1 and num in [8, 52]) or 'precip' in attrs_str or 'accum' in attrs_str:
if 'precip' not in d: assign('precip', np.nan_to_num(val, nan=0.0), is_upper)
elif sName in ['weasd', 'snod', 'snow', 'asnow'] or (disc == 0 and cat == 1 and num in [11, 13, 29, 60]) or 'snow' in attrs_str:
assign('snow', val, is_upper)
elif sName in ['wea', 'nswrs', 'nswrv', 'weather'] or (disc == 0 and cat == 19 and num == 192) or 'weather' in attrs_str:
assign('wea', val, is_upper)
elif sName in ['thund', 'lig', 'ltng', 'thunder', 'prstm', 'tstm'] or (disc == 0 and cat == 19 and num == 193) or 'thunder' in attrs_str:
assign('thund', val, is_upper)
for ds in dss: ds.close()
except Exception: pass
finally: gc.collect()
process_file(f1_mini)
process_file(f2_mini)
try:
if f1_mini != "NONE" and os.path.exists(f1_mini): os.remove(f1_mini)
if f2_mini != "NONE" and os.path.exists(f2_mini): os.remove(f2_mini)
except Exception: pass
for t_ft in target_fts:
d = d_all[t_ft]
if not d: continue
if mode in ["GPV", "ANAL"]:
for lvl in [300, 500, 600, 700, 850, 925, 950, 975]:
tc = d.get(f't{lvl}'); rh = d.get(f'r{lvl}')
if tc is not None and rh is not None:
rh_c = np.clip(rh, 0.1, 100)
e = 6.112 * np.exp((17.67*tc)/(tc+243.5)) * (rh_c/100.0)
td = (243.5*np.log(e/6.112))/(17.67-np.log(e/6.112))
d[f'tddep{lvl}'] = tc - td
tk = tc + 273.15; theta = tk*(1000.0/lvl)**0.2854; w = 0.622*e/(lvl-e)
d[f'ep{lvl}'] = theta * np.exp((2.5e6*w)/(1004.0*tk))
u500 = d.get('u500'); v500 = d.get('v500')
if u500 is not None and v500 is not None:
vort_lon = d.get('lon_pall') if 'lon_pall' in d else d.get('lon_surf')
vort_lat = d.get('lat_pall') if 'lat_pall' in d else d.get('lat_surf')
if vort_lon is not None and vort_lat is not None:
d['vort500'] = calculate_vorticity(u500, v500, vort_lon, vort_lat)
pfx = "GUID_" if mode == "GUID" else ""
final_filepath = os.path.join(cache, f"{model}_{pfx}{init}_FT{t_ft:02d}.npz")
temp_filepath = os.path.join(cache, f"~tmp_{model}_{pfx}{init}_FT{t_ft:02d}.npz")
try:
np.savez_compressed(temp_filepath, **d)
success = False
for _ in range(10):
try:
os.replace(temp_filepath, final_filepath); success = True; break
except PermissionError:
time.sleep(0.5)
if not success:
try: os.remove(final_filepath)
except Exception: pass
np.savez_compressed(final_filepath, **d)
try:
if os.path.exists(temp_filepath): os.remove(temp_filepath)
except Exception: pass
print(f"SUCCESS:{t_ft}", flush=True)
except Exception as e:
print(f"DEBUG_EXCEPTION: Error saving FT={t_ft} - {e}", flush=True)
except BaseException as e:
print(f"CRITICAL_ERROR: {traceback.format_exc()}", flush=True)
sys.exit(1)
if __name__ == '__main__': main()
"""
未分類
