Note

Go to the end to download the full example code.

PointStat: Hurricane Matthew I-WRF

model_applications/tc_and_extra_tc/PointStat_fcstWRF_obsMADIS_hurricane_matthew.conf

Scientific Objective

Perform verification of WRF data generated over Hurricane Matthew from 2016 using MADIS RAOB and METAR observations. This use case was created originally to serve the needs of the NSF-funded I-WRF project (https://i-wrf.org). The goal of I-WRF is to create a series of connected software containers to enable multi-node WRF simulations, verification with METplus, and visualization that can run on a range of platforms (HPC, cloud environments, and laptops), in order to lower the bar for researchers and students to tackle interesting scientific questions with WRF and METplus. This Hurricane Matthew use case borrows the inputs and configuration from the first simulation in the WRF Online Tutorial, and both the WRF and METplus configurations here were designed to provide a simple, initial demo case for a small domain that would run quickly and confirm for users that the containers are working as expected. The I-WRF project provided the funding to add the capability of METplus to read in native WRF output files directly, in order to avoid the extra complicating step of running WRF output through UPP to generate grib files.

Version Added

METplus version 6.1

Datasets

Forecast: WRF, temperature and wind

Observation: MADIS RAOB and METAR

Climatology: None

Location: All of the input data required for this use case can be found in a sample data tarball. Each use case category will have one or more sample data tarballs. It is only necessary to download the tarball with the use case’s dataset and not the entire collection of sample data. Click here to access the METplus releases page and download sample data for the appropriate release: https://github.com/dtcenter/METplus/releases This tarball should be unpacked into the directory that you will set the value of INPUT_BASE. See Running METplus section for more information.

METplus Components

This use case calls 2 instances of MADIS2NC to convert the RAOB and METAR observations to NetCDF. Then 2 instances of PointStat are called to process surface and upper air fields. Next, UserScript is used to call a python script to plot the WRF data. Finally, another UserScript instance calls a METdataio script to reformat the PointStat output and a METplotpy script to produce a line plot.

METplus Workflow

Beginning time (INIT_BEG): 2016-10-06 00Z

End time (INIT_END): 2016-10-06 00Z

Increment between beginning and end times (INIT_INCREMENT): 6 hours

Sequence of forecast leads to process (LEAD_SEQ): 0, 3, 6, 9, 12, 15, 18, 21, 24, 27, 30, 33, 36, 39, 42, 45, 48

The MADIS2NC instances process every forecast lead from 0-48. The RAOB MADIS2NC instance has missing data for a few of the leads, so a missing data threshold of 87.5% is set to prevent errors. The surface instance of PointStat allows a +/- 30 minute window for observations and includes files from the previous hour. The upper air instance of PointStat allows +/- 90 minute window for observations and includes files from +/- 2 hours around the valid time.

METplus Configuration

METplus first loads all of the configuration files found in parm/metplus_config, then it loads any configuration files passed to METplus via the command line: parm/use_cases/model_applications/tc_and_extra_tc/PointStat_fcstWRF_obsMADIS_hurricane_matthew.conf

[config]

# For additional information, please see the METplus Users Guide.
# https://metplus.readthedocs.io/en/latest/Users_Guide

###
# Processes to run
# https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#process-list
###

PROCESS_LIST = MADIS2NC(metar), MADIS2NC(raob), PointStat(surface), PointStat(upper_air), UserScript(wrf_plot), UserScript(metplotpy)

###
# Time Info
# LOOP_BY options are INIT, VALID, RETRO, and REALTIME
# If set to INIT or RETRO:
#   INIT_TIME_FMT, INIT_BEG, INIT_END, and INIT_INCREMENT must also be set
# If set to VALID or REALTIME:
#   VALID_TIME_FMT, VALID_BEG, VALID_END, and VALID_INCREMENT must also be set
# LEAD_SEQ is the list of forecast leads to process
# https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#timing-control
###

LOOP_BY = INIT
INIT_TIME_FMT = %Y%m%d%H
INIT_BEG = 2016100600
INIT_END = 2016100600
INIT_INCREMENT = 6H

LEAD_SEQ = begin_end_incr(0,48,3)

[metar]
LEAD_SEQ = begin_end_incr(0,48,1)

[raob]

LEAD_SEQ = begin_end_incr(0,48,1)

MADIS2NC_ALLOW_MISSING_INPUTS = True
MADIS2NC_MISSING_INPUT_THRESH = 0.875

[config]

###
# File I/O
# https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#directory-and-filename-template-info
###

INPUT_DIR = {INPUT_BASE}/model_applications/tc_and_extra_tc/PointStat_fcstWRF_obsMADIS_hurricane_matthew
CONFIG_DIR = {PARM_BASE}/use_cases/model_applications/tc_and_extra_tc/PointStat_fcstWRF_obsMADIS_hurricane_matthew

MADIS2NC_INPUT_DIR = {INPUT_DIR}/obs
MADIS2NC_INPUT_TEMPLATE = {instance}/{valid?fmt=%Y%m%d_%H%M}

MADIS2NC_OUTPUT_DIR = {OUTPUT_BASE}/madis2nc
MADIS2NC_OUTPUT_TEMPLATE = {instance}/met_{valid?fmt=%Y%m%d_%H%M}.nc

MADIS2NC_TYPE = {instance}

FCST_POINT_STAT_INPUT_DIR = {INPUT_DIR}/wrf

[surface]

FCST_POINT_STAT_INPUT_TEMPLATE = {init?fmt=%Y%m%d_%H}/wrfout_d01_{valid?fmt=%Y-%m-%d_%H:%M:%S}

OBS_POINT_STAT_INPUT_DIR = {MADIS2NC_OUTPUT_DIR}/metar
OBS_POINT_STAT_INPUT_TEMPLATE = met_{valid?fmt=%Y%m%d_%H%M}.nc

POINT_STAT_MESSAGE_TYPE = ADPSFC

FCST_VAR1_NAME = T2
FCST_VAR1_LEVELS = "(0,*,*)"
FCST_VAR1_OPTIONS = set_attr_level = "Z2"

OBS_VAR1_NAME = TMP
OBS_VAR1_LEVELS = Z2

FCST_VAR2_NAME = U10
FCST_VAR2_LEVELS = "(0,*,*)"
FCST_VAR2_OPTIONS = set_attr_level = "Z10"

OBS_VAR2_NAME = UGRD
OBS_VAR2_LEVELS = Z10

FCST_VAR3_NAME = V10
FCST_VAR3_LEVELS = "(0,*,*)"
FCST_VAR3_OPTIONS = set_attr_level = "Z10"

OBS_VAR3_NAME = VGRD
OBS_VAR3_LEVELS = Z10

OBS_POINT_STAT_WINDOW_BEGIN = -1799
OBS_POINT_STAT_WINDOW_END = 1800

OBS_POINT_STAT_FILE_WINDOW_BEGIN = -1H

[upper_air]

FCST_POINT_STAT_INPUT_TEMPLATE = {init?fmt=%Y%m%d_%H}/wrfout_plev_d01_{valid?fmt=%Y-%m-%d_%H:%M:%S}

OBS_POINT_STAT_INPUT_DIR = {MADIS2NC_OUTPUT_DIR}/raob
OBS_POINT_STAT_INPUT_TEMPLATE = met_{valid?fmt=%Y%m%d_%H%M}.nc

POINT_STAT_MESSAGE_TYPE = ADPUPA

# WRF pressure levels:
# 92500,85000,70000,50000,40000,30000,25000,20000,15000,10000

UPPER_AIR_FCST_LEVELS = "(0,@92500,*,*)", "(0,@85000,*,*)" , "(0,@70000,*,*)", "(0,@50000,*,*)", "(0,@40000,*,*)", "(0,@30000,*,*)", "(0,@25000,*,*)", "(0,@20000,*,*)", "(0,@15000,*,*)", "(0,@10000,*,*)"
UPPER_AIR_OBS_LEVELS = P925, P850, P700, P500, P400, P300, P250, P200, P150, P100

FCST_VAR1_NAME = T_PL
FCST_VAR1_LEVELS = {UPPER_AIR_FCST_LEVELS}

OBS_VAR1_NAME = TMP
OBS_VAR1_LEVELS = {UPPER_AIR_OBS_LEVELS}

FCST_VAR2_NAME = U_PL
FCST_VAR2_LEVELS = {UPPER_AIR_FCST_LEVELS}

OBS_VAR2_NAME = UGRD
OBS_VAR2_LEVELS = {UPPER_AIR_OBS_LEVELS}

FCST_VAR3_NAME = V_PL
FCST_VAR3_LEVELS = {UPPER_AIR_FCST_LEVELS}

OBS_VAR3_NAME = VGRD
OBS_VAR3_LEVELS = {UPPER_AIR_OBS_LEVELS}

OBS_POINT_STAT_WINDOW_BEGIN = -5399
OBS_POINT_STAT_WINDOW_END = 5400

OBS_POINT_STAT_FILE_WINDOW_BEGIN = -2H
OBS_POINT_STAT_FILE_WINDOW_END = 2H

[config]

POINT_STAT_ONCE_PER_FIELD = False

POINT_STAT_OUTPUT_DIR = {OUTPUT_BASE}/point_stat

###
# PointStat Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#pointstat
###

POINT_STAT_MET_CONFIG_OVERRIDES = wind_thresh = [ >2 ]; wind_logic = INTERSECTION;

POINT_STAT_DUPLICATE_FLAG = UNIQUE
POINT_STAT_OBS_SUMMARY = NEAREST

POINT_STAT_INTERP_TYPE_METHOD = BILIN
POINT_STAT_INTERP_TYPE_WIDTH = 2

POINT_STAT_OUTPUT_FLAG_CNT = BOTH
POINT_STAT_OUTPUT_FLAG_SL1L2 = STAT
POINT_STAT_OUTPUT_FLAG_VL1L2 = STAT
POINT_STAT_OUTPUT_FLAG_VCNT = BOTH

MODEL = WRF

POINT_STAT_DESC = NA
OBTYPE =

POINT_STAT_OUTPUT_PREFIX = {instance}

POINT_STAT_MASK_GRID = FULL


[user_env_vars]

PLOT_INPUT_DIR = {POINT_STAT_OUTPUT_DIR}
PLOT_OUTPUT_DIR = {OUTPUT_BASE}/met_plot
PLOT_LOG_DIR = {LOG_DIR}

[metplotpy]

METDATAIO_BASE = {METPLUS_BASE}/../METdataio
METPLOTPY_BASE = {METPLUS_BASE}/../METplotpy

USER_SCRIPT_RUNTIME_FREQ = RUN_ONCE
USER_SCRIPT_COMMAND = mkdir -p {OUTPUT_BASE}/met_plot; cd {OUTPUT_BASE}/met_plot; python {METDATAIO_BASE}/METreformat/write_stat_ascii.py {CONFIG_DIR}/reformat_cnt.yaml; python {METPLOTPY_BASE}/metplotpy/plots/line/line.py {CONFIG_DIR}/custom_line.yaml


[wrf_plot]

USER_SCRIPT_RUNTIME_FREQ = RUN_ONCE
USER_SCRIPT_COMMAND={CONFIG_DIR}/plot_wrf.py -w {INPUT_DIR}/wrf -o {OUTPUT_BASE}/wrf_plot

MET Configuration

METplus sets environment variables based on user settings in the METplus configuration file. See How METplus controls MET config file settings for more details.

YOU SHOULD NOT SET ANY OF THESE ENVIRONMENT VARIABLES YOURSELF! THEY WILL BE OVERWRITTEN BY METPLUS WHEN IT CALLS THE MET TOOLS!

If there is a setting in the MET configuration file that is currently not supported by METplus you’d like to control, please refer to: Overriding Unsupported MET config file settings

PointStatConfig_wrapped

////////////////////////////////////////////////////////////////////////////////
//
// Point-Stat configuration file.
//
// For additional information, see the MET_BASE/config/README file.
//
////////////////////////////////////////////////////////////////////////////////

//
// Output model name to be written
//
// model =
${METPLUS_MODEL}

//
// Output description to be written
// May be set separately in each "obs.field" entry
//
// desc =
${METPLUS_DESC}

////////////////////////////////////////////////////////////////////////////////

//
// Verification grid
//
// regrid = {
${METPLUS_REGRID_DICT}

////////////////////////////////////////////////////////////////////////////////

//
// May be set separately in each "field" entry
//
censor_thresh = [];
censor_val    = [];
cat_thresh    = [ NA ];
cnt_thresh    = [ NA ];
cnt_logic     = UNION;
wind_thresh   = [ NA ];
wind_logic    = UNION;
eclv_points   = 0.05;
//hss_ec_value =
${METPLUS_HSS_EC_VALUE}
rank_corr_flag = FALSE;

//
// Forecast and observation fields to be verified
//
fcst = {
  ${METPLUS_FCST_FILE_TYPE}
  //field = [
  ${METPLUS_FCST_FIELD}
  ${METPLUS_FCST_CLIMO_MEAN_DICT}
  ${METPLUS_FCST_CLIMO_STDEV_DICT}
}

obs = {
  ${METPLUS_OBS_FILE_TYPE}
  //field = [
  ${METPLUS_OBS_FIELD}
  ${METPLUS_OBS_CLIMO_MEAN_DICT}
  ${METPLUS_OBS_CLIMO_STDEV_DICT}
}
////////////////////////////////////////////////////////////////////////////////

//
// Point observation filtering options
// May be set separately in each "obs.field" entry
//
// message_type =
${METPLUS_MESSAGE_TYPE}
sid_exc        = [];

//obs_quality_inc =
${METPLUS_OBS_QUALITY_INC}

//obs_quality_exc =
${METPLUS_OBS_QUALITY_EXC}

//duplicate_flag =
${METPLUS_DUPLICATE_FLAG}

//obs_summary =
${METPLUS_OBS_SUMMARY}

//obs_perc_value =
${METPLUS_OBS_PERC_VALUE}

//
// Mapping of message type group name to comma-separated list of values.
//
//message_type_group_map =
${METPLUS_MESSAGE_TYPE_GROUP_MAP}

//obtype_as_group_val_flag =
${METPLUS_OBTYPE_AS_GROUP_VAL_FLAG}

////////////////////////////////////////////////////////////////////////////////

//
// Climatology data
//
//climo_mean = {
${METPLUS_CLIMO_MEAN_DICT}


//climo_stdev = {
${METPLUS_CLIMO_STDEV_DICT}

//
// May be set separately in each "obs.field" entry
//
//climo_cdf = {
${METPLUS_CLIMO_CDF_DICT}

////////////////////////////////////////////////////////////////////////////////

//
// Land/Sea mask
// For LANDSF message types, only use forecast grid points where land = TRUE.
// For WATERSF message types, only use forecast grid points where land = FALSE.
// land_mask.flag may be set separately in each "obs.field" entry.
//
//land_mask = {
${METPLUS_LAND_MASK_DICT}

//
// Topography
// For SURFACE message types, only use observations where the topo - station
// elevation difference meets the use_obs_thresh threshold.
// For the observations kept, when interpolating forecast data to the
// observation location, only use forecast grid points where the topo - station
// difference meets the interp_fcst_thresh threshold.
// topo_mask.flag may be set separately in each "obs.field" entry.
//
//topo_mask = {
${METPLUS_TOPO_MASK_DICT}

////////////////////////////////////////////////////////////////////////////////

//
// Point observation time window
//
// obs_window = {
${METPLUS_OBS_WINDOW_DICT}

////////////////////////////////////////////////////////////////////////////////

//
// Verification masking regions
//
//mask = {
${METPLUS_MASK_DICT}

////////////////////////////////////////////////////////////////////////////////

//
// Confidence interval settings
//
ci_alpha  = [ 0.05 ];

boot = {
   interval = PCTILE;
   rep_prop = 1.0;
   n_rep    = 0;
   rng      = "mt19937";
   seed     = "";
}

////////////////////////////////////////////////////////////////////////////////

//
// Interpolation methods
//
//interp = {
${METPLUS_INTERP_DICT}

////////////////////////////////////////////////////////////////////////////////

//
// HiRA verification method
//
//hira = {
${METPLUS_HIRA_DICT}

////////////////////////////////////////////////////////////////////////////////
// Threshold for SEEPS p1 (Probability of being dry)

//seeps_p1_thresh =
${METPLUS_SEEPS_P1_THRESH}

////////////////////////////////////////////////////////////////////////////////

//
// Statistical output types
//
//output_flag = {
${METPLUS_OUTPUT_FLAG_DICT}

////////////////////////////////////////////////////////////////////////////////

//ugrid_dataset =
${METPLUS_UGRID_DATASET}

//ugrid_max_distance_km =
${METPLUS_UGRID_MAX_DISTANCE_KM}

//ugrid_coordinates_file =
${METPLUS_UGRID_COORDINATES_FILE}

////////////////////////////////////////////////////////////////////////////////

//point_weight_flag =
${METPLUS_POINT_WEIGHT_FLAG}

tmp_dir = "${MET_TMP_DIR}";

// output_prefix =
${METPLUS_OUTPUT_PREFIX}
//version        = "V10.0.0";

////////////////////////////////////////////////////////////////////////////////

${METPLUS_MET_CONFIG_OVERRIDES}

Python Embedding

This use case does not use Python Embedding.

User Scripting

This use case calls files from METdataio and METcalcpy to generate a line plot of the point_stat output. It also calls a Python script to perform plotting of WRF data.

parm/use_cases/model_applications/tc_and_extra_tc/PointStat_fcstWRF_obsMADIS_hurricane_matthew/plot_wrf.py

#! /usr/bin/env python3

"""
plot_wrf.py

Written by: Jared A. Lee (jaredlee@ucar.edu)
Written on: 13 May 2024
"""

import sys
import argparse
import pathlib
import datetime as dt
import numpy as np
import pandas as pd
import netCDF4
import wrf
import matplotlib as mpl

# Import functions from a local file
import map_funcs

ADD_BARBS_STRING = '+barbs'


def setup_plot_configuration():
    """Set up default plotting configuration and options."""
    return {
        'plot_type': 'png',
        'plot_subdomain': False,
        'plot_stations': True,
        'plot_terrain': True,
        'plot_t2': True,
        'plot_rh2': True,
        'plot_slp': True,
        'plot_ws10': True,
        'plot_refl': True,
        'plot_rain': True,
        'plot_ws100': True,
        'plot_wind_barbs_sfc': True,
        'plot_wind_barbs_upr': True,
        'water_color': 'lightblue',
        'suptitle_y': 1.00,
        'plot_fontsize': 13,
        'barb_thin': 10,
        'barb_width': 0.5,
        # Domain plotting ranges
        'i_beg': 0, 'i_end': -1,
        'j_beg': 0, 'j_end': -1,
        # Station data
        'text1_lab': ['Miami', 'Jacksonville', 'Charleston'],
        'mark1_lat': np.array([25.7617, 30.3322, 32.7833]),
        'mark1_lon': np.array([-80.1918, -81.6557, -79.9320]),
        'mark1_size': 36,
        'mark1_color': 'black',
        'lat_labels': [16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40],
        'lon_labels': [-62, -64, -66, -68, -70, -72, -74, -76, -78, -80, -82, -84, -86, -88],
    }


def setup_plot_limits():
    """Set up plotting limits and contour intervals for all variables."""
    return {
        'terrain': {'min': 0.0, 'max': 1500.1, 'int': 100.0},
        'slp': {'min': 980.0, 'max': 1020.1, 'int': 2.0},
        't2': {'min': 0.0, 'max': 40.1, 'int': 2.0},
        'rh2': {'min': 0.0, 'max': 100.1, 'int': 5.0},
        'ws10': {'min': 0.0, 'max': 35.0, 'int': 2.5},
        'rain': {'min': 0.0, 'max': 100.1, 'int': 5.0},
    }


def setup_constants_and_formats():
    """Define constants, format strings, and other static values."""
    constants = {
        'c_to_k': 273.15,
        'missing_val': -9999.0,
        'mpl_ms1': r'm $\mathregular{s^{-1}}$',
        'deg_uni': '\u00B0',
    }

    formats = {
        'fmt_yyyymmdd_hh': '%Y%m%d_%H',
        'fmt_yyyymmdd_hhmm': '%Y%m%d_%H%M',
        'fmt_wrf_date': '%Y-%m-%d',
        'fmt_wrf_time': '%H:%M:%S',
        'fmt_time_plot': '%d %b %Y/%H%M UTC',
    }
    formats['fmt_wrf_dt'] = formats['fmt_wrf_date'] + '_' + formats['fmt_wrf_time']
    formats['fmt_time_file'] = formats['fmt_yyyymmdd_hhmm']

    # Define custom colormap for radar reflectivity
    cmap_radar = np.array([
        [200, 200, 200], [4, 233, 231], [1, 159, 244], [3, 0, 244],
        [2, 253, 2], [1, 197, 1], [0, 142, 0],
        [253, 248, 2], [229, 188, 0], [253, 149, 0],
        [253, 0, 0], [212, 0, 0], [188, 0, 0],
        [248, 0, 253], [152, 84, 198], [228, 199, 243]], np.float32) / 255.0
    bounds_radar = np.arange(0., 75.01, 5.0)

    return constants, formats, {'cmap_radar': cmap_radar, 'bounds_radar': bounds_radar}


def adjust_subdomain_ranges(config):
    """Adjust domain ranges if subdomain plotting is requested."""
    if config['plot_subdomain']:
        config['i_beg'], config['i_end'] = 10, 81
        config['j_beg'], config['j_end'] = 10, 90


def initialize_static_fields(ds_wrf_nc, config):
    """Initialize static WRF fields and mapping objects (run once)."""
    print('Initializing static fields...')

    # Adjust subdomain if needed
    adjust_subdomain_ranges(config)

    # Read latitude and longitude
    da_lat = wrf.getvar(ds_wrf_nc, 'lat', squeeze=False)
    wrf_lats, wrf_lons = wrf.latlon_coords(da_lat)

    # Setup cartopy mapping objects
    print('Getting cartopy mapping objects')
    cart_proj = wrf.get_cartopy(wrfin=ds_wrf_nc)
    cart_bounds = wrf.geo_bounds(var=da_lat[0, config['j_beg']:config['j_end'], config['i_beg']:config['i_end']])
    cart_xlim = wrf.cartopy_xlim(wrfin=ds_wrf_nc, geobounds=cart_bounds)
    cart_ylim = wrf.cartopy_ylim(wrfin=ds_wrf_nc, geobounds=cart_bounds)
    borders, states, oceans, lakes, _, _ = map_funcs.get_cartopy_features()

    # Build base map options dictionary
    map_opts = {
        'cart_proj': cart_proj, 'cart_xlim': cart_xlim, 'cart_ylim': cart_ylim,
        'borders': borders, 'states': states, 'oceans': oceans, 'lakes': lakes,
        'lons': wrf_lons, 'lats': wrf_lats, 'suptitle': 'Hurricane Matthew Test Case',
        'suptitle_y': config['suptitle_y'], 'lat_labels': config['lat_labels'],
        'lon_labels': config['lon_labels'], 'fontsize': config['plot_fontsize'],
        'map_x_thin': config['barb_thin'], 'map_y_thin': config['barb_thin'],
        'barb_width': config['barb_width'],
    }

    # Add station markers if requested
    if config['plot_stations']:
        text1_lat = config['mark1_lat'] + np.array([-0.20, -0.20, -0.40])
        text1_lon = config['mark1_lon'] + np.array([1.50, 3.00, 2.70])
        map_opts.update({
            'mark1_lat': config['mark1_lat'], 'mark1_lon': config['mark1_lon'],
            'text1_lab': config['text1_lab'], 'text1_lat': text1_lat, 'text1_lon': text1_lon,
            'mark1_size': config['mark1_size'], 'mark1_color': config['mark1_color']
        })

    return map_opts, wrf_lats, wrf_lons


def plot_terrain(ds_wrf_nc, map_opts, config, limits, out_dir):
    """Plot terrain height."""
    print('   Reading terrain')
    da_terrain = wrf.getvar(ds_wrf_nc, 'ter', squeeze=False)
    wrf_terrain = da_terrain.values[0, :, :]

    var_file = 'TERRAIN'
    var_name = 'Terrain Height'
    var_unit = 'm'

    min_val = np.nanmin(wrf_terrain[config['j_beg']:config['j_end'], config['i_beg']:config['i_end']])
    max_val = np.nanmax(wrf_terrain[config['j_beg']:config['j_end'], config['i_beg']:config['i_end']])

    extend = 'both'
    cmap = map_funcs.truncate_cmap(mpl.cm.terrain, minval=0.20, maxval=0.95)
    bounds = np.arange(limits['terrain']['min'], limits['terrain']['max'], limits['terrain']['int'])
    norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend=extend)

    title_l = var_name + f'\nMin: {min_val:.1f} ' + var_unit + f', Max: {max_val:.1f} ' + var_unit

    plot_opts = map_opts.copy()
    plot_opts.update({
        'fill_var': wrf_terrain, 'water_color': config['water_color'], 'extend': extend,
        'cmap': cmap, 'bounds': bounds, 'norm': norm,
        'cbar_lab': 'Model ' + var_name + ' [' + var_unit + ']',
        'fname': out_dir.joinpath('map_wrf_' + config['wrf_dom'] + '_' + var_file + '.' + config['plot_type']),
        'title_l': title_l, 'title_r': ''
    })

    map_funcs.map_plot(plot_opts)


