#!/usr/bin/env python3 # # cust_emit.py --- Complex emission scenarios with the EMITIMES file # See HYSPLIT Tutorial Sec. 16.5 # # Change history: # 18 Apr 2025 (BB) - standard-aligned, cross-platform version # 16 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_emit.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(dat_name: str): env_override = os.environ.get("METDIR") candidates = [] if env_override: candidates.append(Path(env_override)) candidates += [ wrk / "volcano", TTR / "volcano", TTR, PGM / "tutorial" / "volcano", PGM / "tutorial", OUT / "volcano", PGM / "met", PGM / "meteorology", ] for base in candidates: if base and (base / dat_name).exists(): return base for root in [TTR, wrk, PGM]: try: for p in Path(root).rglob(dat_name): return p.parent except Exception: pass return None # ---------- Mirror BAT variables ---------- wrk = OUT.parent syr, smo, sda, shr = 10, 4, 14, 0 lat, lon, lvl = 63.63, -19.62, 100.0 hgt = 6000.0 run = 36 top = 32000.0 dat = "apr1420.bin" met = find_met_dir(dat) if met is None: print(f"ERROR: Meteorological file not found: {dat}") print("Set METDIR to the directory containing the file.") print("Example Windows: SET METDIR=C:\\Tutorial\\volcano") sys.exit(1) # Clean files for fname in [ "ASCDATA.CFG", "CONTROL", "CONTROL1", "SETUP.CFG", "LABELS.CFG", "EMITIMES", "volcano.bin", "parxplot.html", "concplot.html", ]: p = OUT / fname if p.exists(): p.unlink() # ASCDATA.CFG with open(OUT / "ASCDATA.CFG", "w", encoding="ascii") as f: f.write("-90.0 -180.0 lat/lon of lower left corner \n") f.write("1.0 1.0 lat/lon spacing in degrees \n") f.write("180 360 lat/lon number of data points \n") f.write("2 default land use category \n") f.write("0.2 default roughness length (m) \n") f.write(f"'{PGM / 'bdyfiles'}{os.sep}' directory of files \n") odir = f".{os.sep}" # CONTROL1 block control1_lines = [ f"{run}", "0", f"{top:.1f}", "1", f"{met}{os.sep}", dat, "4", "p006", "0.0", "0.0", "00 00 00 00 00", "p020", "0.0", "0.0", "00 00 00 00 00", "p060", "0.0", "0.0", "00 00 00 00 00", "p200", "0.0", "0.0", "00 00 00 00 00", "1", "45.0 0.0", "0.25 0.25", "90.0 360.0", odir, "volcano.bin", "1", "10000", "00 00 00 00 00", "00 00 00 00 00", "00 06 00", "4", "0.6 2.5 1.0", "0.0 0.0 0.0 0.0 0.0", "0.0 8.0E-05 8.0E-05", "0.0", "0.0", "2.0 2.5 1.0", "0.0 0.0 0.0 0.0 0.0", "0.0 8.0E-05 8.0E-05", "0.0", "0.0", "6.0 2.5 1.0", "0.0 0.0 0.0 0.0 0.0", "0.0 8.0E-05 8.0E-05", "0.0", "0.0", "20.0 2.5 1.0", "0.0 0.0 0.0 0.0 0.0", "0.0 8.0E-05 8.0E-05", "0.0", "0.0", ] # CONTROL with open(OUT / "CONTROL", "w", encoding="ascii") as f: f.write(f"{syr:02d} {smo:02d} {sda:02d} {shr:02d}\n") f.write("2\n") f.write(f"{lat:.2f} {lon:.2f} {lvl:.1f}\n") f.write(f"{lat:.2f} {lon:.2f} {hgt:.1f}\n") f.write("\n".join(control1_lines) + "\n") # SETUP.CFG with open(OUT / "SETUP.CFG", "w", encoding="ascii") as f: f.write("&SETUP\n") f.write("kmsl = 0,\n") f.write("ninit = 1,\n") f.write("ndump = 12,\n") f.write("ncycl = 12,\n") f.write("numpar = 5000,\n") f.write("maxpar = 50000,\n") f.write("efile = 'EMITIMES',\n") f.write("/\n") # EMITIMES emitimes_lines = [ "YYYY MM DD HH DURATION(hhhh) #RECORDS", "YYYY MM DD HH MM DURATION(hhmm) LAT LON HGT(m) RATE(/h) AREA(m2) HEAT(w)", "2010 04 14 00 0012 8", "2010 04 14 00 00 1200 63.63 -19.62 100.0 0.01 0.0 0.0", "2010 04 14 00 00 1200 63.63 -19.62 100.0 0.07 0.0 0.0", "2010 04 14 00 00 1200 63.63 -19.62 100.0 0.25 0.0 0.0", "2010 04 14 00 00 1200 63.63 -19.62 100.0 0.67 0.0 0.0", "2010 04 14 00 00 1200 63.63 -19.62 6000.0 0.01 0.0 0.0", "2010 04 14 00 00 1200 63.63 -19.62 6000.0 0.07 0.0 0.0", "2010 04 14 00 00 1200 63.63 -19.62 6000.0 0.25 0.0 0.0", "2010 04 14 00 00 1200 63.63 -19.62 6000.0 0.67 0.0 0.0", "2010 04 14 12 0012 8", "2010 04 14 12 00 1200 63.63 -19.62 100.0 0.01 0.0 0.0", "2010 04 14 12 00 1200 63.63 -19.62 100.0 0.07 0.0 0.0", "2010 04 14 12 00 1200 63.63 -19.62 100.0 0.25 0.0 0.0", "2010 04 14 12 00 1200 63.63 -19.62 100.0 0.67 0.0 0.0", "2010 04 14 12 00 1200 63.63 -19.62 12000.0 0.01 0.0 0.0", "2010 04 14 12 00 1200 63.63 -19.62 12000.0 0.07 0.0 0.0", "2010 04 14 12 00 1200 63.63 -19.62 12000.0 0.25 0.0 0.0", "2010 04 14 12 00 1200 63.63 -19.62 12000.0 0.67 0.0 0.0", "2010 04 15 00 0012 8", "2010 04 15 00 00 1200 63.63 -19.62 100.0 0.01 0.0 0.0", "2010 04 15 00 00 1200 63.63 -19.62 100.0 0.07 0.0 0.0", "2010 04 15 00 00 1200 63.63 -19.62 100.0 0.25 0.0 0.0", "2010 04 15 00 00 1200 63.63 -19.62 100.0 0.67 0.0 0.0", "2010 04 15 00 00 1200 63.63 -19.62 6000.0 0.01 0.0 0.0", "2010 04 15 00 00 1200 63.63 -19.62 6000.0 0.07 0.0 0.0", "2010 04 15 00 00 1200 63.63 -19.62 6000.0 0.25 0.0 0.0", "2010 04 15 00 00 1200 63.63 -19.62 6000.0 0.67 0.0 0.0", ] with open(OUT / "EMITIMES", "w", encoding="ascii") as f: f.write("\n".join(emitimes_lines) + "\n") # Run model run_cmd( [str(PGM / "exec" / f"hycs_std{exe_suffix}")], "Running hycs_std" ) # LABELS.CFG with open(OUT / "LABELS.CFG", "w", encoding="ascii") as f: f.write( f"'TITLE&','### {Path(sys.argv[0]).name}: " f"Eruption using EMITIMES ### &'\n" ) # Particle plot run_cmd([ str(PGM / "exec" / f"parxplot{exe_suffix}"), "+g1", "-iPARDUMP", f"-j{PGM / 'graphics' / 'arlmap'}", "-k1", "-z80", ], "Plotting particles") html_file = OUT / "parxplot.html" if html_file.exists(): if DSP == "YES": webbrowser.open(str(html_file.resolve())) else: convert_html_to_img(str(html_file), [("etime014.png",3)]) if DSP == "YES": input("Press Enter to continue...") # Concentration plot run_cmd([ str(PGM / "exec" / f"concplot{exe_suffix}"), "+g1", "-ivolcano.bin", f"-j{PGM / 'graphics' / 'arlmap'}", "-s0", "-z80", "-c4", "-uug", "-v100+50+20+10", "-x1.0E+16", ], "Plotting concentration") 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), [("etime013.png",2)]) print("\nScript completed successfully!")