Note
Go to the end to download the full example code.
Point-Stat: Computing PBLH from AMDAR data using two methods: Theta-increase, Critical Bulk Richardson Number
model_applications/pbl/PointStat_fcstHRRR_obsAMDAR_PBLH_PyEmbed.conf
Scientific Objective
The Planetary Boundary Layer Height (PBLH) arises from a complex interaction of lower atmosphere and surface processes, and is therefore a useful metric to evaluate models. This PointStat use case computes PBLH from AMDAR aircraft data for 71 airports specified in a csv file using two methods: 1) “Theta-increase” method (Nielsen-Gammon et al., 2008, J. App. Met. Clim.), which computes PBLH by finding the lowest altitude in an aircraft profile that exceeds a specified increase in potential temperature from a base value. This theta-increase (pt_delta) value is set to 1.25 K, which is the same as the HRRR model. 2) “Critical Bulk Richardson Number” method (Seidel et al., 2012, JGR), which computes PBLH by finding the lowest altitude in an aircraft profile which exceeds a Critical Bulk Richardson Number. The name of the airport list csv file and the sounding are specified in the configuration script.
Version Added
METplus version 6.0
Datasets
Forecast: NOAA High Resolution Rapid Refresh (HRRR)
Observation: AMDAR hourly 1-d netcdf files
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 utilizes the METplus PointStat tool to compare PBLH from AMDAR data to model output. The python embedding script “calc_amdar_pblh.py” computes PBLH and sends data MET via python embedding.
METplus Workflow
Beginning time (INIT_BEG): 2022070108
End time (INIT_END): 2022070108
Increment between beginning and end times (INIT_INCREMENT): 1H
Sequence of forecast leads to process (LEAD_SEQ): 12
PointStat is called in this example. The following run times are processed:
Valid: 2022-07-01_20Z
Forecast lead: 12 hour
PointStat is run with Python embedding (calc_amdar_pblh.py). The Python embedding script calls an airport csv file.
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, i.e. parm/use_cases/model_applications/pbl/PointStat_fcstHRRR_obsAMDAR_PBLH_PyEmbed.conf
[config]
# Documentation for this use case can be found at
# https://metplus.readthedocs.io/en/latest/generated/model_applications/pbl/PointStat_fcstHRRR_obsAMDAR_PBLH_PyEmbed.html
# 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 = PointStat
###
# 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 = 2022070108
INIT_END = 2022070108
INIT_INCREMENT = 1H
LEAD_SEQ = 12
###
# USER SETTINGS HERE
# PY_VAL_TIME = Valid_Time(YYYYMMDD_HHMMSS)
# PY_SOUNDING_FLAG = ALL, ASCENTS, DESCENTS [only one value here]
# PY_AIRPORT_FILE = name of file containing airport info (without csv suffix)
###
PY_VAL_TIME = {valid?fmt=%Y%m%d_%H0000}
PY_SOUNDING_FLAG = ASCENTS
PY_AIRPORT_FILE = amdar_airports_pblh_usa_71
###
# File I/O
# https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#directory-and-filename-template-info
###
FCST_POINT_STAT_INPUT_DIR = {INPUT_BASE}/model_applications/pbl/PointStat_fcstHRRR_obsAMDAR_PBLH_PyEmbed
FCST_POINT_STAT_INPUT_TEMPLATE = {init?fmt=%Y%m%d%H}/pblh_22{init?fmt=%j}{init?fmt=%H}0000{lead?fmt=%HH}
OBS_POINT_STAT_INPUT_DIR =
OBS_POINT_STAT_INPUT_TEMPLATE = PYTHON_NUMPY= {PARM_BASE}/use_cases/model_applications/pbl/PointStat_fcstHRRR_obsAMDAR_PBLH_PyEmbed/calc_amdar_pblh.py {INPUT_BASE}/model_applications/pbl/PointStat_fcstHRRR_obsAMDAR_PBLH_PyEmbed/22{valid?fmt=%j}{valid?fmt=%H}00q.cdf
POINT_STAT_OUTPUT_DIR = {OUTPUT_BASE}/point_stat_pblh
POINT_STAT_OUTPUT_TEMPLATE = {valid?fmt=%Y%m%d}
POINT_STAT_CLIMO_MEAN_INPUT_DIR =
POINT_STAT_CLIMO_MEAN_INPUT_TEMPLATE =
POINT_STAT_CLIMO_STDEV_INPUT_DIR =
POINT_STAT_CLIMO_STDEV_INPUT_TEMPLATE =
###
# Field Info
# https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#field-info
###
POINT_STAT_ONCE_PER_FIELD = False
FCST_VAR1_NAME = HPBL
FCST_VAR1_LEVELS = L0
FCST_VAR1_THRESH = <=0, >10000
OBS_VAR1_NAME = HPBL_TI
OBS_VAR1_LEVELS = L0
OBS_VAR1_THRESH = <=0, >10000
OBS_VAR1_OPTIONS = set_attr_units = "m"
FCST_VAR2_NAME = HPBL
FCST_VAR2_LEVELS = L0
FCST_VAR2_THRESH = <=0, >10000
OBS_VAR2_NAME = HPBL_BR
OBS_VAR2_LEVELS = L0
OBS_VAR2_THRESH = <=0, >10000
OBS_VAR2_OPTIONS = set_attr_units = "m"
# PointStat Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#pointstat
###
LOG_POINT_STAT_VERBOSITY = 4
POINT_STAT_CONFIG_FILE ={PARM_BASE}/met_config/PointStatConfig_wrapped
POINT_STAT_CLIMO_MEAN_TIME_INTERP_METHOD = NEAREST
POINT_STAT_INTERP_TYPE_METHOD = BILIN
POINT_STAT_INTERP_TYPE_WIDTH = 2
POINT_STAT_OUTPUT_FLAG_SL1L2 = STAT
POINT_STAT_OUTPUT_FLAG_VL1L2 = STAT
POINT_STAT_OUTPUT_FLAG_MPR = STAT
OBS_POINT_STAT_WINDOW_BEGIN = -3600
OBS_POINT_STAT_WINDOW_END = 3600
POINT_STAT_OFFSETS = 0
MODEL = HRRR
POINT_STAT_DESC = {PY_SOUNDING_FLAG}
OBTYPE =
POINT_STAT_REGRID_TO_GRID = NONE
POINT_STAT_REGRID_METHOD = BILIN
POINT_STAT_REGRID_WIDTH = 2
POINT_STAT_OUTPUT_PREFIX = {PY_SOUNDING_FLAG}
POINT_STAT_MESSAGE_TYPE = ADPSFC
# INFO TO PASS TO PYTHON SCRIPT
[user_env_vars]
VAL_TIME = {PY_VAL_TIME}
SOUNDING_FLAG = {PY_SOUNDING_FLAG}
AIRPORT_FILE = {PARM_BASE}/use_cases/model_applications/pbl/PointStat_fcstHRRR_obsAMDAR_PBLH_PyEmbed/{PY_AIRPORT_FILE}
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}
//lapse_rate_correction = {
${METPLUS_LAPSE_RATE_CORRECTION_DICT}
//msl_agl_conversion = {
${METPLUS_MSL_AGL_CONVERSION_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 uses a Python embedding script to read input data.
parm/use_cases/model_applications/pbl/PointStat_fcstHRRR_obsAMDAR_PBLH_PyEmbed/calc_amdar_pblh.py
"""
This script reads AMDAR hourly netcdf files, computes PBLH, and sends 11-column ascii table to MET for point-stat
An airport csv file is read in containing lat, lon, gnd height, and rbox for each airport
See accompanying PointStat_fcstHRRR_obsAMDAR_PBLH_PyEmbed.conf for settings and passing in env variables here
Jason M. English, May 2025
"""
import os
import sys
from pathlib import Path
import pandas as pd
import numpy as np
import netCDF4 as nc
import math
from typing import Tuple
from collections import defaultdict
import warnings
# Suppress non-critical warnings
warnings.filterwarnings("ignore", category=DeprecationWarning)
print("Python Script:\t" + repr(sys.argv[0]))
##
## input file specified on the command line
## load the data into the numpy array
##
# ---------------------------------------------------------------------------
# ENVIRONMENT VARIABLES
# ---------------------------------------------------------------------------
val_time = os.environ.get('VAL_TIME') #'20220701_200000' #os.environ.get('VAL_TIME')
sf_include = os.environ.get('SOUNDING_FLAG') #'ALL' #os.environ.get('SOUNDING_FLAG') # 'ASCENTS', 'DESCENTS', or 'ALL'
airport_file = os.environ.get('AIRPORT_FILE') # airport file template
# ---------------------------------------------------------------------------
# CONFIGURATION & PHYSICS CONSTANTS
# ---------------------------------------------------------------------------
CONFIG = {
"alt_dp": 4,
"alt_adj_flag": True,
"pt_delta": 1.25,
"altmax_sfc": 200.,
"gap_max": 400.,
"gap_max_denom": 20.,
"cbrn": 0.5,
"zs": 40.,
"ustar": 0.3,
"beta": 100.,
}
PHYSICS = {
"SLP": 101325.0,
"GRAV": 9.80665,
"M_MASS": 0.0289644,
"R0": 8.31432,
"LR": 0.0065,
"H0": 44307.694
}
PHYSICS["EXPON"] = (-PHYSICS["GRAV"] * PHYSICS["M_MASS"]) / (PHYSICS["R0"] * -PHYSICS["LR"])
# ---------------------------------------------------------------------------
# UTILITY FUNCTIONS
# ---------------------------------------------------------------------------
def load_hourly_netcdf(filepath: Path) -> dict:
""" Load raw AMDAR netcdf file (1d arrays of all data that hour). """
with nc.Dataset(filepath, 'r') as ncf:
data = {
'tailNumber': np.array(ncf['tailNumber'][:]),
'altitude': np.array(ncf['altitude'][:]),
'latitude': np.array(ncf['latitude'][:]),
'longitude': np.array(ncf['longitude'][:]),
'temperature': np.array(ncf['temperature'][:]),
'windSpeed': np.array(ncf['windSpeed'][:]),
'windDir': np.array(ncf['windDir'][:]),
'timeObs': np.array(ncf['timeObs'][:]),
'sounding_flag': np.array(ncf['sounding_flag'][:])
}
return data
def sf_mask(sf_value: int) -> bool:
""" Return sounding flag mask so we can process the desired flights. """
if sf_value == 0: # discard cruising flights
return False
if sf_include == "ALL":
return True
if sf_include == "ASCENTS":
return sf_value == 1
if sf_include == "DESCENTS":
return sf_value == -1
def get_tail_number_string(tn_array: np.ndarray) -> np.ndarray:
"""Convert 9-character tail number array to a concatenated string."""
tnc = np.char.array(tn_array.astype(str))
return tnc[:, 0] + tnc[:, 1] + tnc[:, 2] + tnc[:, 3] + tnc[:, 4] + tnc[:, 5] + tnc[:, 6] + tnc[:, 7] + tnc[:, 8]
def compute_pressure(alt: np.ndarray, ground_level: float) -> np.ndarray:
"""Compute pressure (Pa) using the hypsometric formula."""
z = alt + ground_level
return PHYSICS["SLP"] * (1 - z / PHYSICS["H0"]) ** PHYSICS["EXPON"]
def compute_potential_temperature(temp: np.ndarray, pressure: np.ndarray) -> np.ndarray:
"""Compute potential temperature (K)."""
return temp * (PHYSICS["SLP"] / pressure) ** 0.286
# ---------------------------------------------------------------------------
# PBLH COMPUTATION METHODS
# ---------------------------------------------------------------------------
def compute_pblh_ti(alt: np.ndarray, pt: np.ndarray, lat: np.ndarray, lon: np.ndarray) -> Tuple[float, float, float]:
"""
Compute PBLH via the theta increase (TI) method.
Returns: (pblh, lat_pblh, lon_pblh)
"""
valid = alt < CONFIG["altmax_sfc"] # To find pt_min consider only alt below altmax_sfc
if not np.any(valid):
return np.nan, np.nan, np.nan
pt_min = np.nanmin(np.where(valid, pt, np.nan))
if pt_min <= 0 or pt_min >= 3040:
return np.nan, np.nan, np.nan
try:
pt_min_index = np.nonzero(pt == pt_min)[0][0]
except IndexError:
return np.nan, np.nan, np.nan
# now mask out the profile below pt_min, and search above it for pt_delta
alt_ti = alt.copy()
alt_ti[:pt_min_index] = np.nan
pt_ti = pt.copy()
pt_ti[:pt_min_index] = np.nan
pt_target = pt_min + CONFIG['pt_delta']
inds = np.nonzero(pt_ti >= pt_target)[0]
if inds.size == 0 or inds[0] == 0:
return np.nan, np.nan, np.nan
pblh_index = inds[0]
# discard this PBLH value if gap between data points is too big
alt_gap = alt_ti[pblh_index] - alt_ti[pblh_index - 1]
if alt_gap > CONFIG["gap_max"] + alt_ti[pblh_index] / CONFIG["gap_max_denom"]:
return np.nan, np.nan, np.nan
i1, i0 = inds[0], inds[0] - 1
# interpolate between data points to get PBLH
pblh = float(np.interp(pt_target, pt_ti[i0:i1+1], alt_ti[i0:i1+1]))
lat_pblh = float(np.interp(pt_target, pt_ti[i0:i1+1], lat[i0:i1+1]))
lon_pblh = float(np.interp(pt_target, pt_ti[i0:i1+1], lon[i0:i1+1]))
return pblh, lat_pblh, lon_pblh
def compute_pblh_br(alt: np.ndarray, pt: np.ndarray, lat: np.ndarray, lon: np.ndarray,
ws: np.ndarray, wd: np.ndarray) -> Tuple[float, float, float]:
"""
Compute PBLH via the Critical Bulk Richardson Number (CBRN) method.
Returns: (pblh, lat_pblh, lon_pblh)
"""
br_sfc_ind = np.argmin(np.abs(alt - CONFIG["zs"])) # Find data point closest to zs
if alt[br_sfc_ind] > CONFIG["altmax_sfc"]: # Make sure it's not too far from zs
return np.nan, np.nan, np.nan
if np.isnan(ws[br_sfc_ind]) or np.isnan(pt[br_sfc_ind]) or np.isnan(wd[br_sfc_ind]):
return np.nan, np.nan, np.nan
wd_math = (270.0 - wd) % 360.0
u = ws * np.cos(np.radians(wd_math))
v = ws * np.sin(np.radians(wd_math))
u_sfc, v_sfc, pt_sfc = u[br_sfc_ind], v[br_sfc_ind], pt[br_sfc_ind]
brn_prev, alt_prev, lat_prev, lon_prev = None, None, None, None
for i in range(br_sfc_ind+1, len(alt)):
if np.isnan(ws[i]) or np.isnan(pt[i]) or np.isnan(alt[i]):
continue
brn = (PHYSICS["GRAV"] / pt_sfc) * (pt[i] - pt_sfc) * (alt[i] - CONFIG["zs"]) / (
(u[i] - u_sfc)**2 + (v[i] - v_sfc)**2 + CONFIG['beta'] * CONFIG['ustar']**2)
if brn > CONFIG["cbrn"]:
if i == 0 or brn_prev is None:
return alt[i], lat[i], lon[i]
# discard this PBLH value if gap between data points is too big
alt_gap = alt[i] - alt_prev
if alt_gap > CONFIG["gap_max"] + alt[i] / CONFIG["gap_max_denom"]:
return np.nan, np.nan, np.nan
# interpolate between data points to get PBLH
pblh = float(np.interp(CONFIG["cbrn"], [brn_prev, brn], [alt_prev, alt[i]]))
lat_pblh = float(np.interp(CONFIG["cbrn"], [brn_prev, brn], [lat_prev, lat[i]]))
lon_pblh = float(np.interp(CONFIG["cbrn"], [brn_prev, brn], [lon_prev, lon[i]]))
return pblh, lat_pblh, lon_pblh
brn_prev, alt_prev, lat_prev, lon_prev = brn, alt[i], lat[i], lon[i]
return np.nan, np.nan, np.nan
# ---------------------------------------------------------------------------
# FLIGHT PROCESSING
# ---------------------------------------------------------------------------
def process_flight(tail: str, indices: np.ndarray, data: dict,
ground_level: float) -> dict:
"""
Process one flight/tail number. Sort the data points to be ascending. Compute PBLH via TI and BR methods.
Returns: PBLH, lat, lon for that tail number.
"""
flight_data = {k: np.take(data[k], indices, axis=0) for k in data}
sort_order = np.argsort(flight_data['altitude'])
for k in flight_data:
flight_data[k] = flight_data[k][sort_order]
alt = flight_data['altitude'] - ground_level
if CONFIG["alt_adj_flag"] and np.nanmin(alt) < 0: # if minimum altitude is negative, adjust to zero
offset = float(np.trunc(np.nanmin(alt)))
alt -= offset
pres = compute_pressure(alt, ground_level)
pt = compute_potential_temperature(flight_data['temperature'], pres)
pblh_ti, lat_ti, lon_ti = compute_pblh_ti(alt, pt, flight_data['latitude'], flight_data['longitude'])
pblh_br, lat_br, lon_br = compute_pblh_br(alt, pt, flight_data['latitude'], flight_data['longitude'],
flight_data['windSpeed'], flight_data['windDir'])
return {'tail_number': tail,
'pblh_ti': pblh_ti, 'lat_ti': lat_ti, 'lon_ti': lon_ti,
'pblh_br': pblh_br, 'lat_br': lat_br, 'lon_br': lon_br}
# ---------------------------------------------------------------------------
# CREATE DATAFRAME
# ---------------------------------------------------------------------------
def create_dataframe(sid_all, var_all, lat_all, lon_all, elev_all, pblh_all, val_time: str) -> pd.DataFrame:
""" Create an 11-column dataframe in the format required by MET. """
# (1) string: Message_Type ('ADPSFC')
# (2) string: Station_ID (AIRPORT)
# (3) string: Valid_Time(YYYYMMDD_HHMMSS)
# (4) numeric: Lat(Deg North)
# (5) numeric: Lon(Deg East)
# (6) numeric: Elevation(msl)
# (7) string: Var_Name(or GRIB_Code)
# (8) numeric: Level
# (9) numeric: Height(msl or agl)
# (10) string: QC_String
# (11) numeric: Observation_Value
df = pd.DataFrame({
'typ': ['ADPSFC']*len(sid_all), 'sid': sid_all, 'vld': [val_time]*len(sid_all),
'lat': lat_all, 'lon': lon_all, 'elv': elev_all, 'var': var_all,
'lvl': [0]*len(sid_all), 'hgt': [0]*len(sid_all), 'qc': ['AMDAR']*len(sid_all), 'obs': pblh_all
})
return df[df['obs'].notna()]
# ---------------------------------------------------------------------------
# MAIN EXECUTION
# ---------------------------------------------------------------------------
# Load hourly AMDAR netcdf file
infile = Path(sys.argv[1]) if len(sys.argv)>1 else Path('221822000q.cdf') # remove the end after debugging
data = load_hourly_netcdf(infile)
# Convert tail 2d char array to 1d string array
tn_str = get_tail_number_string(data['tailNumber'])
data['tailNumber'] = tn_str
# Apply sounding flag mask to the whole file (remove cruising flights, and asc/desc if specified)
sf = data['sounding_flag']
mask_sf = np.array([sf_mask(val) for val in sf])
for k in data:
data[k] = data[k][mask_sf]
# Load airport csv file
airport_df = pd.read_csv(airport_file + ".csv", index_col=0)
airports = airport_df.index.tolist()
# Create empty lists to append
sid_all, var_all, lat_all, lon_all, elev_all, pblh_all = [], [], [], [], [], []
# Loop over airports, masking that lat/lon box for each
for airport in airports:
code = airport_df.loc[airport, 'airport_code']
lat0 = airport_df.loc[airport, 'lat_degN']
lon0 = airport_df.loc[airport, 'lon_degE']
gnd0 = airport_df.loc[airport, 'hgt_m_MSL']
r0 = airport_df.loc[airport, 'airport_radius_deg']
mask_box = ((data['latitude'] > lat0 - r0) & (data['latitude'] < lat0 + r0) &
(data['longitude'] > lon0 - r0) & (data['longitude'] < lon0 + r0))
# Filter tail numbers within box and remove NaNs
filtered_tn = np.where(mask_box, data['tailNumber'], "nan")
valid_tails = np.array([t for t in np.unique(filtered_tn) if isinstance(t, str) and t.lower() != "nan"])
print(f"\n==Processing airport {airport} ({code}): {len(valid_tails)} flights==")
if valid_tails.size == 0:
continue
# Loop through all flights (tail numbers) within the lat/lon box of this airport
for tail in valid_tails:
indices = np.nonzero(filtered_tn == tail)[0]
print(f"Processing flight {tail} with {len(indices)} points.")
if len(indices) < CONFIG["alt_dp"]:
print(f"Flight {tail}: insufficient data points ({len(indices)}). Skipping.")
continue
flight_result = process_flight(str(tail), indices, data, gnd0)
# Append dataframe row with TI method if successful
if not np.isnan(flight_result['pblh_ti']):
sid_all.append(code)
var_all.append('HPBL_TI')
elev_all.append(gnd0)
lat_all.append(flight_result['lat_ti'])
lon_all.append(flight_result['lon_ti'])
pblh_all.append(flight_result['pblh_ti'])
# Append dataframe row with BR method if successful
if not np.isnan(flight_result['pblh_br']):
sid_all.append(code)
var_all.append('HPBL_BR')
elev_all.append(gnd0)
lat_all.append(flight_result['lat_br'])
lon_all.append(flight_result['lon_br'])
pblh_all.append(flight_result['pblh_br'])
point_data = create_dataframe(sid_all, var_all,
lat_all, lon_all,
elev_all, pblh_all,
val_time)
pd.set_option('display.max_rows', None)
print(point_data)
print(" point_data: Data Length:\t" + repr(len(point_data)))
print(" point_data: Data Type:\t" + repr(type(point_data)))
point_data = point_data.values.tolist()
#except NameError:
# print("Can't find the input file")
For more information on the basic requirements to utilize Python Embedding in METplus, please refer to the MET User’s Guide section on Python embedding.
User Scripting
User Scripting is not used in this use case.
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/PointStat_fcstHRRR_obsAMDAR_PBLH_PyEmbed.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 point_stat_pblh (relative to OUTPUT_BASE) with subdirectories for valid time (YYYYMMDD) and will contain .stat files with the following naming convention:
- convention: point_stat_{SOUNDING_FLAG}_{LEADTIME}L_{VALIDTIME}.stat
example: point_stat_ALL_120000L_20220701_200000V.stat
Keywords
Note
PointStatToolUseCase
PythonEmbeddingFileUseCase
PBLAppUseCase
Navigate to the METplus Quick Search for Use Cases page to discover other similar use cases.
sphinx_gallery_thumbnail_path = ‘_static/pbl-PointStat_fcstHRRR_obsAMDAR_PBLH_PyEmbed.png’