def read_surface_winds(ds_wrf_nc):
    """Read and return 10-m wind components and speed."""
    print('   Reading 10-m wind components (rotated to earth-relative)')
    da_uv10 = wrf.getvar(ds_wrf_nc, 'uvmet10', squeeze=False)
    wrf_u10 = da_uv10.values[0, 0, :, :]
    wrf_v10 = da_uv10.values[1, 0, :, :]
    wrf_ws10 = np.sqrt(wrf_u10 ** 2 + wrf_v10 ** 2)
    return wrf_u10, wrf_v10, wrf_ws10


def plot_wind_speed_10m(wrf_u10, wrf_v10, wrf_ws10, map_opts, config, limits, constants, map_suffix, out_dir):
    """Plot 10-m wind speed with optional wind barbs."""
    var_file = 'WS10'
    var_name = '10-m Wind Speed'
    var_unit = constants['mpl_ms1']

    min_val = np.nanmin(wrf_ws10[config['j_beg']:config['j_end'], config['i_beg']:config['i_end']])
    max_val = np.nanmax(wrf_ws10[config['j_beg']:config['j_end'], config['i_beg']:config['i_end']])

    extend = 'max'
    cmap = mpl.cm.BuGn
    bounds = np.arange(limits['ws10']['min'], limits['ws10']['max'], limits['ws10']['int'])
    norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend=extend)

    # Add wind barbs if requested
    u_var, v_var = (wrf_u10, wrf_v10) if config['plot_wind_barbs_sfc'] else (None, None)
    if config['plot_wind_barbs_sfc']:
        var_file += ADD_BARBS_STRING

    title_l = var_name + f'\nMin: {min_val:.1f} ' + var_unit + f', Max: {max_val:.1f} ' + var_unit

    plot_opts = map_opts.copy()
    plot_opts.update({
        'fill_var': wrf_ws10, 'water_color': 'none', 'extend': extend,
        'cmap': cmap, 'bounds': bounds, 'norm': norm, 'u': u_var, 'v': v_var,
        'cbar_lab': var_name + ' [' + var_unit + ']',
        'fname': out_dir.joinpath('map_wrf_' + config['wrf_dom'] + '_' + var_file + map_suffix),
        'title_l': title_l, 'title_r': config['title_r']
    })

    map_funcs.map_plot(plot_opts)


