#! /usr/bin/env python3 import os import datetime import numpy as np import subprocess import shutil import sys import netCDF4 as netcdf import numpy as np import glob import pandas as pd from time import sleep def run_shell_command(command): """! Run shell command Args: command - list of agrument entries (string) Returns: """ print("Running "+' '.join(command)) if any(mark in ' '.join(command) for mark in ['"', "'", '|', '*', '>']): run_command = subprocess.run( ' '.join(command), shell=True ) else: run_command = subprocess.run(command) if run_command.returncode != 0: print("FATAL ERROR: "+' '.join(run_command.args)+" gave return code " +str(run_command.returncode)) sys.exit(1) def metplus_command(conf_file_name): """! Write out full call to METplus Args: conf_file_name - METplus conf file name (string) Returns: metplus_cmd - full call to METplus (string) """ run_metplus = os.path.join(os.environ['METPLUS_PATH'], 'ush', 'run_metplus.py') machine_conf = os.path.join(os.environ['PARMevs'], 'metplus_config', 'machine.conf') conf_file = os.path.join(os.environ['PARMevs'], 'metplus_config', os.environ['COMPONENT'], os.environ['RUN']+'_'+os.environ['VERIF_CASE'], os.environ['STEP'], conf_file_name) if not os.path.exists(conf_file): print("FATAL ERROR: "+conf_file+" DOES NOT EXIST") sys.exit(1) metplus_cmd = run_metplus+' -c '+machine_conf+' -c '+conf_file return metplus_cmd def python_command(python_script_name, script_arg_list): """! Write out full call to python Args: python_script_name - python script name (string) script_arg_list - list of script agruments (strings) Returns: python_cmd - full call to python (string) """ python_script = os.path.join(os.environ['USHevs'], os.environ['COMPONENT'], python_script_name) if not os.path.exists(python_script): print("FATAL ERROR: "+python_script+" DOES NOT EXIST") sys.exit(1) python_cmd = 'python '+python_script for script_arg in script_arg_list: python_cmd = python_cmd+' '+script_arg return python_cmd def check_file_exists_size(file_name): """! Checks to see if file exists and has size greater than 0 Args: file_name - file path (string) Returns: file_good - boolean - True: file exists,file size >0 - False: file doesn't exist OR file size = 0 """ if os.path.exists(file_name): if os.path.getsize(file_name) > 0: file_good = True else: print("WARNING: "+file_name+" empty, 0 sized") file_good = False else: print("WARNING: "+file_name+" does not exist") file_good = False return file_good def copy_file(source_file, dest_file): """! This copies a file from one location to another Args: source_file - source file path (string) dest_file - destination file path (string) Returns: """ if check_file_exists_size(source_file): print("Copying "+source_file+" to "+dest_file) shutil.copy2(source_file, dest_file) def convert_grib1_grib2(grib1_file, grib2_file): """! Converts GRIB1 data to GRIB2 Args: grib1_file - string of the path to the GRIB1 file to convert (string) grib2_file - string of the path to save the converted GRIB2 file (string) Returns: """ print("Converting GRIB1 file "+grib1_file+" " +"to GRIB2 file "+grib2_file) cnvgrib = os.environ['CNVGRIB'] os.system(cnvgrib+' -g12 '+grib1_file+' ' +grib2_file+' > /dev/null 2>&1') def convert_grib2_grib1(grib2_file, grib1_file): """! Converts GRIB2 data to GRIB1 Args: grib2_file - string of the path to the GRIB2 file to convert grib1_file - string of the path to save the converted GRIB1 file Returns: """ print("Converting GRIB2 file "+grib2_file+" " +"to GRIB1 file "+grib1_file) cnvgrib = os.environ['CNVGRIB'] os.system(cnvgrib+' -g21 '+grib2_file+' ' +grib1_file+' > /dev/null 2>&1') def convert_grib2_grib2(grib2_fileA, grib2_fileB): """! Converts GRIB2 data to GRIB2 Args: grib2_fileA - string of the path to the GRIB2 file to convert grib2_fileB - string of the path to save the converted GRIB2 file Returns: """ print("Converting GRIB2 file "+grib2_fileA+" " +"to GRIB2 file "+grib2_fileB) cnvgrib = os.environ['CNVGRIB'] os.system(cnvgrib+' -g22 '+grib2_fileA+' ' +grib2_fileB+' > /dev/null 2>&1') def get_time_info(date_start, date_end, date_type, init_hr_list, valid_hr_list, fhr_list): """! Creates a list of dictionaries containing information on the valid dates and times, the initialization dates and times, and forecast hour pairings Args: date_start - verification start date (string, format:YYYYmmdd) date_end - verification end_date (string, format:YYYYmmdd) date_type - how to treat date_start and date_end (string, values:VALID or INIT) init_hr_list - list of initialization hours (string) valid_hr_list - list of valid hours (string) fhr_list - list of forecasts hours (string) Returns: time_info - list of dictionaries with the valid, initalization, and forecast hour pairings """ valid_hr_zfill2_list = [hr.zfill(2) for hr in valid_hr_list] init_hr_zfill2_list = [hr.zfill(2) for hr in init_hr_list] if date_type == 'VALID': date_type_hr_list = valid_hr_zfill2_list elif date_type == 'INIT': date_type_hr_list = init_hr_zfill2_list date_type_hr_start = date_type_hr_list[0] date_type_hr_end = date_type_hr_list[-1] if len(date_type_hr_list) > 1: date_type_hr_inc = np.min( np.diff(np.array(date_type_hr_list, dtype=int)) ) else: date_type_hr_inc = 24 date_start_dt = datetime.datetime.strptime(date_start+date_type_hr_start, '%Y%m%d%H') date_end_dt = datetime.datetime.strptime(date_end+date_type_hr_end, '%Y%m%d%H') time_info = [] date_dt = date_start_dt while date_dt <= date_end_dt: if date_type == 'VALID': valid_time_dt = date_dt elif date_type == 'INIT': init_time_dt = date_dt for fhr in fhr_list: if fhr == 'anl': forecast_hour = 0 else: forecast_hour = int(fhr) if date_type == 'VALID': init_time_dt = (valid_time_dt - datetime.timedelta(hours=forecast_hour)) elif date_type == 'INIT': valid_time_dt = (init_time_dt + datetime.timedelta(hours=forecast_hour)) if valid_time_dt.strftime('%H') in valid_hr_zfill2_list \ and init_time_dt.strftime('%H') in init_hr_zfill2_list: t = {} t['valid_time'] = valid_time_dt t['init_time'] = init_time_dt t['forecast_hour'] = str(forecast_hour) time_info.append(t) date_dt = date_dt + datetime.timedelta(hours=int(date_type_hr_inc)) return time_info def format_filler(unfilled_file_format, valid_time_dt, init_time_dt, forecast_hour, str_sub_dict): """! Creates a filled file path from a format Args: unfilled_file_format - file naming convention (string) valid_time_dt - valid time (datetime) init_time_dt - initialization time (datetime) forecast_hour - forecast hour (string) str_sub_dict - other strings to substitue (dictionary) Returns: filled_file_format - file_format filled in with verifying time information (string) """ filled_file_format = '/' format_opt_list = ['lead', 'lead_shift', 'valid', 'valid_shift', 'init', 'init_shift'] if len(list(str_sub_dict.keys())) != 0: format_opt_list = format_opt_list+list(str_sub_dict.keys()) for filled_file_format_chunk in unfilled_file_format.split('/'): for format_opt in format_opt_list: nformat_opt = ( filled_file_format_chunk.count('{'+format_opt+'?fmt=') ) if nformat_opt > 0: format_opt_count = 1 while format_opt_count <= nformat_opt: if format_opt in ['lead_shift', 'valid_shift', 'init_shift']: shift = (filled_file_format_chunk \ .partition('shift=')[2] \ .partition('}')[0]) format_opt_count_fmt = ( filled_file_format_chunk \ .partition('{'+format_opt+'?fmt=')[2] \ .rpartition('?')[0] ) else: format_opt_count_fmt = ( filled_file_format_chunk \ .partition('{'+format_opt+'?fmt=')[2] \ .partition('}')[0] ) if format_opt == 'valid': replace_format_opt_count = valid_time_dt.strftime( format_opt_count_fmt ) elif format_opt == 'lead': if format_opt_count_fmt == '%1H': if int(forecast_hour) < 10: replace_format_opt_count = forecast_hour[1] else: replace_format_opt_count = forecast_hour elif format_opt_count_fmt == '%2H': replace_format_opt_count = forecast_hour.zfill(2) elif format_opt_count_fmt == '%3H': replace_format_opt_count = forecast_hour.zfill(3) else: replace_format_opt_count = forecast_hour elif format_opt == 'init': replace_format_opt_count = init_time_dt.strftime( format_opt_count_fmt ) elif format_opt == 'lead_shift': shift = (filled_file_format_chunk.partition('shift=')[2]\ .partition('}')[0]) forecast_hour_shift = str(int(forecast_hour) + int(shift)) if format_opt_count_fmt == '%1H': if int(forecast_hour_shift) < 10: replace_format_opt_count = ( forecast_hour_shift[1] ) else: replace_format_opt_count = forecast_hour_shift elif format_opt_count_fmt == '%2H': replace_format_opt_count = ( forecast_hour_shift.zfill(2) ) elif format_opt_count_fmt == '%3H': replace_format_opt_count = ( forecast_hour_shift.zfill(3) ) else: replace_format_opt_count = forecast_hour_shift elif format_opt == 'init_shift': shift = (filled_file_format_chunk.partition('shift=')[2]\ .partition('}')[0]) init_shift_time_dt = ( init_time_dt + datetime.timedelta(hours=int(shift)) ) replace_format_opt_count = init_shift_time_dt.strftime( format_opt_count_fmt ) elif format_opt == 'valid_shift': shift = (filled_file_format_chunk.partition('shift=')[2]\ .partition('}')[0]) valid_shift_time_dt = ( valid_time_dt + datetime.timedelta(hours=int(shift)) ) replace_format_opt_count = valid_shift_time_dt.strftime( format_opt_count_fmt ) else: replace_format_opt_count = str_sub_dict[format_opt] if format_opt in ['lead_shift', 'valid_shift', 'init_shift']: filled_file_format_chunk = ( filled_file_format_chunk.replace( '{'+format_opt+'?fmt=' +format_opt_count_fmt +'?shift='+shift+'}', replace_format_opt_count ) ) else: filled_file_format_chunk = ( filled_file_format_chunk.replace( '{'+format_opt+'?fmt=' +format_opt_count_fmt+'}', replace_format_opt_count ) ) format_opt_count+=1 filled_file_format = os.path.join(filled_file_format, filled_file_format_chunk) return filled_file_format def prep_prod_gfs_file(source_file, dest_file, forecast_hour, prep_method): """! Do prep work for GFS production files Args: source_file_format - source file format (string) dest_file - destination file (string) forecast_hour - forecast hour (string) prep_method - name of prep method to do (string) Returns: """ # Environment variables and executables WGRIB2 = os.environ['WGRIB2'] EXECevs = os.environ['EXECevs'] # Working file names prepped_file = os.path.join(os.getcwd(), 'atmos.'+dest_file.rpartition('/')[2]) working_file1 = prepped_file+'.tmp1' # Prep file if prep_method == 'full': if forecast_hour == 0: wgrib_fhr = 'anl' else: wgrib_fhr = forecast_hour thin_var_level_list = [ 'CAPE:surface', 'CAPE:90-0 mb above ground', 'CWAT:entire atmosphere (considered as a single layer)', 'DPT:2 m above ground', 'GUST:surface', 'HGT:1000 mb', 'HGT:925 mb', 'HGT:850 mb', 'HGT:700 mb', 'HGT:500 mb', 'HGT:400 mb', 'HGT:300 mb', 'HGT:250 mb', 'HGT:200 mb', 'HGT:150 mb', 'HGT:100 mb', 'HGT:50 mb','HGT:20 mb', 'HGT:10 mb', 'HGT:5 mb', 'HGT:1 mb', 'HGT:cloud ceiling', 'HGT:tropopause', 'HPBL:surface', 'ICEC:surface', 'ICETK:surface', 'LHTFL:surface', 'O3MR:925 mb', 'O3MR:100 mb', 'O3MR:70 mb', 'O3MR:50 mb', 'O3MR:30 mb', 'O3MR:20 mb', 'O3MR:10 mb', 'O3MR:5 mb', 'O3MR:1 mb', 'PRES:surface', 'PRES:tropopause', 'PRMSL:mean sea level', 'PWAT:entire atmosphere (considered as a single layer)', 'RH:1000 mb', 'RH:925 mb', 'RH:850 mb', 'RH:700 mb', 'RH:500 mb', 'RH:400 mb', 'RH:300 mb', 'RH:250 mb', 'RH:200 mb', 'RH:150 mb', 'RH:100 mb', 'RH:50 mb','RH:20 mb', 'RH:10 mb', 'RH:5 mb', 'RH:1 mb', 'RH:2 m above ground', 'SHTFL:surface', 'SNOD:surface', 'SPFH:1000 mb', 'SPFH:925 mb', 'SPFH:850 mb', 'SPFH:700 mb', 'SPFH:500 mb', 'SPFH:400 mb', 'SPFH:300 mb', 'SPFH:250 mb', 'SPFH:200 mb', 'SPFH:150 mb', 'SPFH:100 mb', 'SPFH:50 mb', 'SPFH:20 mb', 'SPFH:10 mb', 'SPFH:5 mb', 'SPFH:1 mb', 'SPFH:2 m above ground', 'SOILW:0-0.1 m below ground', 'TCDC:entire atmosphere:'+wgrib_fhr, 'TMP:1000 mb', 'TMP:925 mb', 'TMP:850 mb', 'TMP:700 mb', 'TMP:500 mb', 'TMP:400 mb', 'TMP:300 mb', 'TMP:250 mb', 'TMP:200 mb', 'TMP:150 mb', 'TMP:100 mb', 'TMP:50 mb', 'TMP:20 mb', 'TMP:10 mb', 'TMP:5 mb', 'TMP:1 mb', 'TMP:2 m above ground', 'TMP:surface', 'TMP:tropopause', 'TOZNE:entire atmosphere (considered as a single layer)', 'TSOIL:0-0.1 m below ground', 'UGRD:1000 mb', 'UGRD:925 mb', 'UGRD:850 mb', 'UGRD:700 mb', 'UGRD:500 mb', 'UGRD:400 mb', 'UGRD:300 mb', 'UGRD:250 mb', 'UGRD:200 mb', 'UGRD:150 mb', 'UGRD:100 mb', 'UGRD:50 mb', 'UGRD:20 mb', 'UGRD:10 mb', 'UGRD:5 mb', 'UGRD:1 mb', 'UGRD:10 m above ground', 'VGRD:1000 mb', 'VGRD:925 mb', 'VGRD:850 mb', 'VGRD:700 mb', 'VGRD:500 mb', 'VGRD:400 mb', 'VGRD:300 mb', 'VGRD:250 mb', 'VGRD:200 mb', 'VGRD:150 mb', 'VGRD:100 mb', 'VGRD:50 mb', 'VGRD:20 mb', 'VGRD:10 mb', 'VGRD:5 mb', 'VGRD:1 mb', 'VGRD:10 m above ground', 'VIS:surface', 'WEASD:surface' ] # Missing in GFS files: Sea Ice Drift (Velocity) - SICED?? # Sea Ice Extent - Derived from ICEC? # Sea Ice Volume if check_file_exists_size(source_file): run_shell_command(['>', prepped_file]) for thin_var_level in thin_var_level_list: run_shell_command([WGRIB2, '-match', '"'+thin_var_level+'"', source_file+'|'+WGRIB2, '-i', source_file, '-grib', working_file1]) run_shell_command(['cat', working_file1, '>>', prepped_file]) os.remove(working_file1) elif 'precip' in prep_method: if int(forecast_hour) % 24 == 0: thin_var_level = ('APCP:surface:0-' +str(int(int(forecast_hour)/24))) else: thin_var_level = ('APCP:surface:0-'+forecast_hour) if check_file_exists_size(source_file): run_shell_command([WGRIB2, '-match', '"'+thin_var_level+'"', source_file+'|'+WGRIB2, '-i', source_file, '-grib', prepped_file]) copy_file(prepped_file, dest_file) def prep_prod_fnmoc_file(source_file, dest_file, forecast_hour, prep_method): """! Do prep work for FNMOC production files Args: source_file - source file format (string) dest_file - destination file (string) forecast_hour - forecast hour (string) prep_method - name of prep method to do (string) Returns: """ # Environment variables and executables # Working file names prepped_file = os.path.join(os.getcwd(), 'atmos.'+dest_file.rpartition('/')[2]) # Prep file if check_file_exists_size(source_file): convert_grib2_grib2(source_file, prepped_file) copy_file(prepped_file, dest_file) def prep_prod_jma_file(source_file_format, dest_file, forecast_hour, prep_method): """! Do prep work for JMA production files Args: source_file_format - source file format (string) dest_file - destination file (string) forecast_hour - forecast hour (string) prep_method - name of prep method to do (string) Returns: """ # Environment variables and executables WGRIB = os.environ['WGRIB'] EXECevs = os.environ['EXECevs'] JMAMERGE = os.path.join(EXECevs, 'jma_merge') # Working file names prepped_file = os.path.join(os.getcwd(), 'atmos.'+dest_file.rpartition('/')[2]) working_file1 = prepped_file+'.tmp1' working_file2 = prepped_file+'.tmp2' # Prep file if prep_method == 'full': if forecast_hour == 'anl': wgrib_fhr = ':anl' elif int(forecast_hour) == 0: wgrib_fhr = ':anl' else: wgrib_fhr = ':'+forecast_hour+'hr' for hem in ['n', 's']: hem_source_file = source_file_format.replace('{hem?fmt=str}', hem) if hem == 'n': working_file = working_file1 elif hem == 's': working_file = working_file2 if check_file_exists_size(hem_source_file): run_shell_command( [WGRIB+' '+hem_source_file+' | grep "'+wgrib_fhr+'" | ' +WGRIB+' '+hem_source_file+' -i -grib -o ' +working_file] ) if check_file_exists_size(working_file1) \ and check_file_exists_size(working_file2): run_shell_command( [JMAMERGE, working_file1, working_file2, prepped_file] ) elif 'precip' in prep_method: source_file = source_file_format if check_file_exists_size(source_file): run_shell_command( [WGRIB+' '+source_file+' | grep "0-' +forecast_hour+'hr" | '+WGRIB+' '+source_file +' -i -grib -o '+prepped_file] ) copy_file(prepped_file, dest_file) def prep_prod_ecmwf_file(source_file, dest_file, forecast_hour, prep_method): """! Do prep work for ECMWF production files Args: source_file - source file format (string) dest_file - destination file (string) forecast_hour - forecast hour (string) prep_method - name of prep method to do (string) Returns: """ # Environment variables and executables EXECevs = os.environ['EXECevs'] ECMGFSLOOKALIKENEW = os.path.join(EXECevs, 'ecm_gfs_look_alike_new') PCPCONFORM = os.path.join(EXECevs, 'pcpconform') WGRIB = os.environ['WGRIB'] # Working file names prepped_file = os.path.join(os.getcwd(), 'atmos.'+dest_file.rpartition('/')[2]) working_file1 = prepped_file+'.tmp1' # Prep file if prep_method == 'full': if forecast_hour == 'anl': wgrib_fhr = ':anl' elif int(forecast_hour) == 0: wgrib_fhr = ':anl' else: wgrib_fhr = ':'+forecast_hour+'hr' if check_file_exists_size(source_file): run_shell_command( [WGRIB+' '+source_file+' | grep "'+wgrib_fhr+'" | ' +WGRIB+' '+source_file+' -i -grib -o ' +working_file1] ) if check_file_exists_size(working_file1): run_shell_command(['chmod', '750', working_file1]) run_shell_command(['chgrp', 'rstprod', working_file1]) run_shell_command( [ECMGFSLOOKALIKENEW, working_file1, prepped_file] ) elif 'precip' in prep_method: if check_file_exists_size(source_file): run_shell_command( [PCPCONFORM, 'ecmwf', source_file, prepped_file] ) if os.path.exists(prepped_file): run_shell_command(['chmod', '750', prepped_file]) run_shell_command(['chgrp', 'rstprod', prepped_file]) copy_file(prepped_file, dest_file) def prep_prod_ukmet_file(source_file_format, dest_file, forecast_hour, prep_method): """! Do prep work for UKMET production files Args: source_file_format - source file format (string) dest_file - destination file (string) forecast_hour - forecast hour (string) prep_method - name of prep method to do (string) Returns: """ # Environment variables and executables EXECevs = os.environ['EXECevs'] WGRIB = os.environ['WGRIB'] WGRIB2 = os.environ['WGRIB2'] UKMHIRESMERGE = os.path.join(EXECevs, 'ukm_hires_merge') # Working file names prepped_file = os.path.join(os.getcwd(), 'atmos.'+dest_file.rpartition('/')[2]) working_file1 = prepped_file+'.tmp1' working_file2 = prepped_file+'.tmp2' # Prep file if prep_method == 'full': ukmet_fhr_id_dict = { 'anl': 'AAT', '0': 'AAT', '6': 'BBT', '12': 'CCT', '18': 'DDT', '24': 'EET', '30': 'FFT', '36': 'GGT', '42': 'HHT', '48': 'IIT', '54': 'JJT', '60': 'JJT', '66': 'KKT', '72': 'KKT', '78': 'QQT', '84': 'LLT', '90': 'TTT', '96': 'MMT', '102': 'UUT', '108': 'NNT', '114': 'VVT', '120': 'OOT', '126': '11T', '132': 'PPA', '138': '22T', '144': 'PPA' } if forecast_hour in list(ukmet_fhr_id_dict.keys()): if forecast_hour == 'anl': fhr_id = ukmet_fhr_id_dict['anl'] fhr_str = '0' wgrib_fhr = 'anl' else: fhr_id = ukmet_fhr_id_dict[forecast_hour] fhr_str = forecast_hour if forecast_hour == '0': wgrib_fhr = 'anl' else: wgrib_fhr = forecast_hour+'hr' source_file = source_file_format.replace('{letter?fmt=str}', fhr_id) if check_file_exists_size(source_file): run_shell_command( [WGRIB+' '+source_file+' | grep "'+wgrib_fhr +'" | '+WGRIB+' '+source_file+' -i -grib -o ' +working_file1] ) if check_file_exists_size(working_file1): run_shell_command([UKMHIRESMERGE, working_file1, prepped_file, fhr_str]) elif 'precip' in prep_method: source_file = source_file_format source_file_accum = 12 if check_file_exists_size(source_file): run_shell_command( [WGRIB2+' '+source_file+' -if ":TWATP:" -set_var "APCP" ' +'-fi -grib '+working_file1] ) if check_file_exists_size(working_file1): convert_grib2_grib1(working_file1, working_file2) if check_file_exists_size(working_file2): source_file_accum_fhr_start = ( int(forecast_hour) - source_file_accum ) run_shell_command( [WGRIB+' '+working_file2+' | grep "' +str(source_file_accum_fhr_start)+'-' +forecast_hour+'hr" | '+WGRIB+' '+working_file2 +' -i -grib -o '+prepped_file] ) copy_file(prepped_file, dest_file) def prep_prod_dwd_file(source_file, dest_file, forecast_hour, prep_method): """! Do prep work for DWD production files Args: source_file_format - source file format (string) dest_file - destination file (string) forecast_hour - forecast hour (string) prep_method - name of prep method to do (string) Returns: """ # Environment variables and executables EXECevs = os.environ['EXECevs'] PCPCONFORM = os.path.join(EXECevs, 'pcpconform') # Working file names prepped_file = os.path.join(os.getcwd(), 'atmos.'+dest_file.rpartition('/')[2]) #working_file1 = prepped_file+'.tmp1' #### For DWD to run through pcpconform, file name must be #### dwd_YYYYMMDDHH_(hhh)_(hhh).tmp working_file1 = os.path.join(os.getcwd(), source_file.rpartition('/')[2]+'.tmp') # Prep file if 'precip' in prep_method: if check_file_exists_size(source_file): convert_grib2_grib1(source_file, working_file1) if check_file_exists_size(working_file1): run_shell_command( [PCPCONFORM, 'dwd', working_file1, prepped_file] ) copy_file(prepped_file, dest_file) def prep_prod_metfra_file(source_file, dest_file, forecast_hour, prep_method): """! Do prep work for METRFRA production files Args: source_file - source file(string) dest_file - destination file (string) forecast_hour - forecast hour (string) prep_method - name of prep method to do (string) Returns: """ # Environment variables and executables EXECevs = os.environ['EXECevs'] WGRIB = os.environ['WGRIB'] # Temporary file names prepped_file = os.path.join(os.getcwd(), 'atmos.'+dest_file.rpartition('/')[2]) # Prep file if 'precip' in prep_method: file_accum = 24 fhr_accum_start = int(forecast_hour)-file_accum if check_file_exists_size(source_file): run_shell_command( [WGRIB+' '+source_file+' | grep "' +str(fhr_accum_start)+'-' +forecast_hour+'hr" | '+WGRIB+' '+source_file +' -i -grib -o '+prepped_file] ) copy_file(prepped_file, dest_file) def prep_prod_osi_saf_file(daily_source_file_format, daily_dest_file, weekly_source_file_list, weekly_dest_file, weekly_dates): """! Do prep work for OSI-SAF production files Args: daily_source_file_format - daily source file format (string) daily_dest_file - daily destination file (string) weekly_source_file_list - list of daily files to make up weekly average file weekly_dest_file - weekly destination file (string) weekly_dates - date span for weekly dates (tuple of datetimes) Returns: """ # Environment variables and executables FIXevs = os.environ['FIXevs'] CDO_ROOT = os.environ['CDO_ROOT'] # Temporary file names daily_prepped_file = os.path.join(os.getcwd(), 'atmos.' +daily_dest_file.rpartition('/')[2]) weekly_prepped_file = os.path.join(os.getcwd(), 'atmos.' +weekly_dest_file.rpartition('/')[2]) # Prep daily file for hem in ['nh', 'sh']: hem_source_file = daily_source_file_format.replace('{hem?fmt=str}', hem) hem_dest_file = daily_dest_file.replace('multi.', 'multi.'+hem+'.') hem_prepped_file = os.path.join(os.getcwd(), 'atmos.' +hem_dest_file.rpartition('/')[2]) if check_file_exists_size(hem_source_file): run_shell_command( [os.path.join(CDO_ROOT, 'bin', 'cdo'), 'remapbil,' +os.path.join(FIXevs, 'cdo_grids', 'G004.grid'), hem_source_file, hem_prepped_file] ) if hem == 'nh': nh_prepped_file = hem_prepped_file elif hem == 'sh': sh_prepped_file = hem_prepped_file if check_file_exists_size(nh_prepped_file) \ and check_file_exists_size(sh_prepped_file): nh_data = netcdf.Dataset(nh_prepped_file) sh_data = netcdf.Dataset(sh_prepped_file) merged_data = netcdf.Dataset(daily_prepped_file, 'w', format='NETCDF3_CLASSIC') for attr in nh_data.ncattrs(): if attr == 'history': merged_data.setncattr( attr, nh_data.getncattr(attr)+' ' +sh_data.getncattr(attr) ) elif attr == 'southernmost_latitude': merged_data.setncattr(attr, '-90') elif attr == 'area': merged_data.setncattr(attr, 'Global') else: merged_data.setncattr(attr, nh_data.getncattr(attr)) for dim in list(nh_data.dimensions.keys()): merged_data.createDimension(dim, len(nh_data.dimensions[dim])) for var in ['time', 'time_bnds', 'lat', 'lon']: merged_var = merged_data.createVariable( var, nh_data.variables[var].datatype, nh_data.variables[var].dimensions ) for k in nh_data.variables[var].ncattrs(): merged_var.setncatts( {k: nh_data.variables[var].getncattr(k)} ) if var == 'time': merged_var[:] = nh_data.variables[var][:] + 43200 else: merged_var[:] = nh_data.variables[var][:] for var in ['ice_conc', 'status_flag', 'total_uncertainty', 'smearing_uncertainty', 'algorithm_uncertainty']: merged_var = merged_data.createVariable( var, nh_data.variables[var].datatype, ('lat', 'lon') ) for k in nh_data.variables[var].ncattrs(): if k == 'long_name': merged_var.setncatts( {k: nh_data.variables[var].getncattr(k)\ .replace('northern hemisphere', 'globe')} ) else: merged_var.setncatts( {k: nh_data.variables[var].getncattr(k)} ) merged_var_vals = np.ma.masked_equal( np.vstack((sh_data.variables[var][0,:180,:], nh_data.variables[var][0,180:,:])) ,nh_data.variables[var]._FillValue) merged_var[:] = merged_var_vals copy_file(daily_prepped_file, daily_dest_file) def get_model_file(valid_time_dt, init_time_dt, forecast_hour, source_file_format, dest_file_format): """! This get a model file and saves it in the specificed destination Args: valid_time_dt - valid time (datetime) init_time_dt - initialization time (datetime) forecast_hour - forecast hour (string) source_file_format - source file format (string) dest_file_format - destination file format (string) Returns: """ dest_file = format_filler(dest_file_format, valid_time_dt, init_time_dt, forecast_hour, {}) if not os.path.exists(dest_file): source_file = format_filler(source_file_format, valid_time_dt, init_time_dt, forecast_hour, {}) if 'dcom/navgem' in source_file: prep_prod_fnmoc_file(source_file, dest_file, forecast_hour, 'full') elif 'wgrbbul/jma_' in source_file: prep_prod_jma_file(source_file, dest_file, forecast_hour, 'full') elif 'wgrbbul/ecmwf' in source_file: prep_prod_ecmwf_file(source_file, dest_file, forecast_hour, 'full') elif 'wgrbbul/ukmet_hires' in source_file: prep_prod_ukmet_file(source_file, dest_file, forecast_hour, 'full') elif 'qpf_verif/jma' in source_file: prep_prod_jma_file(source_file, dest_file, forecast_hour, 'precip') elif 'qpf_verif/UWD' in source_file: prep_prod_ecmwf_file(source_file, dest_file, forecast_hour, 'precip') elif 'qpf_verif/ukmo' in source_file: prep_prod_ukmet_file(source_file, dest_file, forecast_hour, 'precip') elif 'qpf_verif/dwd' in source_file: prep_prod_dwd_file(source_file, dest_file, forecast_hour, 'precip') elif 'qpf_verif/METFRA' in source_file: prep_prod_metfra_file(source_file, dest_file, forecast_hour, 'precip') else: if os.path.exists(source_file): print("Linking "+source_file+" to "+dest_file) os.symlink(source_file, dest_file) else: print("WARNING: "+source_file+" DOES NOT EXIST") def get_truth_file(valid_time_dt, source_file_format, dest_file_format): """! This get a model file and saves it in the specificed destination Args: valid_time_dt - valid time (datetime) source_file_format - source file format (string) dest_file_format - destination file format (string) Returns: """ dest_file = format_filler(dest_file_format, valid_time_dt, valid_time_dt, ['anl'], {}) if not os.path.exists(dest_file): source_file = format_filler(source_file_format, valid_time_dt, valid_time_dt, ['anl'], {}) if os.path.exists(source_file): print("Linking "+source_file+" to "+dest_file) os.symlink(source_file, dest_file) else: print("WARNING: "+source_file+" DOES NOT EXIST") def get_obs_valid_hrs(obs): """! This returns the valid hour start, end, and increment information for a given observation Args: obs - observation name (string) Returns: valid_hr_start - starting valid hour (integer) valid_hr_end - ending valid hour (integer) valid_hr_inc - valid hour increment (integer) """ obs_valid_hr_dict = { '24hrCCPA': {'valid_hr_start': 12, 'valid_hr_end': 12, 'valid_hr_inc': 24}, '24hrNOHRSC': {'valid_hr_start': 12, 'valid_hr_end': 12, 'valid_hr_inc': 24}, 'OSI-SAF': {'valid_hr_start': 00, 'valid_hr_end': 00, 'valid_hr_inc': 24}, } if obs in list(obs_valid_hr_dict.keys()): valid_hr_start = obs_valid_hr_dict[obs]['valid_hr_start'] valid_hr_end = obs_valid_hr_dict[obs]['valid_hr_end'] valid_hr_inc = obs_valid_hr_dict[obs]['valid_hr_inc'] else: print(f"FATAL ERROR: Cannot get {obs} valid hour information") sys.exit(1) return valid_hr_start, valid_hr_end, valid_hr_inc def get_off_machine_data(job_file, job_name, job_output, machine, user, queue, account): """! This submits a job to the transfer queue to get data that does not reside on current machine Args: job_file - path to job submission file (string) job_name - job submission name (string) job_output - path to write job output (string) machine - machine name (string) user - user name (string) queue - submission queue name (string) account - submission account name (string) Returns: """ # Set up job wall time information walltime = '60' walltime_seconds = ( datetime.timedelta(minutes=int(walltime)).total_seconds() ) walltime = (datetime.datetime.min + datetime.timedelta(minutes=int(walltime))).time() # Submit job print("Submitting "+job_file+" to "+queue) print("Output sent to "+job_output) os.chmod(job_file, 0o755) if machine == 'WCOSS2': os.system('qsub -V -l walltime='+walltime.strftime('%H:%M:%S')+' ' +'-q '+queue+' -A '+account+' -o '+job_output+' ' +'-e '+job_output+' -N '+job_name+' ' +'-l select=1:ncpus=1 '+job_file) job_check_cmd = ('qselect -s QR -u '+user+' '+'-N ' +job_name+' | wc -l') elif machine in ['HERA', 'ORION', 'S4', 'JET']: os.system('sbatch --ntasks=1 --time=' +walltime.strftime('%H:%M:%S')+' --partition='+queue+' ' +'--account='+account+' --output='+job_output+' ' +'--job-name='+job_name+' '+job_file) job_check_cmd = ('squeue -u '+user+' -n '+job_name+' ' +'-t R,PD -h | wc -l') sleep_counter, sleep_checker = 1, 10 while (sleep_counter*sleep_checker) <= walltime_seconds: sleep(sleep_checker) print("Walltime checker: "+str(sleep_counter*sleep_checker)+" " +"out of "+str(int(walltime_seconds))+" seconds") check_job = subprocess.check_output(job_check_cmd, shell=True, encoding='UTF-8') if check_job[0] == '0': break sleep_counter+=1 def initalize_job_env_dict(verif_type, group, verif_case_step_abbrev_type, job): """! This initializes a dictionary of environment variables and their values to be set for the job pulling from environment variables already set previously Args: verif_type - string of the use case name group - string of the group name verif_case_step_abbrev_type - string of reference name in config and environment variables job - string of job name Returns: job_env_dict - dictionary of job settings """ job_env_var_list = [ 'machine', 'evs_ver', 'HOMEevs', 'FIXevs', 'USHevs', 'DATA', 'COMOUT', 'NET', 'RUN', 'VERIF_CASE', 'STEP', 'COMPONENT' ] if group in ['reformat', 'generate', 'gather']: os.environ['MET_TMP_DIR'] = os.path.join( os.environ['DATA'], os.environ['VERIF_CASE']+'_'+os.environ['STEP'], 'METplus_output', 'tmp' ) if not os.path.exists(os.environ['MET_TMP_DIR']): os.makedirs(os.environ['MET_TMP_DIR']) job_env_var_list.extend( ['METPLUS_PATH','log_met_output_to_metplus', 'metplus_verbosity', 'MET_ROOT', 'MET_bin_exec', 'met_verbosity', 'MET_TMP_DIR', 'COMOUT'] ) elif group == 'plot': job_env_var_list.extend(['MET_ROOT', 'met_ver']) job_env_dict = {} for env_var in job_env_var_list: job_env_dict[env_var] = os.environ[env_var] if group == 'plot': job_env_dict['plot_verbosity'] = ( os.environ['metplus_verbosity'] ) job_env_dict['JOB_GROUP'] = group if group in ['reformat', 'generate', 'plot']: job_env_dict['VERIF_TYPE'] = verif_type if group == 'plot': job_env_dict['job_var'] = job else: job_env_dict['job_name'] = job job_env_dict['fhr_start'] = os.environ[ verif_case_step_abbrev_type+'_fhr_min' ] job_env_dict['fhr_end'] = os.environ[ verif_case_step_abbrev_type+'_fhr_max' ] job_env_dict['fhr_inc'] = os.environ[ verif_case_step_abbrev_type+'_fhr_inc' ] if verif_type in ['pres_levs', 'means', 'sfc']: verif_type_valid_hr_list = ( os.environ[verif_case_step_abbrev_type+'_valid_hr_list']\ .split(' ') ) job_env_dict['valid_hr_start'] = ( verif_type_valid_hr_list[0].zfill(2) ) job_env_dict['valid_hr_end'] = ( verif_type_valid_hr_list[-1].zfill(2) ) if len(verif_type_valid_hr_list) > 1: verif_type_valid_hr_inc = np.min( np.diff(np.array(verif_type_valid_hr_list, dtype=int)) ) else: verif_type_valid_hr_inc = 24 job_env_dict['valid_hr_inc'] = str(verif_type_valid_hr_inc) else: if verif_type == 'precip': valid_hr_start, valid_hr_end, valid_hr_inc = ( get_obs_valid_hrs('24hrCCPA') ) elif verif_type == 'snow': valid_hr_start, valid_hr_end, valid_hr_inc = ( get_obs_valid_hrs('24hrNOHRSC') ) elif verif_type == 'sea_ice': valid_hr_start, valid_hr_end, valid_hr_inc = ( get_obs_valid_hrs('OSI-SAF') ) else: valid_hr_start, valid_hr_end, valid_hr_inc = 12, 12, 23 job_env_dict['valid_hr_start'] = str(valid_hr_start) job_env_dict['valid_hr_end'] = str(valid_hr_end) job_env_dict['valid_hr_inc'] = str(valid_hr_inc) verif_type_init_hr_list = ( os.environ[verif_case_step_abbrev_type+'_init_hr_list']\ .split(' ') ) job_env_dict['init_hr_start'] = ( verif_type_init_hr_list[0].zfill(2) ) job_env_dict['init_hr_end'] = ( verif_type_init_hr_list[-1].zfill(2) ) if len(verif_type_init_hr_list) > 1: verif_type_init_hr_inc = np.min( np.diff(np.array(verif_type_init_hr_list, dtype=int)) ) else: verif_type_init_hr_inc = 24 job_env_dict['init_hr_inc'] = str(verif_type_init_hr_inc) return job_env_dict def get_plot_dates(logger, date_type, start_date, end_date, valid_hr_start, valid_hr_end, valid_hr_inc, init_hr_start, init_hr_end, init_hr_inc, forecast_hour): """! This builds the dates to include in plotting based on user configurations Args: logger - logger object date_type - type of date to plot (string: VALID or INIT) start_date - plotting start date (string, format: YYYYmmdd) end_date - plotting end date (string, format: YYYYmmdd) valid_hr_start - starting valid hour (string) valid_hr_end - ending valid hour (string) valid_hr_inc - valid hour increment (string) init_hr_start - starting initialization hour (string) init_hr_end - ending initialization hour (string) init_hr_inc - initialization hour incrrement (string) forecast_hour - forecast hour (string) Returns: valid_dates - array of valid dates (datetime) init_dates - array of initalization dates (datetime) """ # Build date_type date array if date_type == 'VALID': start_date_dt = datetime.datetime.strptime(start_date+valid_hr_start, '%Y%m%d%H') end_date_dt = datetime.datetime.strptime(end_date+valid_hr_end, '%Y%m%d%H') dt_inc = datetime.timedelta(hours=int(valid_hr_inc)) elif date_type == 'INIT': start_date_dt = datetime.datetime.strptime(start_date+init_hr_start, '%Y%m%d%H') end_date_dt = datetime.datetime.strptime(end_date+init_hr_end, '%Y%m%d%H') dt_inc = datetime.timedelta(hours=int(init_hr_inc)) date_type_dates = (np.arange(start_date_dt, end_date_dt+dt_inc, dt_inc)\ .astype(datetime.datetime)) # Build valid and init date arrays if date_type == 'VALID': valid_dates = date_type_dates init_dates = (valid_dates - datetime.timedelta(hours=(int(forecast_hour)))) elif date_type == 'INIT': init_dates = date_type_dates valid_dates = (init_dates + datetime.timedelta(hours=(int(forecast_hour)))) # Check if unrequested hours exist in arrays, and remove valid_remove_idx_list = [] valid_hr_list = [ str(hr).zfill(2) for hr in range(int(valid_hr_start), int(valid_hr_end)+int(valid_hr_inc), int(valid_hr_inc)) ] for d in range(len(valid_dates)): if valid_dates[d].strftime('%H') \ not in valid_hr_list: valid_remove_idx_list.append(d) valid_dates = np.delete(valid_dates, valid_remove_idx_list) init_dates = np.delete(init_dates, valid_remove_idx_list) init_remove_idx_list = [] init_hr_list = [ str(hr).zfill(2) for hr in range(int(init_hr_start), int(init_hr_end)+int(init_hr_inc), int(init_hr_inc)) ] for d in range(len(init_dates)): if init_dates[d].strftime('%H') \ not in init_hr_list: init_remove_idx_list.append(d) valid_dates = np.delete(valid_dates, init_remove_idx_list) init_dates = np.delete(init_dates, init_remove_idx_list) return valid_dates, init_dates def get_met_line_type_cols(logger, met_root, met_version, met_line_type): """! Get the MET columns for a specific line type and MET verison Args: logger - logger object met_root - path to MET (string) met_version - MET version number (string) met_line_type - MET line type (string) Returns: met_version_line_type_col_list - list of MET versoin line type colums (strings) """ if met_version.count('.') == 2: met_minor_version = met_version.rpartition('.')[0] elif met_version.count('.') == 1: met_minor_version = met_version met_minor_version_col_file = os.path.join( met_root, 'share', 'met', 'table_files', 'met_header_columns_V'+met_minor_version+'.txt' ) if os.path.exists(met_minor_version_col_file): with open(met_minor_version_col_file) as f: for line in f: if met_line_type in line: line_type_cols = line.split(' : ')[-1] break else: logger.error(f"FATAL ERROR: {met_minor_version_col_file} DOES NOT EXIST, " +"cannot determine MET data column structure") sys.exit(1) met_version_line_type_col_list = ( line_type_cols.replace('\n', '').split(' ') ) return met_version_line_type_col_list def format_thresh(thresh): """! Format threshold with letter and symbol options Args: thresh - the threshold (string) Return: thresh_symbol - threshold with symbols (string) thresh_letters - treshold with letters (string) """ for opt in ['>=', '>', '==', '!=', '<=', '<', 'ge', 'gt', 'eq', 'ne', 'le', 'lt']: if opt in thresh: thresh_opt = opt thresh_value = thresh.replace(opt, '') if thresh_opt in ['>', 'gt']: thresh_symbol = '>'+thresh_value thresh_letter = 'gt'+thresh_value elif thresh_opt in ['>=', 'ge']: thresh_symbol = '>='+thresh_value thresh_letter = 'ge'+thresh_value elif thresh_opt in ['<', 'lt']: thresh_symbol = '<'+thresh_value thresh_letter = 'lt'+thresh_value elif thresh_opt in ['<=', 'le']: thresh_symbol = '<='+thresh_value thresh_letter = 'le'+thresh_value elif thresh_opt in ['==', 'eq']: thresh_symbol = '=='+thresh_value thresh_letter = 'eq'+thresh_value elif thresh_opt in ['!=', 'ne']: thresh_symbol = '!='+thresh_value thresh_letter = 'ne'+thresh_value return thresh_symbol, thresh_letter def condense_model_stat_files(logger, input_dir, output_file, model, obs, grid, vx_mask, fcst_var_name, obs_var_name, line_type): """! Condense the individual date model stat file and thin out unneeded data Args: logger - logger object input_dir - path to input directory (string) output_file - path to output file (string) model - model name (string) obs - observation name (string) grid - verification grid (string) vx_mask - verification masking region (string) fcst_var_name - forecast variable name (string) obs_var_name - observation variable name (string) line_type - MET line type (string) Returns: """ model_stat_files_wildcard = os.path.join(input_dir, model, model+'_*.stat') model_stat_files = glob.glob(model_stat_files_wildcard, recursive=True) if len(model_stat_files) == 0: logger.warning(f"NO STAT FILES IN MATCHING " +f"{model_stat_files_wildcard}") else: if not os.path.exists(output_file): logger.debug(f"Condensing down stat files matching " +f"{model_stat_files_wildcard}") with open(model_stat_files[0]) as msf: met_header_cols = msf.readline() all_grep_output = '' grep_opts = ( ' | grep "'+obs+' "' +' | grep "'+grid+' "' +' | grep "'+vx_mask+' "' +' | grep "'+fcst_var_name+' "' +' | grep "'+obs_var_name+' "' +' | grep "'+line_type+' "' ) for model_stat_file in model_stat_files: logger.debug(f"Getting data from {model_stat_file}") ps = subprocess.Popen( 'grep -R "'+model+' " '+model_stat_file+grep_opts, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, encoding='UTF-8' ) logger.debug(f"Ran {ps.args}") all_grep_output = all_grep_output+ps.communicate()[0] logger.debug(f"Condensed {model} .stat file at " +f"{output_file}") with open(output_file, 'w') as f: f.write(met_header_cols+all_grep_output) def build_df(logger, input_dir, output_dir, model_info_dict, plot_info_dict, date_info_dict, met_info_dict, dates, met_format_valid_dates): """! Build the data frame for all model stats, Read the model parse file, if doesn't exist parse the model file for need information, and write file Args: logger - logger object input_dir - path to input directory (string) output_dir - path to output directory (string) model_info_dict - model infomation dictionary (strings) plot_info_dict - plot information dictionary (strings) date_info_dict - date information dictionary (strings) met_info_dict - MET information dictionary (strings) dates - array of dates (datetime) met_format_valid_dates - list of valid dates formatted like they are in MET stat files Returns: """ met_version_line_type_col_list = get_met_line_type_cols( logger, met_info_dict['root'], met_info_dict['version'], plot_info_dict['line_type'] ) df_dtype_dict = {} float_idx = met_version_line_type_col_list.index('TOTAL') for col in met_version_line_type_col_list: col_idx = met_version_line_type_col_list.index(col) if col_idx < float_idx: df_dtype_dict[col] = str else: df_dtype_dict[col] = np.float64 for model_num in list(model_info_dict.keys()): model_num_name = ( model_num+'/'+model_info_dict[model_num]['name'] +'/'+model_info_dict[model_num]['plot_name'] ) model_num_df_index = pd.MultiIndex.from_product( [[model_num_name], met_format_valid_dates], names=['model', 'valid_dates'] ) model_num_df = pd.DataFrame(np.nan, index=model_num_df_index, columns=met_version_line_type_col_list) model_dict = model_info_dict[model_num] parsed_model_stat_file = os.path.join( output_dir, 'fcst'+model_dict['name']+'_' +plot_info_dict['fcst_var_name'] +plot_info_dict['fcst_var_level'] +plot_info_dict['fcst_var_thresh']+'_' +'obs'+model_dict['obs_name']+'_' +plot_info_dict['obs_var_name'] +plot_info_dict['obs_var_level'] +plot_info_dict['obs_var_thresh']+'_' +'linetype'+plot_info_dict['line_type']+'_' +'grid'+plot_info_dict['grid']+'_' +'vxmask'+plot_info_dict['vx_mask']+'_' +'interp'+plot_info_dict['interp_method'] +plot_info_dict['interp_points']+'_' +date_info_dict['date_type'].lower() +dates[0].strftime('%Y%m%d%H%M%S')+'to' +dates[-1].strftime('%Y%m%d%H%M%S')+'_' +'fhr'+date_info_dict['forecast_hour'].zfill(3) +'.stat' ) if not os.path.exists(parsed_model_stat_file): condensed_model_file = os.path.join( input_dir, model_num+'_'+model_dict['name']+'.stat' ) if not os.path.exists(condensed_model_file): condense_model_stat_files( logger, input_dir, condensed_model_file, model_dict['name'], model_dict['obs_name'], plot_info_dict['grid'], plot_info_dict['vx_mask'], plot_info_dict['fcst_var_name'], plot_info_dict['obs_var_name'], plot_info_dict['line_type'] ) if plot_info_dict['fcst_var_thresh'] != 'NA': fcst_var_thresh_symbol, fcst_vat_thresh_letter = ( format_thresh(plot_info_dict['fcst_var_thresh']) ) else: fcst_var_thresh_symbol = plot_info_dict['fcst_var_thresh'] fcst_vat_thresh_letter = plot_info_dict['fcst_var_thresh'] if plot_info_dict['obs_var_thresh'] != 'NA': obs_var_thresh_symbol, obs_vat_thresh_letter = ( format_thresh(plot_info_dict['obs_var_thresh']) ) else: obs_var_thresh_symbol = plot_info_dict['obs_var_thresh'] obs_vat_thresh_letter = plot_info_dict['obs_var_thresh'] if os.path.exists(condensed_model_file): logger.debug(f"Parsing file {condensed_model_file}") condensed_model_df = pd.read_csv( condensed_model_file, sep=" ", skiprows=1, skipinitialspace=True, names=met_version_line_type_col_list, keep_default_na=False, dtype='str', header=None ) parsed_model_df = condensed_model_df[ (condensed_model_df['MODEL'] == model_dict['name']) & (condensed_model_df['DESC'] == plot_info_dict['grid']) & (condensed_model_df['FCST_LEAD'] \ == date_info_dict['forecast_hour'].zfill(2)+'0000') & (condensed_model_df['FCST_VAR'] \ == plot_info_dict['fcst_var_name']) & (condensed_model_df['FCST_LEV'] \ == plot_info_dict['fcst_var_level']) & (condensed_model_df['OBS_VAR'] \ == plot_info_dict['obs_var_name']) & (condensed_model_df['OBS_LEV'] \ == plot_info_dict['obs_var_level']) & (condensed_model_df['OBTYPE'] == model_dict['obs_name']) & (condensed_model_df['VX_MASK'] \ == plot_info_dict['vx_mask']) & (condensed_model_df['INTERP_MTHD'] \ == plot_info_dict['interp_method']) & (condensed_model_df['INTERP_PNTS'] \ == plot_info_dict['interp_points']) & (condensed_model_df['FCST_THRESH'] \ == fcst_var_thresh_symbol) & (condensed_model_df['OBS_THRESH'] \ == obs_var_thresh_symbol) & (condensed_model_df['LINE_TYPE'] \ == plot_info_dict['line_type']) ] parsed_model_df = parsed_model_df[ parsed_model_df['FCST_VALID_BEG'].isin(met_format_valid_dates) ] parsed_model_df['FCST_VALID_BEG'] = pd.to_datetime( parsed_model_df['FCST_VALID_BEG'], format='%Y%m%d_%H%M%S' ) parsed_model_df = parsed_model_df.sort_values(by='FCST_VALID_BEG') parsed_model_df['FCST_VALID_BEG'] = ( parsed_model_df['FCST_VALID_BEG'].dt.strftime('%Y%m%d_%H%M%S') ) parsed_model_df.to_csv( parsed_model_stat_file, header=met_version_line_type_col_list, index=None, sep=' ', mode='w' ) if os.path.exists(parsed_model_stat_file): logger.debug(f"Parsed {model_dict['name']} file " +f"at {parsed_model_stat_file}") else: logger.debug(f"Could not create {parsed_model_stat_file}") if os.path.exists(parsed_model_stat_file): logger.debug(f"Reading {parsed_model_stat_file} for " +f"{model_dict['name']}") model_stat_file_df = pd.read_csv( parsed_model_stat_file, sep=" ", skiprows=1, skipinitialspace=True, names=met_version_line_type_col_list, na_values=['NA'], header=None ) model_stat_file_df = model_stat_file_df.astype(df_dtype_dict) for valid_date in met_format_valid_dates: model_stat_file_df_valid_date_idx_list = ( model_stat_file_df.index[ model_stat_file_df['FCST_VALID_BEG'] == valid_date ] ).tolist() if len(model_stat_file_df_valid_date_idx_list) == 0: logger.debug(f"No data matching valid date {valid_date} " +f"in {parsed_model_stat_file}") continue elif len(model_stat_file_df_valid_date_idx_list) > 1: logger.debug(f"Multiple lines matching valid date " +f"{valid_date} in {parsed_model_stat_file}" +f"using first one") else: logger.debug(f"One line matching valid date " +f"{valid_date} in {parsed_model_stat_file}") model_num_df.loc[(model_num_name, valid_date)] = ( model_stat_file_df.loc\ [model_stat_file_df_valid_date_idx_list[0]]\ [:] ) else: logger.warning(f"{parsed_model_stat_file} does not exist") if model_num == 'model1': all_model_df = model_num_df else: all_model_df = pd.concat([all_model_df, model_num_df]) return all_model_df def calculate_stat(logger, data_df, line_type, stat): """! Calculate the statistic from the data from the read in MET .stat file(s) Args: data_df - dataframe containing the model(s) information from the MET .stat files line_type - MET line type (string) stat - statistic to calculate (string) Returns: stat_df - dataframe of the statistic stat_array - array of the statistic """ if line_type == 'SL1L2': FBAR = data_df.loc[:]['FBAR'] OBAR = data_df.loc[:]['OBAR'] FOBAR = data_df.loc[:]['FOBAR'] FFBAR = data_df.loc[:]['FFBAR'] OOBAR = data_df.loc[:]['OOBAR'] elif line_type == 'SAL1L2': FABAR = data_df.loc[:]['FABAR'] OABAR = data_df.loc[:]['OABAR'] FOABAR = data_df.loc[:]['FOABAR'] FFABAR = data_df.loc[:]['FFABAR'] OOABAR = data_df.loc[:]['OOABAR'] elif line_type == 'CNT': FBAR = data_df.loc[:]['FBAR'] FBAR_NCL = data_df.loc[:]['FBAR_NCL'] FBAR_NCU = data_df.loc[:]['FBAR_NCU'] FBAR_BCL = data_df.loc[:]['FBAR_BCL'] FBAR_BCU = data_df.loc[:]['FBAR_BCU'] FSTDEV = data_df.loc[:]['FSTDEV'] FSTDEV_NCL = data_df.loc[:]['FSTDEV_NCL'] FSTDEV_NCU = data_df.loc[:]['FSTDEV_NCU'] FSTDEV_BCL = data_df.loc[:]['FSTDEV_BCL'] FSTDEV_BCU = data_df.loc[:]['FSTDEV_BCU'] OBAR = data_df.loc[:]['OBAR'] OBAR_NCL = data_df.loc[:]['OBAR_NCL'] OBAR_NCU = data_df.loc[:]['OBAR_NCU'] OBAR_BCL = data_df.loc[:]['OBAR_BCL'] OBAR_BCU = data_df.loc[:]['OBAR_BCU'] OSTDEV = data_df.loc[:]['OSTDEV'] OSTDEV_NCL = data_df.loc[:]['OSTDEV_NCL'] OSTDEV_NCU = data_df.loc[:]['OSTDEV_NCU'] OSTDEV_BCL = data_df.loc[:]['OSTDEV_BCL'] OSTDEV_BCU = data_df.loc[:]['OSTDEV_BCU'] PR_CORR = data_df.loc[:]['PR_CORR'] PR_CORR_NCL = data_df.loc[:]['PR_CORR_NCL'] PR_CORR_NCU = data_df.loc[:]['PR_CORR_NCU'] PR_CORR_BCL = data_df.loc[:]['PR_CORR_BCL'] PR_CORR_BCU = data_df.loc[:]['PR_CORR_BCU'] SP_CORR = data_df.loc[:]['SP_CORR'] KT_CORR = data_df.loc[:]['KT_CORR'] RANKS = data_df.loc[:]['RANKS'] FRANKS_TIES = data_df.loc[:]['FRANKS_TIES'] ORANKS_TIES = data_df.loc[:]['ORANKS_TIES'] ME = data_df.loc[:]['ME'] ME_NCL = data_df.loc[:]['ME_NCL'] ME_NCU = data_df.loc[:]['ME_NCU'] ME_BCL = data_df.loc[:]['ME_BCL'] ME_BCU = data_df.loc[:]['ME_BCU'] ESTDEV = data_df.loc[:]['ESTDEV'] ESTDEV_NCL = data_df.loc[:]['ESTDEV_NCL'] ESTDEV_NCU = data_df.loc[:]['ESTDEV_NCU'] ESTDEV_BCL = data_df.loc[:]['ESTDEV_BCL'] ESTDEV_BCU = data_df.loc[:]['ESTDEV_BCU'] MBIAS = data_df.loc[:]['MBIAS'] MBIAS_BCL = data_df.loc[:]['MBIAS_BCL'] MBIAS_BCU = data_df.loc[:]['MBIAS_BCU'] MAE = data_df.loc[:]['MAE'] MAE_BCL = data_df.loc[:]['MAE_BCL'] MAE_BCU = data_df.loc[:]['MAE_BCU'] MSE = data_df.loc[:]['MSE'] MSE_BCL = data_df.loc[:]['MSE_BCL'] MSE_BCU = data_df.loc[:]['MSE_BCU'] BCRMSE = data_df.loc[:]['BCRMSE'] BCRMSE_BCL = data_df.loc[:]['BCRMSE_BCL'] BCRMSE_BCU = data_df.loc[:]['BCRMSE_BCU'] RMSE = data_df.loc[:]['RMSE'] RMSE_BCL = data_df.loc[:]['RMSE_BCL'] RMSE_BCU = data_df.loc[:]['RMSE_BCU'] E10 = data_df.loc[:]['E10'] E10_BCL = data_df.loc[:]['E10_BCL'] E10_BCU = data_df.loc[:]['E10_BCU'] E25 = data_df.loc[:]['E25'] E25_BCL = data_df.loc[:]['E25_BCL'] E25_BCU = data_df.loc[:]['E25_BCU'] E50 = data_df.loc[:]['E50'] E50_BCL = data_df.loc[:]['E50_BCL'] E50_BCU = data_df.loc[:]['E50_BCU'] E75 = data_df.loc[:]['E75'] E75_BCL = data_df.loc[:]['E75_BCL'] E75_BCU = data_df.loc[:]['E75_BCU'] E90 = data_df.loc[:]['E90'] E90_BCL = data_df.loc[:]['E90_BCL'] E90_BCU = data_df.loc[:]['E90_BCU'] IQR = data_df.loc[:]['IQR'] IQR_BCL = data_df.loc[:]['IQR_BCL'] IQR_BCU = data_df.loc[:]['IQR_BCU'] MAD = data_df.loc[:]['MAD'] MAD_BCL = data_df.loc[:]['MAD_BCL'] MAD_BCU = data_df.loc[:]['MAD_BCU'] ANOM_CORR_NCL = data_df.loc[:]['ANOM_CORR_NCL'] ANOM_CORR_NCU = data_df.loc[:]['ANOM_CORR_NCU'] ANOM_CORR_BCL = data_df.loc[:]['ANOM_CORR_BCL'] ANOM_CORR_BCU = data_df.loc[:]['ANOM_CORR_BCU'] ME2 = data_df.loc[:]['ME2'] ME2_BCL = data_df.loc[:]['ME2_BCL'] ME2_BCU = data_df.loc[:]['ME2_BCU'] MSESS = data_df.loc[:]['MSESS'] MSESS_BCL = data_df.loc[:]['MSESS_BCL'] MSESS_BCU = data_df.loc[:]['MSESS_BCU'] RMSFA = data_df.loc[:]['RMSFA'] RMSFA_BCL = data_df.loc[:]['RMSFA_BCL'] RMSFA_BCU = data_df.loc[:]['RMSFA_BCU'] RMSOA = data_df.loc[:]['RMSOA'] RMSOA_BCL = data_df.loc[:]['RMSOA_BCL'] RMSOA_BCU = data_df.loc[:]['RMSOA_BCU'] ANOM_CORR_UNCNTR = data_df.loc[:]['ANOM_CORR_UNCNTR'] ANOM_CORR_UNCNTR_BCL = data_df.loc[:]['ANOM_CORR_UNCNTR_BCL'] ANOM_CORR_UNCNTR_BCU = data_df.loc[:]['ANOM_CORR_UNCNTR_BCU'] SI = data_df.loc[:]['SI'] SI_BCL = data_df.loc[:]['SI_BCL'] SI_BCU = data_df.loc[:]['SI_BCU'] elif line_type == 'GRAD': FGBAR = data_df.loc[:]['FGBAR'] OGBAR = data_df.loc[:]['OGBAR'] MGBAR = data_df.loc[:]['MGBAR'] EGBAR = data_df.loc[:]['EGBAR'] S1 = data_df.loc[:]['S1'] S1_OG = data_df.loc[:]['S1_OG'] FGOG_RATIO = data_df.loc[:]['FGOG_RATIO'] DX = data_df.loc[:]['DX'] DY = data_df.loc[:]['DY'] elif line_type == 'FHO': F_RATE = data_df.loc[:]['F_RATE'] H_RATE = data_df.loc[:]['H_RATE'] O_RATE = data_df.loc[:]['O_RATE'] elif line_type in ['CTC', 'NBRCTC']: FY_OY = data_df.loc[:]['FY_OY'] FY_ON = data_df.loc[:]['FY_ON'] FN_OY = data_df.loc[:]['FN_OY'] FN_ON = data_df.loc[:]['FN_ON'] elif line_type in ['CTS', 'NBRCTS']: BASER = data_df.loc[:]['BASER'] BASER_NCL = data_df.loc[:]['BASER_NCL'] BASER_NCU = data_df.loc[:]['BASER_NCU'] BASER_BCL = data_df.loc[:]['BASER_BCL'] BASER_BCU = data_df.loc[:]['BASER_BCU'] FMEAN = data_df.loc[:]['FMEAN'] FMEAN_NCL = data_df.loc[:]['FMEAN_NCL'] FMEAN_NCU = data_df.loc[:]['FMEAN_NCU'] FMEAN_BCL = data_df.loc[:]['FMEAN_BCL'] FMEAN_BCU = data_df.loc[:]['FMEAN_BCU'] ACC = data_df.loc[:]['ACC'] ACC_NCL = data_df.loc[:]['ACC_NCL'] ACC_NCU = data_df.loc[:]['ACC_NCU'] ACC_BCL = data_df.loc[:]['ACC_BCL'] ACC_BCU = data_df.loc[:]['ACC_BCU'] FBIAS = data_df.loc[:]['FBIAS'] FBIAS_BCL = data_df.loc[:]['FBIAS_BCL'] FBIAS_BCU = data_df.loc[:]['FBIAS_BCU'] PODY = data_df.loc[:]['PODY'] PODY_NCL = data_df.loc[:]['PODY_NCL'] PODY_NCU = data_df.loc[:]['PODY_NCU'] PODY_BCL = data_df.loc[:]['PODY_BCL'] PODY_BCU = data_df.loc[:]['PODY_BCU'] PODN = data_df.loc[:]['PODN'] PODN_NCL = data_df.loc[:]['PODN_NCL'] PODN_NCU = data_df.loc[:]['PODN_NCU'] PODN_BCL = data_df.loc[:]['PODN_BCL'] PODN_BCU = data_df.loc[:]['PODN_BCU'] POFD = data_df.loc[:]['POFD'] POFD_NCL = data_df.loc[:]['POFD_NCL'] POFD_NCU = data_df.loc[:]['POFD_NCU'] POFD_BCL = data_df.loc[:]['POFD_BCL'] POFD_BCU = data_df.loc[:]['POFD_BCU'] FAR = data_df.loc[:]['FAR'] FAR_NCL = data_df.loc[:]['FAR_NCL'] FAR_NCU = data_df.loc[:]['FAR_NCU'] FAR_BCL = data_df.loc[:]['FAR_BCL'] FAR_BCU = data_df.loc[:]['FAR_BCU'] CSI = data_df.loc[:]['CSI'] CSI_NCL = data_df.loc[:]['CSI_NCL'] CSI_NCU = data_df.loc[:]['CSI_NCU'] CSI_BCL = data_df.loc[:]['CSI_BCL'] CSI_BCU = data_df.loc[:]['CSI_BCU'] GSS = data_df.loc[:]['GSS'] GSS_BCL = data_df.loc[:]['GSS_BCL'] GSS_BCU = data_df.loc[:]['GSS_BCU'] HK = data_df.loc[:]['HK'] HK_NCL = data_df.loc[:]['HK_NCL'] HK_NCU = data_df.loc[:]['HK_NCU'] HK_BCL = data_df.loc[:]['HK_BCL'] HK_BCU = data_df.loc[:]['HK_BCU'] HSS = data_df.loc[:]['HSS'] HSS_BCL = data_df.loc[:]['HSS_BCL'] HSS_BCU = data_df.loc[:]['HSS_BCU'] ODDS = data_df.loc[:]['ODDS'] ODDS_NCL = data_df.loc[:]['ODDS_NCL'] ODDS_NCU = data_df.loc[:]['ODDS_NCU'] ODDS_BCL = data_df.loc[:]['ODDS_BCL'] ODDS_BCU = data_df.loc[:]['ODDS_BCU'] LODDS = data_df.loc[:]['LODDS'] LODDS_NCL = data_df.loc[:]['LODDS_NCL'] LODDS_NCU = data_df.loc[:]['LODDS_NCU'] LODDS_BCL = data_df.loc[:]['LODDS_BCL'] LODDS_BCU = data_df.loc[:]['LODDS_BCU'] ORSS = data_df.loc[:]['ORSS'] ORSS_NCL = data_df.loc[:]['ORSS_NCL'] ORSS_NCU = data_df.loc[:]['ORSS_NCU'] ORSS_BCL = data_df.loc[:]['ORSS_BCL'] ORSS_BCU = data_df.loc[:]['ORSS_BCU'] EDS = data_df.loc[:]['EDS'] EDS_NCL = data_df.loc[:]['EDS_NCL'] EDS_NCU = data_df.loc[:]['EDS_NCU'] EDS_BCL = data_df.loc[:]['EDS_BCL'] EDS_BCU = data_df.loc[:]['EDS_BCU'] SEDS = data_df.loc[:]['SEDS'] SEDS_NCL = data_df.loc[:]['SEDS_NCL'] SEDS_NCU = data_df.loc[:]['SEDS_NCU'] SEDS_BCL = data_df.loc[:]['SEDS_BCL'] SEDS_BCU = data_df.loc[:]['SEDS_BCU'] EDI = data_df.loc[:]['EDI'] EDI_NCL = data_df.loc[:]['EDI_NCL'] EDI_NCU = data_df.loc[:]['EDI_NCU'] EDI_BCL = data_df.loc[:]['EDI_BCL'] EDI_BCU = data_df.loc[:]['EDI_BCU'] SEDI = data_df.loc[:]['SEDI'] SEDI_NCL = data_df.loc[:]['SEDI_NCL'] SEDI_NCU = data_df.loc[:]['SEDI_NCU'] SEDI_BCL = data_df.loc[:]['SEDI_BCL'] SEDI_BCU = data_df.loc[:]['SEDI_BCU'] BAGSS = data_df.loc[:]['BAGSS'] BAGSS_BCL = data_df.loc[:]['BAGSS_BCL'] BAGSS_BCU = data_df.loc[:]['BAGSS_BCU'] elif line_type == 'NBRCNT': FBS = data_df.loc[:]['FBS'] FBS_BCL = data_df.loc[:]['FBS_BCL'] FBS_BCU = data_df.loc[:]['FBS_BCU'] FSS = data_df.loc[:]['FSS'] FSS_BCL = data_df.loc[:]['FSS_BCL'] FSS_BCU = data_df.loc[:]['FSS_BCU'] AFSS = data_df.loc[:]['AFSS'] AFSS_BCL = data_df.loc[:]['AFSS_BCL'] AFSS_BCU = data_df.loc[:]['AFSS_BCU'] UFSS = data_df.loc[:]['UFSS'] UFSS_BCL = data_df.loc[:]['UFSS_BCL'] UFSS_BCU = data_df.loc[:]['UFSS_BCU'] F_RATE = data_df.loc[:]['F_RATE'] F_RATE_BCL = data_df.loc[:]['F_RATE_BCL'] F_RATE_BCU = data_df.loc[:]['F_RATE_BCU'] O_RATE = data_df.loc[:]['O_RATE'] O_RATE_BCL = data_df.loc[:]['O_RATE_BCL'] O_RATE_BCU = data_df.loc[:]['O_RATE_BCU'] elif line_type == 'VL1L2': UFBAR = data_df.loc[:]['UFBAR'] VFBAR = data_df.loc[:]['VFBAR'] UOBAR = data_df.loc[:]['UOBAR'] VOBAR = data_df.loc[:]['VOBAR'] UVFOBAR = data_df.loc[:]['UVFOBAR'] UVFFBAR = data_df.loc[:]['UVFFBAR'] UVOOBAR = data_df.loc[:]['UVOOBAR'] elif line_type == 'VAL1L2': UFABAR = data_df.loc[:]['UFABAR'] VFABAR = data_df.loc[:]['VFABAR'] UOABAR = data_df.loc[:]['UOABAR'] VOABAR = data_df.loc[:]['VOABAR'] UVFOABAR = data_df.loc[:]['UVFOABAR'] UVFFABAR = data_df.loc[:]['UVFFABAR'] UVOOABAR = data_df.loc[:]['UVOOABAR'] elif line_type == 'VCNT': FBAR = data_df.loc[:]['FBAR'] OBAR = data_df.loc[:]['OBAR'] FS_RMS = data_df.loc[:]['FS_RMS'] OS_RMS = data_df.loc[:]['OS_RMS'] MSVE = data_df.loc[:]['MSVE'] RMSVE = data_df.loc[:]['RMSVE'] FSTDEV = data_df.loc[:]['FSTDEV'] OSTDEV = data_df.loc[:]['OSTDEV'] FDIR = data_df.loc[:]['FDIR'] ORDIR = data_df.loc[:]['ODIR'] FBAR_SPEED = data_df.loc[:]['FBAR_SPEED'] OBAR_SPEED = data_df.loc[:]['OBAR_SPEED'] VDIFF_SPEED = data_df.loc[:]['VDIFF_SPEED'] VDIFF_DIR = data_df.loc[:]['VDIFF_DIR'] SPEED_ERR = data_df.loc[:]['SPEED_ERR'] SPEED_ABSERR = data_df.loc[:]['SPEED_ABSERR'] DIR_ERR = data_df.loc[:]['DIR_ERR'] DIR_ABSERR = data_df.loc[:]['DIR_ABSERR'] if stat == 'ACC': # Anomaly Correlation Coefficient if line_type == 'SAL1L2': stat_df = (FOABAR - FABAR*OABAR) \ /np.sqrt((FFABAR - FABAR*FABAR)* (OOABAR - OABAR*OABAR)) elif line_type == 'CNT': stat_df = ANOM_CORR elif line_type == 'VAL1L2': stat_df = UVFOABAR/np.sqrt(UVFFABAR*UVOOABAR) elif stat == 'BIAS': # Bias if line_type == 'SL1L2': stat_df = FBAR - OBAR elif line_type == 'CNT': stat_df = ME elif line_type == 'VL1L2': stat_df = np.sqrt(UVFFBAR) - np.sqrt(UVOOBAR) elif stat == 'CSI': # Critical Success Index' if line_type == 'CTC': stat_df = FY_OY/(FY_OY + FY_ON + FN_OY) elif stat in ['ETS', 'GSS']: # Equitable Threat Score/Gilbert Skill Score if line_type == 'CTC': TOTAL = FY_OY + FY_ON + FN_OY + FN_ON C = ((FY_OY + FY_ON)*(FY_OY + FN_OY))/TOTAL stat_df = (FY_OY - C)/(FY_OY + FY_ON + FN_OY - C) elif line_type == 'CTS': stat_df = GSS elif stat == 'FBAR': # Forecast Mean if line_type == 'SL1L2': stat_df = FBAR elif stat == 'FBIAS': # Frequency Bias if line_type == 'CTC': stat_df = (FY_OY + FY_ON)/(FY_OY + FN_OY) elif line_type == 'CTS': stat_df = FBIAS elif stat == 'FSS': # Fraction Skill Score if line_type == 'NBRCNT': stat_df = FSS elif stat == 'FY_OY': # Forecast Yes/Obs Yes if line_type == 'CTC': stat_df = FY_OY elif stat == 'HSS': # Heidke Skill Score if line_type == 'CTC': TOTAL = FY_OY + FY_ON + FN_OY + FN_ON CA = (FY_OY+FY_ON)*(FY_OY+FN_OY) CB = (FN_OY+FN_ON)*(FY_ON+FN_ON) C = (CA + CB)/TOTAL stat_df = (FY_OY + FN_ON - C)/(TOTAL - C) elif stat == 'POD': # Probability of Detection if line_type == 'CTC': stat_df = FY_OY/(FY_OY + FN_OY) elif stat == 'RMSE': # Root Mean Square Error if line_type == 'SL1L2': stat_df = np.sqrt(FFBAR + OOBAR - 2*FOBAR) elif line_type == 'CNT': stat_df = RMSE elif line_type == 'VL1L2': stat_df = np.sqrt(UVFFBAR + UVOOBAR - 2*UVFOBAR) elif stat == 'S1': # S1 if line_type == 'GRAD': stat_df = S1 elif stat == 'SRATIO': # Success Ratio if line_type == 'CTC': stat_df = 1 - (FY_ON/(FY_ON + FY_OY)) else: logger.error("FATAL ERROR: "+stat+" IS NOT AN OPTION") sys.exit(1) idx = 0 idx_dict = {} while idx < stat_df.index.nlevels: idx_dict['index'+str(idx)] = len( stat_df.index.get_level_values(idx).unique() ) idx+=1 if stat_df.index.nlevels == 2: stat_array = stat_df.values.reshape( idx_dict['index0'], idx_dict['index1'] ) return stat_df, stat_array