Add driver script to convert GSI bias coefficient files for UFO

ush/run_satbias_conv.py

@@ -0,0 +1,135 @@
+#!/usr/bin/env python3
+import os
+import datetime as dt
+import logging
+import subprocess
+import csv
+import shutil
+import glob
+import argparse
+import yaml
+
+logging.basicConfig(format='%(asctime)s:%(levelname)s:%(message)s',
+                    level=logging.INFO, datefmt='%Y-%m-%d %H:%M:%S')
+
+predictors = [
+    'constant',
+    'zenith_angle',
+    'cloud_liquid_water',
+    'lapse_rate_order_2',
+    'lapse_rate',
+    'cosine_of_latitude_times_orbit_node',
+    'sine_of_latitude',
+    'emissivity',
+    'scan_angle_order_4',
+    'scan_angle_order_3',
+    'scan_angle_order_2',
+    'scan_angle',
+]
+
+
+def run_satbias_conv(config):
+    # run ioda-converters satbias converter to take GSI satbias file
+    # and create UFO compatible input files
+    # get times from config
+    startTime = config['start time']
+    endTime = config['end time']
+    assim_freq = int(config['assim_freq'])
+    # get paths from config
+    gsi_bc_root = config['gsi_bc_root']
+    ufo_bc_root = config['ufo_bc_root']
+    work_root = config['work_root']
+    converter_exe = config['satbias2ioda']
+    cdump = config.get('dump', 'gdas')
+    # loop through all cycles
+    nowTime = startTime
+    while nowTime <= endTime:
+        cdate = nowTime.strftime("%Y%m%d%H")
+        pdy = nowTime.strftime("%Y%m%d")
+        cyc = nowTime.strftime("%H")
+        logging.info(f'Processing sat bias files for cycle: {cdate}')
+        # create working directory
+        workdir = os.path.join(work_root, f'{cdump}.{pdy}', cyc, 'atmos')
+        os.makedirs(workdir, exist_ok=True)
+        # link the GSI files to the working directory
+        cycle_in_dir = os.path.join(gsi_bc_root, f'{cdump}.{pdy}', cyc, 'atmos')
+        prefix = f"{cdump}.t{cyc}z"
+        abias_path = os.path.join(cycle_in_dir, prefix + '.abias')
+        abias_pc_path = os.path.join(cycle_in_dir, prefix + '.abias_pc')
+        orig_paths = [converter_exe, abias_path, abias_pc_path]
+        new_paths = [
+            os.path.join(workdir, 'satbias2ioda.x'),
+            os.path.join(workdir, 'satbias_crtm_in'),
+            os.path.join(workdir, 'satbias_crtm_pc'),
+        ]
+        for src, dest in zip(orig_paths, new_paths):
+            if os.path.exists(dest):
+                os.unlink(dest)
+            os.symlink(src, dest)
+        # open the text file to get a list of satellites/sensors to process
+        satlist = []
+        with open(new_paths[1]) as csvfile:
+            reader = csv.reader(csvfile)
+            for row in reader:
+                splitrow = row[0].split()
+                if splitrow[1] not in satlist:
+                    try:
+                        a = float(splitrow[1])
+                    except ValueError:
+                        if len(splitrow[1]) > 0:
+                            satlist.append(splitrow[1])
+        # loop through satellites/sensors to write tlapmean txt file
+        for sat in satlist:
+            outstr = ''
+            outfile = os.path.join(workdir, f'{sat}_tlapmean.txt')
+            with open(new_paths[1]) as csvfile:
+                reader = csv.reader(csvfile)
+                for row in reader:
+                    splitrow = row[0].split()
+                    if splitrow[1] == sat:
+                        outstr = outstr + f'{sat} {splitrow[2]} {splitrow[3]}\n'
+            with open(outfile, 'w') as f:
+                f.write(outstr)
+        # create YAML for input to converter
+        outyaml = os.path.join(workdir, 'satbias_converter.yaml')
+        with open(outyaml, 'w') as f:
+            f.write('input coeff file: satbias_crtm_in\n')
+            f.write('input err file: satbias_crtm_pc\n')
+            f.write('default predictors: &default_preds\n')
+            for pred in predictors:
+                f.write(f'- {pred}\n')
+            f.write('output:\n')
+            for sat in satlist:
+                f.write(f'- sensor: {sat}\n')
+                f.write(f'  output file: {sat}_satbias.nc4\n')
+                f.write('  predictors: *default_preds\n')
+        # run executable
+        runcmd = f'./satbias2ioda.x satbias_converter.yaml'
+        p = subprocess.Popen(runcmd, shell=True, cwd=workdir)
+        p.communicate()
+        # copy output files to output directory
+        outdir = os.path.join(ufo_bc_root, f'{cdump}.{pdy}', cyc, 'atmos')
+        os.makedirs(outdir, exist_ok=True)
+        ncfiles = glob.glob(os.path.join(workdir, '*.nc4'))
+        txtfiles = glob.glob(os.path.join(workdir, '*.txt'))
+        allfiles = ncfiles + txtfiles
+        for f in allfiles:
+            shutil.move(f, os.path.join(outdir, os.path.basename(f)))
+        # remove temp directory
+        shutil.rmtree(workdir)
+        # advance to the next cycle
+        nowTime = nowTime + dt.timedelta(hours=assim_freq)
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-c', '--config', type=str, help='Input YAML Configuration', required=True)
+    args = parser.parse_args()
+    # read YAML for config
+    try:
+        with open(args.config, 'r') as yaml_opened:
+            config = yaml.safe_load(yaml_opened)
+    except Exception as e:
+        logging.error(f'Error occurred when attempting to load: {args.config}, error: {e}')
+
+    run_satbias_conv(config)