def plot_sea_level_pressure(ds_wrf_nc, wrf_u10, wrf_v10, map_opts, config, limits, map_suffix, out_dir):
    """Plot sea level pressure with optional wind barbs."""
    print('   Reading sea level pressure')
    da_slp = wrf.getvar(ds_wrf_nc, 'slp', squeeze=False)
    wrf_slp = da_slp.values[0, :, :]

    var_file = 'SLP'
    var_name = 'Sea-Level Pressure'
    var_unit = 'hPa'

    min_val = np.nanmin(wrf_slp[config['j_beg']:config['j_end'], config['i_beg']:config['i_end']])
    max_val = np.nanmax(wrf_slp[config['j_beg']:config['j_end'], config['i_beg']:config['i_end']])

    extend = 'both'
    cmap = mpl.cm.viridis
    bounds = np.arange(limits['slp']['min'], limits['slp']['max'], limits['slp']['int'])
    norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend=extend)

    # Add wind barbs if requested
    u_var, v_var = (wrf_u10, wrf_v10) if config['plot_wind_barbs_sfc'] else (None, None)
    if config['plot_wind_barbs_sfc']:
        var_file += ADD_BARBS_STRING

    title_l = var_name + f'\nMin: {min_val:.1f} ' + var_unit + f', Max: {max_val:.1f} ' + var_unit

    plot_opts = map_opts.copy()
    plot_opts.update({
        'fill_var': wrf_slp, 'water_color': 'none', 'extend': extend,
        'cmap': cmap, 'bounds': bounds, 'norm': norm, 'u': u_var, 'v': v_var,
        'cbar_lab': var_name + ' [' + var_unit + ']',
        'fname': out_dir.joinpath('map_wrf_' + config['wrf_dom'] + '_' + var_file + map_suffix),
        'title_l': title_l, 'title_r': config['title_r']
    })

    map_funcs.map_plot(plot_opts)


