#!/usr/bin/env python3 # # cust_smoke.py --- Smoke Example with Emissions File # See HYSPLIT Tutorial Sec. 14.1 # # Change history: # 24 Feb 2026 (BB) - standard-aligned, cross-platform version # 15 Jun 2026 (SZ) - create image files 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 "cust_smoke.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) OUT.mkdir(parents=True, exist_ok=True) 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") # constants wrk = TTR met_dir = wrk / "smoke" emit_file = met_dir / "EMITIMES.txt" # robust met file discovery preferred_met = met_dir / "WRF200704.bin" if preferred_met.exists(): met_file = preferred_met.name else: candidates = sorted(met_dir.glob("*.bin")) if not candidates: print(f"ERROR: No .bin meteorology files found in {met_dir}") print("Check your tutorial folder. Expected something like WRF200704.bin") sys.exit(1) met_file = candidates[0].name print(f"WARNING: WRF200704.bin not found. Using {met_file} instead.") # input file checks if not met_dir.exists(): print(f"ERROR: Smoke tutorial directory not found: {met_dir}") sys.exit(1) if not emit_file.exists(): print(f"ERROR: EMITIMES file not found: {emit_file}") sys.exit(1) if not (met_dir / "mayo_pm25.txt").exists(): print(f"ERROR: Observation file not found: {met_dir / 'mayo_pm25.txt'}") sys.exit(1) syr, smo, sda, shr = "07", "04", "16", "00" lat, lon, lvl = 31.15, -82.40, 100.0 run_hours = 96 top = 10000.0 # ASCDATA.CFG with open("ASCDATA.CFG", "w", encoding="ascii") as f: f.write( "-90.0 -180.0\n" "1.0 1.0\n" "180 360\n" "2\n" "0.2\n" f"'{PGM / 'bdyfiles'}{os.sep}'\n" ) # CONTROL with open("CONTROL", "w", encoding="ascii") as f: f.write(f"{syr} {smo} {sda} {shr}\n") f.write("1\n") f.write(f"{lat} {lon} {lvl}\n") f.write(f"{run_hours}\n") f.write("0\n") f.write(f"{top}\n") f.write("1\n") f.write(f"{met_dir}{os.sep}\n") f.write(f"{met_file}\n") f.write("1\n") f.write("PM25\n") f.write("1.0\n") f.write("72.0\n") f.write("00 00 00 18 00\n") f.write("1\n") f.write("0.0 0.0\n") f.write("0.10 0.10\n") f.write("20.0 20.0\n") f.write(f".{os.sep}\n") f.write("smoke.bin\n") f.write("1\n") f.write("100\n") f.write("00 00 00 00 00\n") f.write("00 00 00 00 00\n") f.write("00 01 00\n") f.write("1\n") f.write("0.0 0.0 0.0\n") f.write("0.0 0.0 0.0 0.0 0.0\n") f.write("0.0 0.0 0.0\n") f.write("0.0\n") f.write("0.0\n") # SETUP.CFG with open("SETUP.CFG", "w", encoding="ascii") as f: f.write("&SETUP\n") f.write(" initd = 3,\n") f.write(" khmax = 24,\n") f.write(" numpar = 10000,\n") f.write(" maxpar = 50000,\n") f.write(f" efile = '{emit_file}',\n") f.write("/\n") # run model Path("smoke.bin").unlink(missing_ok=True) run_cmd( [str(PGM / "exec" / f"hycs_std{exe_suffix}")], "Running hycs_std" ) # concplot Path("concplot.html").unlink(missing_ok=True) with open("LABELS.CFG", "w", encoding="ascii") as f: f.write("'TITLE&','### cust_smoke ### &'\n") run_cmd([ str(PGM / "exec" / f"concplot{exe_suffix}"), "+g1", "-ismoke.bin", f"-j{PGM / 'graphics' / 'arlmap'}", "-c4", "-z80", "-uug", "-n40:40", "-v500+200+100+50+20+10" ], "Generating concentration plot") html_file = OUT / "concplot.html" if html_file.exists(): if DSP == "YES": webbrowser.open(str(html_file.resolve())) else: convert_html_to_img(str(html_file), [("smoke_plot.png",1)]) if DSP == "YES": input("Press Enter to continue...") # station extraction with open("samplers.txt", "w", encoding="ascii") as f: f.write("Mayo 30.2556 -81.4533\n") run_cmd([ str(PGM / "exec" / f"con2stn{exe_suffix}"), "-ismoke.bin", "-ocon2stn.txt", "-r1", "-ssamplers.txt" ], "Extracting station concentrations") if Path("con2stn.txt").exists(): print(Path("con2stn.txt").read_text(encoding="ascii")) # timeplot Path("timeplot.html").unlink(missing_ok=True) run_cmd([ str(PGM / "exec" / f"timeplot{exe_suffix}"), "+g1", "-icon2stn.txt", f"-s{met_dir / 'mayo_pm25.txt'}" ], "Generating time series plot", allow_fail=True) # display html_file = OUT / "timeplot.html" if html_file.exists(): if DSP == "YES": webbrowser.open(str(html_file.resolve())) else: convert_html_to_img(str(html_file), [("smoke_times.png",1)])