#!/usr/bin/env python3 # # ens_phys.py --- Physics Ensemble # See HYSPLIT Tutorial Sec. 12.3 # # Change history: # 13 Nov 2025 (BB) - standard-aligned, cross-platform version # 09 Jan 2025 (SZ) - minor edits for style consistency # - create images when DSP=NO # 18 Jun 2026 (SZ) - standardize code for run_cmd(), convert_html_to_img(), # etc. across Python scripts # - ensure encoding is specified when opening files import os import sys import platform from pathlib import Path import subprocess import webbrowser if sys.version_info < (3, 8): print("Please use Python 3.8 or higher.") sys.exit(1) system = platform.system().lower() if system == "windows": PGM = Path(sys.argv[1]) if len(sys.argv) > 1 else Path("C:\\hysplit") OUT = Path(sys.argv[2]) if len(sys.argv) > 2 else PGM / "working" DSP = sys.argv[3] if len(sys.argv) > 3 else "YES" TTR = Path(os.environ.get("TTR", "C:\\Tutorial")) exe_suffix = ".exe" else: PGM = Path(sys.argv[1]) if len(sys.argv) > 1 else Path.home() / "hysplit" OUT = Path(sys.argv[2]) if len(sys.argv) > 2 else PGM / "working" DSP = sys.argv[3] if len(sys.argv) > 3 else "YES" TTR = Path(os.environ.get("TTR", Path.home() / "tutorial")) exe_suffix = "" if not PGM.exists(): script = Path(__file__).name if "__file__" in globals() else "ens_phys.py" print(f"ERROR: HYSPLIT not found at {PGM}") print(f""" Please specify the path to your HYSPLIT installation using a command-line argument: python {script} YOUR_HYSPLIT_DIRECTORY [OUT] [DSP (YES|NO)] Example: python {script} C:\\hysplit C:\\hysplit\\working YES """) sys.exit(1) if not TTR.exists(): print(f"ERROR: HYSPLIT tutorial not found at {TTR}") print(""" Please set the TTR environment variable to point to your HYSPLIT tutorial directory. For Windows users (Command Prompt): SET TTR=YOUR_HYSPLIT_TUTORIAL_DIRECTORY For Linux users (bash shell): export TTR='YOUR_HYSPLIT_TUTORIAL_DIRECTORY' For Linux users (csh or tcsh shell): setenv TTR 'YOUR_HYSPLIT_TUTORIAL_DIRECTORY' """) sys.exit(1) os.chdir(OUT) # ------------------------------------------------------------------ # Helper functions # ------------------------------------------------------------------ def run_cmd(cmd, desc="", allow_fail=False, stdout_file=None, append=False): if desc: print(f"{desc}...") exitcode = 0 with subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True) as p: if stdout_file: mode="a" if append else "w" with open(stdout_file, mode, encoding="utf-8") as f_out: for ln in p.stdout: print(ln, end="") f_out.write(ln) else: with p.stdout as f: for ln in iter(f.readline, b'\n'): if len(ln) == 0: break print(ln.rstrip()) done = False while not done: try: exitcode = p.wait(timeout=1) done = True except subprocess.TimeoutExpired: print("WARN: subprocess timed out.") if exitcode != 0: if allow_fail: print(f"WARN: command failed with exit code {exitcode} -> {' '.join(cmd)}") else: print(f"ERROR: command failed with exit code {exitcode} -> {' '.join(cmd)}") sys.exit(1) def convert_html_to_img(html_file, img_name_frame_nums): # Split SVG into Linux-consistent naming: F00-sec71_001.svg, etc. temp_img_name = "temp" run_cmd([str(PGM / "exec" / f"splitsvg{exe_suffix}"), f"-i{html_file}", f"-o{temp_img_name}.svg"]) # Locate ImageMagick for cross-platform conversion if system == "windows": convert_exe = None imgk_dir = next( (k for k in os.environ.get('PATH', '').split(os.pathsep) if "ImageMagick" in k), None ) if imgk_dir: # Use magick.exe for ImageMagick 7.x and later. Use convert.exe for earlier versions. for cmd in ('magick', 'convert'): candidate = os.path.join(imgk_dir, f'{cmd}{exe_suffix}') if os.path.exists(candidate): convert_exe = candidate break else: convert_exe = "convert" if convert_exe: # Match Linux naming for DSP=NO mode exactly for img_name, frame_num in img_name_frame_nums: svg = Path(f"F{frame_num:02d}-{temp_img_name}.svg") run_cmd([convert_exe, str(svg), img_name]) print(f"Converted {svg} -> {img_name}") svg.unlink() else: print("ImageMagick not found — skipping .svg conversion") def find_met_dir(): """ Try tutorial CAPTEX directory first, then fall back to PGM/captex. Mirrors %WRK%\\captex logic but more robust. """ candidates = [ TTR / "captex", PGM / "captex" ] for c in candidates: if c.exists(): return c print("ERROR: Could not find CAPTEX met directory in:") for c in candidates: print(f" - {c}") sys.exit(1) def write_ascdata_cfg(): asc_path = OUT / "ASCDATA.CFG" if asc_path.exists(): asc_path.unlink() lines = [ "-90.0 -180.0 lat/lon of lower left corner", "1.0 1.0 lat/lon spacing in degrees", "180 360 lat/lon number of data points", "2 default land use category", "0.2 default roughness length (m)", f"'{PGM / 'bdyfiles'}{os.sep}' directory of files" ] asc_path.write_text("\n".join(lines) + "\n", encoding="ascii") def write_control(met_dir_str, met_file, inp_name): """ Direct port of the CONTROL content from the 12.3 ens_phys batch script. """ control_path = OUT / "CONTROL" lines = [ "83 09 25 17", # SYR SMO SDA SHR "1", # Number of starting locations "39.90 -84.22 10.0", # LAT LON LVL "49", # RUN (hours) "0", # Vertical motion method "10000.0", # TOP "1", # Number of meteorological files f"{met_dir_str}{os.sep}",# MET path (directory) met_file, # MET filename "1", # Number of pollutants "PMCH", # Pollutant ID "67000.0", # Emission rate "3.0", # Emission hours "00 00 00 00 00", # Release start (offset) "1", # Number of concentration grids "42.0 -78.0", # Grid center lat/lon "0.25 0.25", # Grid spacing "15.0 25.0", # Grid span (x,y) f".{os.sep}", # Output path inp_name, # Output file base name (ensphys) "1", # Ensemble starting index "100", # Number of ensemble members "83 09 26 06 00", # Start of ensemble window "83 09 28 15 00", # End of ensemble window "00 03 00", # Time step "1", # Probability type "0.0 0.0 0.0", # Reserved params "0.0 0.0 0.0 0.0 0.0", "0.0 0.0 0.0", "0.0", "0.0" ] control_path.write_text("\n".join(lines) + "\n", encoding="ascii") def write_setup_cfg(member_code: str): """ Rebuild SETUP.CFG for each physics ensemble member (001–015), matching the IF logic from the batch file. FAST_RUN=YES -> lighter particle settings for dev testing. """ setup_path = OUT / "SETUP.CFG" if setup_path.exists(): setup_path.unlink() fast = os.environ.get("FAST_RUN", "NO").upper() == "YES" if fast: numpar = 1500 maxpar = 20000 else: numpar = 5000 maxpar = 10000 lines = [ "&SETUP", f" numpar = {numpar},", f" maxpar = {maxpar},", " initd = 0," ] # Map the 001–015 logic if member_code == "001": lines.append(" idsp = 2,") elif member_code == "002": lines.append(" kmixd = 1,") elif member_code == "003": lines.append(" kmixd = 3,") elif member_code == "004": lines.append(" kmix0 = 50,") elif member_code == "005": lines.append(" kzmix = 1,") elif member_code == "006": lines.append(" kdef = 1,") elif member_code == "007": lines.append(" kbls = 2,") elif member_code == "008": lines.append(" kblt = 1,") elif member_code == "009": lines.append(" kblt = 3,") elif member_code == "010": lines.append(" kblt = 5,") elif member_code == "011": lines.append(" vscales = 200.0,") elif member_code == "012": lines.append(" vscales = -1.0,") elif member_code == "013": lines.append(" kblt = 1,") lines.append(" vscales = -1.0,") elif member_code == "014": lines.append(" kblt = 3,") lines.append(" vscales = -1.0,") elif member_code == "015": lines.append(" kblt = 5,") lines.append(" vscales = -1.0,") lines.append("/") setup_path.write_text("\n".join(lines) + "\n", encoding="ascii") def write_labels_cfg(): labels_path = OUT / "LABELS.CFG" lines = [ "TITLE& ### ens_phys.py ### &", "UNITS& g&", "MAPID& Air Concentration&" ] labels_path.write_text("\n".join(lines) + "\n", encoding="ascii") def handle_html_output(html_name: str, img_name_frame_nums: list): html_path = OUT / html_name if DSP.upper() == "YES": if html_path.exists(): webbrowser.open(html_path.as_uri()) input("Press Enter to continue ...") else: print(f"WARNING: {html_path} does not exist to display.") else: convert_html_to_img(html_path, img_name_frame_nums) # --------------------------------------------------------------------------- # Main workflow – physical parameter ensemble (ens_phys) # --------------------------------------------------------------------------- def main(): print(f"Using HYSPLIT at: {PGM}") print(f"Working directory: {OUT}") print(f"Tutorial directory: {TTR}") print(f"Display mode (DSP): {DSP}") OUT.mkdir(parents=True, exist_ok=True) inp_name = "ensphys" map_dir = PGM / "graphics" / "arlmap" # 1. ASCDATA.CFG write_ascdata_cfg() print("ASCDATA.CFG written.") # 2. CONTROL met_dir = find_met_dir() met_file = "captex2_wrf27uw.bin" met_dir_str = str(met_dir) write_control(met_dir_str, met_file, inp_name) print("CONTROL file written.") # 3. Clear old ensphys.* files if present for f in OUT.glob(f"{inp_name}.*"): try: f.unlink() except Exception as exc: print(f"WARNING: Could not delete {f}: {exc}") # 4. Physics ensemble loop – recreate SETUP.CFG, run hycs_std, rename output hycs_std_exe = PGM / "exec" / f"hycs_std{exe_suffix}" member_codes = [ "001", "002", "003", "004", "005", "006", "007", "008", "009", "010", "011", "012", "013", "014", "015" ] env_ens = os.environ.get("ENSEMBLE_MAX") if env_ens: try: max_member = int(env_ens) except ValueError: print("WARNING: Invalid ENSEMBLE_MAX value; using full 15 members.") max_member = len(member_codes) max_member = max(1, min(max_member, len(member_codes))) codes_to_run = member_codes[:max_member] else: codes_to_run = member_codes print(f"Running physical ensemble members: {', '.join(codes_to_run)}") for code in codes_to_run: write_setup_cfg(code) run_cmd([str(hycs_std_exe)], desc=f"Running hycs_std for member {code}") base_file = OUT / inp_name # e.g. ensphys member_file = OUT / f"{inp_name}.{code}" # ensphys.001 etc. if base_file.exists(): if member_file.exists(): try: member_file.unlink() except Exception as exc: print(f"WARNING: Could not remove existing {member_file}: {exc}") try: base_file.rename(member_file) except Exception as exc: print(f"WARNING: Could not rename {base_file} to {member_file}: {exc}") else: print(f"WARNING: Expected output {base_file} not found after hycs_std.") # 5. Probability post-processing: conprob conprob_exe = PGM / "exec" / f"conprob{exe_suffix}" run_cmd([str(conprob_exe), "-bensphys", "-t1", "-p1", "-z1"], desc="Running conprob for physical ensemble probability") # 6. LABELS.CFG for plots write_labels_cfg() print("LABELS.CFG written.") # 7. concplot – 90% probability map concplot_exe = PGM / "exec" / f"concplot{exe_suffix}" concplot_html = OUT / "concplot.html" if concplot_html.exists(): concplot_html.unlink() run_cmd([ str(concplot_exe), "+g1", "-iprob90", f"-j{map_dir}", "-m0", "-k1", "-z80" ], desc="Generating concplot.html") handle_html_output("concplot.html", [("ephys003.png", 4)]) # 8. boxplots – time series boxplots at receptor point (n=10) boxplots_exe = PGM / "exec" / f"boxplots{exe_suffix}" boxplots_html = OUT / "boxplots.html" if boxplots_html.exists(): boxplots_html.unlink() run_cmd([ str(boxplots_exe), "+g1", "-n10", "-x-78.80", "-y42.25" ], desc="Generating boxplots.html") handle_html_output("boxplots.html", [("ephys004.png", 1)]) # 9. ensplots – physical ensemble spread plots (n=10) ensplots_exe = PGM / "exec" / f"ensplots{exe_suffix}" ensplots_html = OUT / "ensplots.html" if ensplots_html.exists(): ensplots_html.unlink() run_cmd([ str(ensplots_exe), "+g1", "-n10", "-bensphys", "-x-78.80", "-y42.25" ], desc="Generating ensplots.html") handle_html_output("ensplots.html", [("ephys005.png", 1)]) print("Physical ensemble probability and plot generation complete.") if __name__ == "__main__": main()