def plot_temperature_2m(ds_wrf_nc, wrf_u10, wrf_v10, map_opts, config, limits, constants, map_suffix, out_dir):
    """Plot 2-m air temperature with optional wind barbs."""
    print('   Reading 2-m air temperature')
    da_t2 = wrf.getvar(ds_wrf_nc, 'T2', squeeze=False)
    if da_t2.attrs['units'] == 'K':
        da_t2 = da_t2 - constants['c_to_k']
        da_t2.attrs['units'] = 'degC'
    wrf_t2 = da_t2.values[0, :, :]

    var_file = 'T2'
    var_name = '2-m Air Temperature'
    var_unit = constants['deg_uni'] + 'C'

    min_val = np.nanmin(wrf_t2[config['j_beg']:config['j_end'], config['i_beg']:config['i_end']])
    max_val = np.nanmax(wrf_t2[config['j_beg']:config['j_end'], config['i_beg']:config['i_end']])

    extend = 'both'
    cmap = mpl.cm.rainbow
    bounds = np.arange(limits['t2']['min'], limits['t2']['max'], limits['t2']['int'])
    norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend=extend)

    # Add wind barbs if requested
    u_var, v_var = (wrf_u10, wrf_v10) if config['plot_wind_barbs_sfc'] else (None, None)
    if config['plot_wind_barbs_sfc']:
        var_file += ADD_BARBS_STRING

    title_l = var_name + f'\nMin: {min_val:.1f} ' + var_unit + f', Max: {max_val:.1f} ' + var_unit

    plot_opts = map_opts.copy()
    plot_opts.update({
        'fill_var': wrf_t2, 'water_color': 'none', 'extend': extend,
        'cmap': cmap, 'bounds': bounds, 'norm': norm, 'u': u_var, 'v': v_var,
        'cbar_lab': var_name + ' [' + var_unit + ']',
        'fname': out_dir.joinpath('map_wrf_' + config['wrf_dom'] + '_' + var_file + map_suffix),
        'title_l': title_l, 'title_r': config['title_r']
    })

    map_funcs.map_plot(plot_opts)


def plot_humidity_2m(ds_wrf_nc, wrf_u10, wrf_v10, map_opts, config, limits, map_suffix, out_dir):
    """Plot 2-m relative humidity with optional wind barbs."""
    print('   Reading 2-m relative humidity')
    da_rh2 = wrf.getvar(ds_wrf_nc, 'rh2', squeeze=False)
    wrf_rh2 = da_rh2.values[0, :, :]

    var_file = 'RH2'
    var_name = '2-m Relative Humidity'
    var_unit = '%'

    min_val = np.nanmin(wrf_rh2[config['j_beg']:config['j_end'], config['i_beg']:config['i_end']])
    max_val = np.nanmax(wrf_rh2[config['j_beg']:config['j_end'], config['i_beg']:config['i_end']])

    extend = 'max'
    cmap = mpl.cm.YlGnBu
    bounds = np.arange(limits['rh2']['min'], limits['rh2']['max'], limits['rh2']['int'])
    norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend=extend)

    # Add wind barbs if requested
    u_var, v_var = (wrf_u10, wrf_v10) if config['plot_wind_barbs_sfc'] else (None, None)
    if config['plot_wind_barbs_sfc']:
        var_file += ADD_BARBS_STRING

    title_l = var_name + f'\nMin: {min_val:.1f}' + var_unit + f', Max: {max_val:.1f}' + var_unit

    plot_opts = map_opts.copy()
    plot_opts.update({
        'fill_var': wrf_rh2, 'water_color': 'none', 'extend': extend,
        'cmap': cmap, 'bounds': bounds, 'norm': norm, 'u': u_var, 'v': v_var,
        'cbar_lab': var_name + ' [' + var_unit + ']',
        'fname': out_dir.joinpath('map_wrf_' + config['wrf_dom'] + '_' + var_file + map_suffix),
        'title_l': title_l, 'title_r': config['title_r']
    })

    map_funcs.map_plot(plot_opts)


def plot_precipitation(ds_wrf_nc, wrf_u10, wrf_v10, map_opts, config, limits, constants, map_suffix, out_dir):
    """Plot accumulated precipitation."""
    print('   Reading accumulated rainfall')
    da_rainc = wrf.getvar(ds_wrf_nc, 'RAINC', squeeze=False)
    da_rainnc = wrf.getvar(ds_wrf_nc, 'RAINNC', squeeze=False)
    wrf_rain = da_rainc.values[0, :, :] + da_rainnc.values[0, :, :]

    # Mask RAIN=0.0 for plotting
    wrf_rain_plot = np.ma.masked_equal(np.where(wrf_rain == 0.0, constants['missing_val'], wrf_rain),
                                       constants['missing_val'])

    var_file = 'RAIN'
    var_name = 'Accumulated Precipitation'
    var_unit = 'mm'

    min_val = np.nanmin(wrf_rain[config['j_beg']:config['j_end'], config['i_beg']:config['i_end']])
    max_val = np.nanmax(wrf_rain[config['j_beg']:config['j_end'], config['i_beg']:config['i_end']])

    extend = 'max'
    cmap = mpl.cm.GnBu
    bounds = np.arange(limits['rain']['min'], limits['rain']['max'], limits['rain']['int'])
    norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend=extend)

    # Add wind barbs if requested
    u_var, v_var = (wrf_u10, wrf_v10) if config['plot_wind_barbs_sfc'] else (None, None)
    if config['plot_wind_barbs_sfc']:
        var_file += ADD_BARBS_STRING

    title_l = var_name + f'\nMin: {min_val:.1f} ' + var_unit + f', Max: {max_val:.1f} ' + var_unit

    plot_opts = map_opts.copy()
    plot_opts.update({
        'fill_var': wrf_rain_plot, 'water_color': 'none', 'extend': extend,
        'cmap': cmap, 'bounds': bounds, 'norm': norm, 'u': u_var, 'v': v_var,
        'cbar_lab': var_name + ' [' + var_unit + ']',
        'fname': out_dir.joinpath('map_wrf_' + config['wrf_dom'] + '_' + var_file + map_suffix),
        'title_l': title_l, 'title_r': config['title_r']
    })

    map_funcs.map_plot(plot_opts)


def plot_radar_reflectivity(ds_wrf_nc, wrf_u10, wrf_v10, map_opts, config, constants, radar_data, map_suffix, out_dir):
    """Plot radar reflectivity with optional wind barbs."""
    print('   Reading radar reflectivity')
    da_refl = wrf.getvar(ds_wrf_nc, 'dbz', squeeze=False)
    wrf_refl = da_refl.values[0, 0, :, :]

    # Mask REFL <= 0.0 for plotting
    wrf_refl_plot = np.ma.masked_equal(np.where(wrf_refl <= 0.0, constants['missing_val'], wrf_refl),
                                       constants['missing_val'])

    var_file = 'REFL'
    var_name = 'Radar Reflectivity'
    var_unit = 'dBZ'

    min_val = np.nanmin(wrf_refl[config['j_beg']:config['j_end'], config['i_beg']:config['i_end']])
    max_val = np.nanmax(wrf_refl[config['j_beg']:config['j_end'], config['i_beg']:config['i_end']])

    extend = 'max'
    refl_rgb = radar_data['cmap_radar']
    bounds = radar_data['bounds_radar']
    cmap, norm = mpl.colors.from_levels_and_colors(bounds, refl_rgb, extend=extend)

    cbar_lab = var_name + ' [' + var_unit + ']'

    # Add wind barbs if requested
    u_var, v_var = (wrf_u10, wrf_v10) if config['plot_wind_barbs_sfc'] else (None, None)
    if config['plot_wind_barbs_sfc']:
        var_file += ADD_BARBS_STRING
        var_name += '; 10-m Barbs'

    title_l = var_name + f'\nMin: {min_val:.1f} ' + var_unit + f', Max: {max_val:.1f} ' + var_unit

    plot_opts = map_opts.copy()
    plot_opts.update({
        'fill_var': wrf_refl_plot, 'water_color': 'none', 'extend': extend,
        'cmap': cmap, 'bounds': bounds, 'norm': norm, 'u': u_var, 'v': v_var,
        'cbar_lab': cbar_lab,
        'fname': out_dir.joinpath('map_wrf_' + config['wrf_dom'] + '_' + var_file + map_suffix),
        'title_l': title_l, 'title_r': config['title_r']
    })

    map_funcs.map_plot(plot_opts)


def check_and_open_zlev_file(wrf_fname_zlev):
    """Check if z-level file exists and open it."""
    try:
        if not wrf_fname_zlev.is_file():
            print(f'WARNING: File {wrf_fname_zlev} does not exist. Skipping upper-level winds.')
            return None
    except FileNotFoundError:
        print(f'WARNING: File {wrf_fname_zlev} does not exist. Skipping upper-level winds.')
        return None

    print(f'Reading {wrf_fname_zlev}')
    return netCDF4.Dataset(wrf_fname_zlev, mode='r')


def plot_upper_level_winds(wrf_fname_zlev, map_opts, config, limits, constants, map_suffix, out_dir):
    """Plot 100-m wind speed from z-level file with optional wind barbs."""
    ds_wrf_zlev_nc = check_and_open_zlev_file(wrf_fname_zlev)
    if ds_wrf_zlev_nc is None:
        return

    try:
        wrf_z_zlev = wrf.getvar(ds_wrf_zlev_nc, 'Z_ZL', squeeze=False)

        # Find 100-m level
        ind_z = np.where(wrf_z_zlev == -100)[0][0]
        wrf_ws100 = wrf.getvar(ds_wrf_zlev_nc, 'S_ZL', squeeze=False)[0, ind_z, :, :]

        var_file = 'WS100'
        var_name = '100-m Wind Speed'
        var_unit = constants['mpl_ms1']

        min_val = np.nanmin(wrf_ws100[config['j_beg']:config['j_end'], config['i_beg']:config['i_end']])
        max_val = np.nanmax(wrf_ws100[config['j_beg']:config['j_end'], config['i_beg']:config['i_end']])

        extend = 'max'
        cmap = mpl.cm.BuGn
        bounds = np.arange(limits['ws10']['min'], limits['ws10']['max'], limits['ws10']['int'])
        norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend=extend)

        # Add wind barbs if requested
        u_var, v_var = None, None
        if config['plot_wind_barbs_upr']:
            var_file += ADD_BARBS_STRING
            var_name += '; Barbs'
            wrf_u100 = wrf.getvar(ds_wrf_zlev_nc, 'U_ZL', squeeze=False).values[0, ind_z, :, :]
            wrf_v100 = wrf.getvar(ds_wrf_zlev_nc, 'V_ZL', squeeze=False).values[0, ind_z, :, :]
            u_var, v_var = wrf_u100, wrf_v100

        title_l = var_name + f'\nMin: {min_val:.1f} ' + var_unit + f', Max: {max_val:.1f} ' + var_unit

        plot_opts = map_opts.copy()
        plot_opts.update({
            'fill_var': wrf_ws100, 'water_color': 'none', 'extend': extend,
            'cmap': cmap, 'bounds': bounds, 'norm': norm, 'u': u_var, 'v': v_var,
            'cbar_lab': var_name + ' [' + var_unit + ']',
            'fname': out_dir.joinpath('map_wrf_' + config['wrf_dom'] + '_' + var_file + map_suffix),
            'title_l': title_l, 'title_r': config['title_r']
        })

        map_funcs.map_plot(plot_opts)

    finally:
        ds_wrf_zlev_nc.close()


def determine_read_zlev(config):
    """Determine if z-level file needs to be read."""
    return config['plot_ws100']


def parse_time_ranges(script_config_opts, formats):
    """Parse and build time ranges for processing."""
    cycle_dt_first = pd.to_datetime(script_config_opts['cycle_dt_first'], format=formats['fmt_yyyymmdd_hh'])
    cycle_dt_last = pd.to_datetime(script_config_opts['cycle_dt_last'], format=formats['fmt_yyyymmdd_hh'])
    cycle_dt_all = pd.date_range(start=cycle_dt_first, end=cycle_dt_last,
                                 freq=str(script_config_opts['cycle_stride_h']) + 'h')
    return cycle_dt_all


def build_valid_times(cycle_dt, script_config_opts):
    """Build array of valid times for a forecast cycle."""
    beg_lead_h = int(script_config_opts['beg_lead_time'].split(':')[0])
    beg_lead_m = int(script_config_opts['beg_lead_time'].split(':')[1])
    end_lead_h = int(script_config_opts['end_lead_time'].split(':')[0])
    end_lead_m = int(script_config_opts['end_lead_time'].split(':')[1])
    valid_dt_beg = cycle_dt + dt.timedelta(hours=beg_lead_h, minutes=beg_lead_m)
    valid_dt_end = cycle_dt + dt.timedelta(hours=end_lead_h, minutes=end_lead_m)
    valid_dt_all = pd.date_range(start=valid_dt_beg, end=valid_dt_end,
                                 freq=str(script_config_opts['str_lead_time']) + 'min')
    return valid_dt_all


def process_forecast_cycle(cycle_dt, script_config_opts, static_data, config, constants, formats, limits, radar_data):
    """Process a single forecast cycle (all valid times)."""
    cycle_dt_str = cycle_dt.strftime(formats['fmt_yyyymmdd_hh'])
    cycle_dt_plot = cycle_dt.strftime(formats['fmt_time_plot'])
    start_time_plot = 'Start: ' + cycle_dt_plot
    wrf_dir = script_config_opts['wrf_dir_parent'].joinpath(cycle_dt_str)
    out_dir = script_config_opts['out_dir_parent'].joinpath(cycle_dt_str, 'plots')

    # Build valid time array
    valid_dt_all = build_valid_times(cycle_dt, script_config_opts)
    n_valid = len(valid_dt_all)

    # Loop over valid times (this assumes one output time per file)
    for vv in range(n_valid):
        valid_dt = valid_dt_all[vv]
        process_valid_time(valid_dt, vv, wrf_dir, out_dir, start_time_plot, static_data,
                           config, constants, formats, limits, radar_data, script_config_opts)


def process_valid_time(valid_dt, time_index, wrf_dir, out_dir, start_time_plot, static_data,
                       config, constants, formats, limits, radar_data, script_config_opts):
    """Process a single valid time within a forecast cycle."""
    valid_dt_wrf = valid_dt.strftime(formats['fmt_wrf_dt'])
    valid_dt_plot = valid_dt.strftime(formats['fmt_time_plot'])
    valid_dt_file = valid_dt.strftime(formats['fmt_time_file'])
    valid_time_plot = 'Valid: ' + valid_dt_plot

    config['title_r'] = start_time_plot + '\n' + valid_time_plot
    config['wrf_dom'] = 'd0' + script_config_opts['domain']
    map_suffix = '_' + valid_dt_file + '.' + config['plot_type']

    wrf_fname = wrf_dir.joinpath('wrfout_' + config['wrf_dom'] + '_' + valid_dt_wrf)
    wrf_fname_zlev = wrf_dir.joinpath('wrfout_zlev_' + config['wrf_dom'] + '_' + valid_dt_wrf)

    if not check_wrf_file_exists(wrf_fname):
        return

    print(f'Reading {wrf_fname}')
    ds_wrf_nc = netCDF4.Dataset(wrf_fname, mode='r')

    try:
        # Initialize static fields on first iteration
        if not static_data:
            map_opts, wrf_lats, wrf_lons = initialize_static_fields(ds_wrf_nc, config)
            static_data['map_opts'] = map_opts
            static_data['wrf_lats'] = wrf_lats
            static_data['wrf_lons'] = wrf_lons

            # Plot terrain (only once)
            if config['plot_terrain']:
                plot_terrain(ds_wrf_nc, map_opts, config, limits, out_dir)

        # Make the water color transparent for all subsequent plots
        static_data['map_opts']['water_color'] = 'none'

        plot_all_variables(ds_wrf_nc, wrf_fname_zlev, static_data['map_opts'], config,
                           limits, constants, radar_data, map_suffix, out_dir, time_index)

    finally:
        ds_wrf_nc.close()


def check_wrf_file_exists(wrf_fname):
    """Check if WRF file exists."""
    try:
        if not wrf_fname.is_file():
            print(f'WARNING: File {wrf_fname} does not exist. Continuing to the next valid time.')
            return False
    except FileNotFoundError:
        print(f'WARNING: File {wrf_fname} does not exist. Continuing to the next valid time.')
        return False
    return True


def plot_all_variables(ds_wrf_nc, wrf_fname_zlev, map_opts, config, limits, constants,
                       radar_data, map_suffix, out_dir, time_index):
    """Coordinate plotting of all requested variables."""
    # Read surface winds if needed for any plots
    wrf_u10, wrf_v10, wrf_ws10 = None, None, None
    if config['plot_wind_barbs_sfc'] or config['plot_ws10']:
        wrf_u10, wrf_v10, wrf_ws10 = read_surface_winds(ds_wrf_nc)

    # Plot 10-m wind speed
    if config['plot_ws10'] and wrf_ws10 is not None:
        plot_wind_speed_10m(wrf_u10, wrf_v10, wrf_ws10, map_opts, config, limits, constants, map_suffix, out_dir)

    # Plot sea level pressure
    if config['plot_slp']:
        plot_sea_level_pressure(ds_wrf_nc, wrf_u10, wrf_v10, map_opts, config, limits, map_suffix, out_dir)

    # Plot 2-m temperature
    if config['plot_t2']:
        plot_temperature_2m(ds_wrf_nc, wrf_u10, wrf_v10, map_opts, config, limits, constants, map_suffix, out_dir)

    # Plot 2-m relative humidity
    if config['plot_rh2']:
        plot_humidity_2m(ds_wrf_nc, wrf_u10, wrf_v10, map_opts, config, limits, map_suffix, out_dir)

    # Plot precipitation (skip first time step)
    if config['plot_rain'] and time_index > 0:
        plot_precipitation(ds_wrf_nc, wrf_u10, wrf_v10, map_opts, config, limits, constants, map_suffix, out_dir)

    # Plot radar reflectivity (skip first time step)
    if config['plot_refl'] and time_index > 0:
        plot_radar_reflectivity(ds_wrf_nc, wrf_u10, wrf_v10, map_opts, config, constants, radar_data, map_suffix,
                                out_dir)

    # Plot upper-level winds if requested
    if determine_read_zlev(config):
        plot_upper_level_winds(wrf_fname_zlev, map_opts, config, limits, constants, map_suffix, out_dir)


def main(script_config_opts):
    """Main function with reduced complexity."""
    # Setup configuration and constants
    config = setup_plot_configuration()
    limits = setup_plot_limits()
    constants, formats, radar_data = setup_constants_and_formats()

    # Parse time ranges
    cycle_dt_all = parse_time_ranges(script_config_opts, formats)
    n_cycles = len(cycle_dt_all)

    # Initialize static data storage
    static_data = {}

    # Loop over forecast cycles/initializations
    for cc in range(n_cycles):
        cycle_dt = cycle_dt_all[cc]
        process_forecast_cycle(cycle_dt, script_config_opts, static_data, config,
                               constants, formats, limits, radar_data)


def parse_args():
    """Parse command line arguments."""
    parser = argparse.ArgumentParser()
    parser.add_argument('-w', '--wrf_dir_parent', default='/data/input/wrf',
                        help='string specifying the directory path to the parent WRF output directories, '
                             'above any experiment or cycle datetime subdirectories (default: /data/input/wrf)')
    parser.add_argument('-o', '--out_dir_parent', default='/data/output/wrf',
                        help='string specifying the directory path to the parent plot directories (default: /data/output/wrf)')
    parser.add_argument('-f', '--cycle_dt_first', default='20161006_00',
                        help='beginning date/time of first WRF simulation [YYYYMMDD_HH] (default: 20161006_00)')
    parser.add_argument('-l', '--cycle_dt_last', default=None,
                        help='beginning date/time of last WRF simulation [YYYYMMDD_HH]')
    parser.add_argument('-i', '--cycle_stride_h', default=24, type=int,
                        help='stride in hours between cycles (default: 24)')
    parser.add_argument('-b', '--beg_lead_time', default='00:00',
                        help='beginning lead time for plotting WRF simulations [HH:MM] (default: 00:00)')
    parser.add_argument('-e', '--end_lead_time', default='48:00',
                        help='ending lead time for plotting WRF simulations [HH:MM] (default: 48:00)')
    parser.add_argument('-s', '--str_lead_time', default=180, type=int,
                        help='stride to create plots every N minutes (default: 180)')
    parser.add_argument('-d', '--domain', default='1', help='WRF domain number to be plotted (default: 1)')

    args = parser.parse_args()
    wrf_dir_parent = args.wrf_dir_parent
    out_dir_parent = args.out_dir_parent
    cycle_dt_first = args.cycle_dt_first
    cycle_dt_last = args.cycle_dt_last
    cycle_stride_h = args.cycle_stride_h
    beg_lead_time = args.beg_lead_time
    end_lead_time = args.end_lead_time
    str_lead_time = args.str_lead_time
    domain = args.domain

    if out_dir_parent is None:
        out_dir_parent = wrf_dir_parent

    # Make both paths into pathlib objects
    wrf_dir_parent = pathlib.Path(wrf_dir_parent)
    out_dir_parent = pathlib.Path(out_dir_parent)

    # Validate input formats
    validate_datetime_input(cycle_dt_first, 'init_dt_first', parser)

    if cycle_dt_last is not None:
        validate_datetime_input(cycle_dt_last, 'init_dt_last', parser)
    else:
        cycle_dt_last = cycle_dt_first

    validate_time_input(beg_lead_time, '-b (beg_lead_time)', parser)
    validate_time_input(end_lead_time, '-e (end_lead_time)', parser)

    # Put all these configuration options into a dictionary
    script_config_opts = {
        'wrf_dir_parent': wrf_dir_parent,
        'out_dir_parent': out_dir_parent,
        'cycle_dt_first': cycle_dt_first,
        'cycle_dt_last': cycle_dt_last,
        'cycle_stride_h': cycle_stride_h,
        'beg_lead_time': beg_lead_time,
        'end_lead_time': end_lead_time,
        'str_lead_time': str_lead_time,
        'domain': domain,
    }

    return script_config_opts


def validate_datetime_input(dt_str, arg_name, parser):
    """Validate datetime string format."""
    if len(dt_str) != 11:
        print(f'ERROR! Incorrect length for positional argument {arg_name}. Exiting!')
        parser.print_help()
        sys.exit()
    elif dt_str[8] != '_':
        print(f'ERROR! Incorrect format for positional argument {arg_name}. Exiting!')
        parser.print_help()
        sys.exit()


def validate_time_input(time_str, arg_name, parser):
    """Validate time string format."""
    if len(time_str) != 5:
        print(f'ERROR! Incorrect length for optional argument {arg_name}. Exiting!')
        parser.print_help()
        sys.exit()
    elif time_str[2] != ':':
        print(f'ERROR! Incorrect format for optional argument {arg_name}. Exiting!')
        parser.print_help()
        sys.exit()


if __name__ == '__main__':
    now_time_beg = dt.datetime.now(dt.UTC)

    script_config_opts = parse_args()
    main(script_config_opts)

    now_time_end = dt.datetime.now(dt.UTC)
    run_time_tot = now_time_end - now_time_beg
    now_time_beg_str = now_time_beg.strftime('%Y-%m-%d %H:%M:%S')
    now_time_end_str = now_time_end.strftime('%Y-%m-%d %H:%M:%S')
    print('\nScript completed successfully.')
    print('   Beg time: ' + now_time_beg_str)
    print('   End time: ' + now_time_end_str)
    print('   Run time: ' + str(run_time_tot) + '\n')

Running METplus

Pass the use case configuration file to the run_metplus.py script along with any user-specific system configuration files if desired:

run_metplus.py /path/to/METplus/parm/use_cases/model_applications/tc_and_extra_tc/PointStat_fcstWRF_obsMADIS_hurricane_matthew.conf /path/to/user_system.conf

See Running METplus for more information.

Expected Output

A successful run will output the following both to the screen and to the logfile:

INFO: METplus has successfully finished running.

Refer to the value set for OUTPUT_BASE to find where the output data was generated. Output for this use case will be found in {OUTPUT_BASE} and will contain the following files:

met_plot/cnt_reformatted.data
met_plot/line_T2_RMSE_BCMSE.png
point_stat/point_stat_surface_000000L_20161006_000000V.stat
point_stat/point_stat_surface_000000L_20161006_000000V_cnt.txt
point_stat/point_stat_surface_000000L_20161006_000000V_vcnt.txt
point_stat/point_stat_surface_030000L_20161006_030000V.stat
point_stat/point_stat_surface_030000L_20161006_030000V_cnt.txt
point_stat/point_stat_surface_030000L_20161006_030000V_vcnt.txt
point_stat/point_stat_surface_060000L_20161006_060000V.stat
point_stat/point_stat_surface_060000L_20161006_060000V_cnt.txt
point_stat/point_stat_surface_060000L_20161006_060000V_vcnt.txt
point_stat/point_stat_surface_090000L_20161006_090000V.stat
point_stat/point_stat_surface_090000L_20161006_090000V_cnt.txt
point_stat/point_stat_surface_090000L_20161006_090000V_vcnt.txt
point_stat/point_stat_surface_120000L_20161006_120000V.stat
point_stat/point_stat_surface_120000L_20161006_120000V_cnt.txt
point_stat/point_stat_surface_120000L_20161006_120000V_vcnt.txt
point_stat/point_stat_surface_150000L_20161006_150000V.stat
point_stat/point_stat_surface_150000L_20161006_150000V_cnt.txt
point_stat/point_stat_surface_150000L_20161006_150000V_vcnt.txt
point_stat/point_stat_surface_180000L_20161006_180000V.stat
point_stat/point_stat_surface_180000L_20161006_180000V_cnt.txt
point_stat/point_stat_surface_180000L_20161006_180000V_vcnt.txt
point_stat/point_stat_surface_210000L_20161006_210000V.stat
point_stat/point_stat_surface_210000L_20161006_210000V_cnt.txt
point_stat/point_stat_surface_210000L_20161006_210000V_vcnt.txt
point_stat/point_stat_surface_240000L_20161007_000000V.stat
point_stat/point_stat_surface_240000L_20161007_000000V_cnt.txt
point_stat/point_stat_surface_240000L_20161007_000000V_vcnt.txt
point_stat/point_stat_surface_270000L_20161007_030000V.stat
point_stat/point_stat_surface_270000L_20161007_030000V_cnt.txt
point_stat/point_stat_surface_270000L_20161007_030000V_vcnt.txt
point_stat/point_stat_surface_300000L_20161007_060000V.stat
point_stat/point_stat_surface_300000L_20161007_060000V_cnt.txt
point_stat/point_stat_surface_300000L_20161007_060000V_vcnt.txt
point_stat/point_stat_surface_330000L_20161007_090000V.stat
point_stat/point_stat_surface_330000L_20161007_090000V_cnt.txt
point_stat/point_stat_surface_330000L_20161007_090000V_vcnt.txt
point_stat/point_stat_surface_360000L_20161007_120000V.stat
point_stat/point_stat_surface_360000L_20161007_120000V_cnt.txt
point_stat/point_stat_surface_360000L_20161007_120000V_vcnt.txt
point_stat/point_stat_surface_390000L_20161007_150000V.stat
point_stat/point_stat_surface_390000L_20161007_150000V_cnt.txt
point_stat/point_stat_surface_390000L_20161007_150000V_vcnt.txt
point_stat/point_stat_surface_420000L_20161007_180000V.stat
point_stat/point_stat_surface_420000L_20161007_180000V_cnt.txt
point_stat/point_stat_surface_420000L_20161007_180000V_vcnt.txt
point_stat/point_stat_surface_450000L_20161007_210000V.stat
point_stat/point_stat_surface_450000L_20161007_210000V_cnt.txt
point_stat/point_stat_surface_450000L_20161007_210000V_vcnt.txt
point_stat/point_stat_surface_480000L_20161008_000000V.stat
point_stat/point_stat_surface_480000L_20161008_000000V_cnt.txt
point_stat/point_stat_surface_480000L_20161008_000000V_vcnt.txt
point_stat/point_stat_upper_air_000000L_20161006_000000V.stat
point_stat/point_stat_upper_air_000000L_20161006_000000V_cnt.txt
point_stat/point_stat_upper_air_000000L_20161006_000000V_vcnt.txt
point_stat/point_stat_upper_air_030000L_20161006_030000V.stat
point_stat/point_stat_upper_air_030000L_20161006_030000V_cnt.txt
point_stat/point_stat_upper_air_030000L_20161006_030000V_vcnt.txt
point_stat/point_stat_upper_air_060000L_20161006_060000V.stat
point_stat/point_stat_upper_air_060000L_20161006_060000V_cnt.txt
point_stat/point_stat_upper_air_060000L_20161006_060000V_vcnt.txt
point_stat/point_stat_upper_air_090000L_20161006_090000V.stat
point_stat/point_stat_upper_air_090000L_20161006_090000V_cnt.txt
point_stat/point_stat_upper_air_090000L_20161006_090000V_vcnt.txt
point_stat/point_stat_upper_air_120000L_20161006_120000V.stat
point_stat/point_stat_upper_air_120000L_20161006_120000V_cnt.txt
point_stat/point_stat_upper_air_120000L_20161006_120000V_vcnt.txt
point_stat/point_stat_upper_air_150000L_20161006_150000V.stat
point_stat/point_stat_upper_air_150000L_20161006_150000V_cnt.txt
point_stat/point_stat_upper_air_150000L_20161006_150000V_vcnt.txt
point_stat/point_stat_upper_air_180000L_20161006_180000V.stat
point_stat/point_stat_upper_air_180000L_20161006_180000V_cnt.txt
point_stat/point_stat_upper_air_180000L_20161006_180000V_vcnt.txt
point_stat/point_stat_upper_air_210000L_20161006_210000V.stat
point_stat/point_stat_upper_air_210000L_20161006_210000V_cnt.txt
point_stat/point_stat_upper_air_210000L_20161006_210000V_vcnt.txt
point_stat/point_stat_upper_air_240000L_20161007_000000V.stat
point_stat/point_stat_upper_air_240000L_20161007_000000V_cnt.txt
point_stat/point_stat_upper_air_240000L_20161007_000000V_vcnt.txt
point_stat/point_stat_upper_air_270000L_20161007_030000V.stat
point_stat/point_stat_upper_air_270000L_20161007_030000V_cnt.txt
point_stat/point_stat_upper_air_270000L_20161007_030000V_vcnt.txt
point_stat/point_stat_upper_air_300000L_20161007_060000V.stat
point_stat/point_stat_upper_air_300000L_20161007_060000V_cnt.txt
point_stat/point_stat_upper_air_300000L_20161007_060000V_vcnt.txt
point_stat/point_stat_upper_air_330000L_20161007_090000V.stat
point_stat/point_stat_upper_air_330000L_20161007_090000V_cnt.txt
point_stat/point_stat_upper_air_330000L_20161007_090000V_vcnt.txt
point_stat/point_stat_upper_air_360000L_20161007_120000V.stat
point_stat/point_stat_upper_air_360000L_20161007_120000V_cnt.txt
point_stat/point_stat_upper_air_360000L_20161007_120000V_vcnt.txt
point_stat/point_stat_upper_air_390000L_20161007_150000V.stat
point_stat/point_stat_upper_air_390000L_20161007_150000V_cnt.txt
point_stat/point_stat_upper_air_390000L_20161007_150000V_vcnt.txt
point_stat/point_stat_upper_air_420000L_20161007_180000V.stat
point_stat/point_stat_upper_air_420000L_20161007_180000V_cnt.txt
point_stat/point_stat_upper_air_420000L_20161007_180000V_vcnt.txt
point_stat/point_stat_upper_air_450000L_20161007_210000V.stat
point_stat/point_stat_upper_air_450000L_20161007_210000V_cnt.txt
point_stat/point_stat_upper_air_450000L_20161007_210000V_vcnt.txt
point_stat/point_stat_upper_air_480000L_20161008_000000V.stat
point_stat/point_stat_upper_air_480000L_20161008_000000V_cnt.txt
point_stat/point_stat_upper_air_480000L_20161008_000000V_vcnt.txt
wrf_plot/20161006_00/plots/map_wrf_d01_RAIN+barbs_20161006_0300.png
wrf_plot/20161006_00/plots/map_wrf_d01_RAIN+barbs_20161006_0600.png
wrf_plot/20161006_00/plots/map_wrf_d01_RAIN+barbs_20161006_0900.png
wrf_plot/20161006_00/plots/map_wrf_d01_RAIN+barbs_20161006_1200.png
wrf_plot/20161006_00/plots/map_wrf_d01_RAIN+barbs_20161006_1500.png
wrf_plot/20161006_00/plots/map_wrf_d01_RAIN+barbs_20161006_1800.png
wrf_plot/20161006_00/plots/map_wrf_d01_RAIN+barbs_20161006_2100.png
wrf_plot/20161006_00/plots/map_wrf_d01_RAIN+barbs_20161007_0000.png
wrf_plot/20161006_00/plots/map_wrf_d01_RAIN+barbs_20161007_0300.png
wrf_plot/20161006_00/plots/map_wrf_d01_RAIN+barbs_20161007_0600.png
wrf_plot/20161006_00/plots/map_wrf_d01_RAIN+barbs_20161007_0900.png
wrf_plot/20161006_00/plots/map_wrf_d01_RAIN+barbs_20161007_1200.png
wrf_plot/20161006_00/plots/map_wrf_d01_RAIN+barbs_20161007_1500.png
wrf_plot/20161006_00/plots/map_wrf_d01_RAIN+barbs_20161007_1800.png
wrf_plot/20161006_00/plots/map_wrf_d01_RAIN+barbs_20161007_2100.png
wrf_plot/20161006_00/plots/map_wrf_d01_RAIN+barbs_20161008_0000.png
wrf_plot/20161006_00/plots/map_wrf_d01_REFL+barbs_20161006_0300.png
wrf_plot/20161006_00/plots/map_wrf_d01_REFL+barbs_20161006_0600.png
wrf_plot/20161006_00/plots/map_wrf_d01_REFL+barbs_20161006_0900.png
wrf_plot/20161006_00/plots/map_wrf_d01_REFL+barbs_20161006_1200.png
wrf_plot/20161006_00/plots/map_wrf_d01_REFL+barbs_20161006_1500.png
wrf_plot/20161006_00/plots/map_wrf_d01_REFL+barbs_20161006_1800.png
wrf_plot/20161006_00/plots/map_wrf_d01_REFL+barbs_20161006_2100.png
wrf_plot/20161006_00/plots/map_wrf_d01_REFL+barbs_20161007_0000.png
wrf_plot/20161006_00/plots/map_wrf_d01_REFL+barbs_20161007_0300.png
wrf_plot/20161006_00/plots/map_wrf_d01_REFL+barbs_20161007_0600.png
wrf_plot/20161006_00/plots/map_wrf_d01_REFL+barbs_20161007_0900.png
wrf_plot/20161006_00/plots/map_wrf_d01_REFL+barbs_20161007_1200.png
wrf_plot/20161006_00/plots/map_wrf_d01_REFL+barbs_20161007_1500.png
wrf_plot/20161006_00/plots/map_wrf_d01_REFL+barbs_20161007_1800.png
wrf_plot/20161006_00/plots/map_wrf_d01_REFL+barbs_20161007_2100.png
wrf_plot/20161006_00/plots/map_wrf_d01_REFL+barbs_20161008_0000.png
wrf_plot/20161006_00/plots/map_wrf_d01_RH2+barbs_20161006_0000.png
wrf_plot/20161006_00/plots/map_wrf_d01_RH2+barbs_20161006_0300.png
wrf_plot/20161006_00/plots/map_wrf_d01_RH2+barbs_20161006_0600.png
wrf_plot/20161006_00/plots/map_wrf_d01_RH2+barbs_20161006_0900.png
wrf_plot/20161006_00/plots/map_wrf_d01_RH2+barbs_20161006_1200.png
wrf_plot/20161006_00/plots/map_wrf_d01_RH2+barbs_20161006_1500.png
wrf_plot/20161006_00/plots/map_wrf_d01_RH2+barbs_20161006_1800.png
wrf_plot/20161006_00/plots/map_wrf_d01_RH2+barbs_20161006_2100.png
wrf_plot/20161006_00/plots/map_wrf_d01_RH2+barbs_20161007_0000.png
wrf_plot/20161006_00/plots/map_wrf_d01_RH2+barbs_20161007_0300.png
wrf_plot/20161006_00/plots/map_wrf_d01_RH2+barbs_20161007_0600.png
wrf_plot/20161006_00/plots/map_wrf_d01_RH2+barbs_20161007_0900.png
wrf_plot/20161006_00/plots/map_wrf_d01_RH2+barbs_20161007_1200.png
wrf_plot/20161006_00/plots/map_wrf_d01_RH2+barbs_20161007_1500.png
wrf_plot/20161006_00/plots/map_wrf_d01_RH2+barbs_20161007_1800.png
wrf_plot/20161006_00/plots/map_wrf_d01_RH2+barbs_20161007_2100.png
wrf_plot/20161006_00/plots/map_wrf_d01_RH2+barbs_20161008_0000.png
wrf_plot/20161006_00/plots/map_wrf_d01_SLP+barbs_20161006_0000.png
wrf_plot/20161006_00/plots/map_wrf_d01_SLP+barbs_20161006_0300.png
wrf_plot/20161006_00/plots/map_wrf_d01_SLP+barbs_20161006_0600.png
wrf_plot/20161006_00/plots/map_wrf_d01_SLP+barbs_20161006_0900.png
wrf_plot/20161006_00/plots/map_wrf_d01_SLP+barbs_20161006_1200.png
wrf_plot/20161006_00/plots/map_wrf_d01_SLP+barbs_20161006_1500.png
wrf_plot/20161006_00/plots/map_wrf_d01_SLP+barbs_20161006_1800.png
wrf_plot/20161006_00/plots/map_wrf_d01_SLP+barbs_20161006_2100.png
wrf_plot/20161006_00/plots/map_wrf_d01_SLP+barbs_20161007_0000.png
wrf_plot/20161006_00/plots/map_wrf_d01_SLP+barbs_20161007_0300.png
wrf_plot/20161006_00/plots/map_wrf_d01_SLP+barbs_20161007_0600.png
wrf_plot/20161006_00/plots/map_wrf_d01_SLP+barbs_20161007_0900.png
wrf_plot/20161006_00/plots/map_wrf_d01_SLP+barbs_20161007_1200.png
wrf_plot/20161006_00/plots/map_wrf_d01_SLP+barbs_20161007_1500.png
wrf_plot/20161006_00/plots/map_wrf_d01_SLP+barbs_20161007_1800.png
wrf_plot/20161006_00/plots/map_wrf_d01_SLP+barbs_20161007_2100.png
wrf_plot/20161006_00/plots/map_wrf_d01_SLP+barbs_20161008_0000.png
wrf_plot/20161006_00/plots/map_wrf_d01_T2+barbs_20161006_0000.png
wrf_plot/20161006_00/plots/map_wrf_d01_T2+barbs_20161006_0300.png
wrf_plot/20161006_00/plots/map_wrf_d01_T2+barbs_20161006_0600.png
wrf_plot/20161006_00/plots/map_wrf_d01_T2+barbs_20161006_0900.png
wrf_plot/20161006_00/plots/map_wrf_d01_T2+barbs_20161006_1200.png
wrf_plot/20161006_00/plots/map_wrf_d01_T2+barbs_20161006_1500.png
wrf_plot/20161006_00/plots/map_wrf_d01_T2+barbs_20161006_1800.png
wrf_plot/20161006_00/plots/map_wrf_d01_T2+barbs_20161006_2100.png
wrf_plot/20161006_00/plots/map_wrf_d01_T2+barbs_20161007_0000.png
wrf_plot/20161006_00/plots/map_wrf_d01_T2+barbs_20161007_0300.png
wrf_plot/20161006_00/plots/map_wrf_d01_T2+barbs_20161007_0600.png
wrf_plot/20161006_00/plots/map_wrf_d01_T2+barbs_20161007_0900.png
wrf_plot/20161006_00/plots/map_wrf_d01_T2+barbs_20161007_1200.png
wrf_plot/20161006_00/plots/map_wrf_d01_T2+barbs_20161007_1500.png
wrf_plot/20161006_00/plots/map_wrf_d01_T2+barbs_20161007_1800.png
wrf_plot/20161006_00/plots/map_wrf_d01_T2+barbs_20161007_2100.png
wrf_plot/20161006_00/plots/map_wrf_d01_T2+barbs_20161008_0000.png
wrf_plot/20161006_00/plots/map_wrf_d01_TERRAIN.png
wrf_plot/20161006_00/plots/map_wrf_d01_WS10+barbs_20161006_0000.png
wrf_plot/20161006_00/plots/map_wrf_d01_WS10+barbs_20161006_0300.png
wrf_plot/20161006_00/plots/map_wrf_d01_WS10+barbs_20161006_0600.png
wrf_plot/20161006_00/plots/map_wrf_d01_WS10+barbs_20161006_0900.png
wrf_plot/20161006_00/plots/map_wrf_d01_WS10+barbs_20161006_1200.png
wrf_plot/20161006_00/plots/map_wrf_d01_WS10+barbs_20161006_1500.png
wrf_plot/20161006_00/plots/map_wrf_d01_WS10+barbs_20161006_1800.png
wrf_plot/20161006_00/plots/map_wrf_d01_WS10+barbs_20161006_2100.png
wrf_plot/20161006_00/plots/map_wrf_d01_WS10+barbs_20161007_0000.png
wrf_plot/20161006_00/plots/map_wrf_d01_WS10+barbs_20161007_0300.png
wrf_plot/20161006_00/plots/map_wrf_d01_WS10+barbs_20161007_0600.png
wrf_plot/20161006_00/plots/map_wrf_d01_WS10+barbs_20161007_0900.png
wrf_plot/20161006_00/plots/map_wrf_d01_WS10+barbs_20161007_1200.png
wrf_plot/20161006_00/plots/map_wrf_d01_WS10+barbs_20161007_1500.png
wrf_plot/20161006_00/plots/map_wrf_d01_WS10+barbs_20161007_1800.png
wrf_plot/20161006_00/plots/map_wrf_d01_WS10+barbs_20161007_2100.png
wrf_plot/20161006_00/plots/map_wrf_d01_WS10+barbs_20161008_0000.png

Keywords

Note

TCandExtraTCAppUseCase
PointStatToolUseCase
METplotpyUseCase
UserScriptUseCase
WRFFileUseCase
MADIS2NCToolUseCase
GRIB2FileUseCase
MADISFileUseCase

Navigate to the METplus Quick Search for Use Cases page to discover other similar use cases.

Gallery generated by Sphinx-Gallery