Note
Go to the end to download the full example code.
PointStat: Verify UFS Soil Moisture and Temperature with ISMN Observations
model_applications/land_surface/PointStat_fcstUFS_obsISMN_SoilMoistureTemp.conf
Scientific Objective
This use case examines soil moisture and temperature biases in the UFS global forecast system (GFS). The specific configuration is a GFSv17 pre-release version (HR1). Correct representation of land states are important for many processes, such as prediction of convection, near surface sensible weather (e.g. winds, humidity, temperature), and subsequent hydrological forecasting applications such as snowpack evolution and runoff generation. Here we use the International Soil Moisture Network (ISMN) soil moisture and temperature data to assess UFS forecast errors across CONUS for an example summer 60 hour forecast. We plot spatial differences (forecast-observation) to diagnose regional variations as well as a 2D probability density plot of forecast difference versus forecast value to assess the frequency of conditional differences. These diagnostics are useful for process understanding and may aid in land model development and forecast error improvement.
Version Added
METplus version 6.1
Datasets
Forecast: Global Forecast System (GFS) version 17 prototype. Global 1-degree grid including soil moisture content and soil temperature.
Observation: International Soil Moisture Network (ISMN) observations of soil moisture content and soil temperature.
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 uses ASCII2NC to process ISMN observations from their raw format and then calls PointStat to verify the forecast data using those observations. After PointStat, StatAnalysis is called to compute an aggregate mean error (ME, Bias) across all times and forecast leads. Then METplus UserScripts are used to create various graphics from the output.
METplus Workflow
Beginning time (INIT_BEG): 2020-08-05 12:00
End time (INIT_END): 2020-08-05 12:00
Increment between beginning and end times (INIT_INCREMENT): 12 Hours
Sequence of forecast leads to process (LEAD_SEQ): 60
Only a single time is used to demonstrate the workflow for this use case. For each time, ASCII2NC is used to convert the ISMN observations to NetCDF. Then PointStat is called to create matched forecast and observation pairs for each variable and level defined by the user. After PointStat, StatAnalysis is used to compute an aggregate mean error (ME, bias) for all forecast/observation pairs across all valid times and forecast leads. This bias value is used to create a plot of the ME at each ISMN site on a map, color-coded by the bias value. In addition, 2-D histogram plots are created showing the relationship between the forecast minus observation values and the forecast values.
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/land_surface/PointStat_fcstUFS_obsISMN_SoilMoistureTemp.conf
[config]
# Documentation for this use case can be found at
# https://metplus.readthedocs.io/en/latest/generated/model_applications/land_surface/PointStat_fcstUFS_obsISMN_SoilMoistureTemp.html
# For additional information, please see the METplus Users Guide.
# https://metplus.readthedocs.io/en/latest/Users_Guide
###########################################################
# General METplus Wrappers Settings
###########################################################
ISMN_NETWORK_LIST = USCRN,ARM,SNOTEL,SCAN
############################################################
# Processes to run
# https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#process-list
############################################################
PROCESS_LIST =
ASCII2NC,
PointStat,
StatAnalysis,
UserScript(plot_soilw_0_0.1_bias_map),
UserScript(plot_tsoil_0_0.1_bias_map),
UserScript(plot_soilw_0_0.1_2d_hist),
UserScript(plot_tsoil_0_0.1_2d_hist)
############################################################
# 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 = VALID
VALID_BEG = 2020080512
VALID_END = 2020080512
VALID_INCREMENT = 6H
VALID_TIME_FMT = %Y%m%d%H
LEAD_SEQ = 60
############################################################
# ASCII2NC settings for processing GDAS observations
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#pb2nc
############################################################
ASCII2NC_INPUT_DIR = {INPUT_BASE}/model_applications/land_surface/PointStat_fcstUFS_obsISMN_SoilMoistureTemp/obs/ismn/HeadVals
ASCII2NC_INPUT_TEMPLATE = SNOTEL, SCAN, USCRN, ARM
ASCII2NC_INPUTRX = .stm$
ASCII2NC_INPUT_FORMAT = ismn
ASCII2NC_OUTPUT_DIR = {OUTPUT_BASE}/ascii2nc
ASCII2NC_OUTPUT_TEMPLATE = {valid?fmt=%Y%m%d%H}_ismn.nc
LOG_ASCII2NC_VERBOSITY = 3
############################################################
# PointStat Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#pointstat
############################################################
POINT_STAT_CONFIG_FILE = {PARM_BASE}/met_config/PointStatConfig_wrapped
POINT_STAT_MESSAGE_TYPE = {ISMN_NETWORK_LIST}
POINT_STAT_OUTPUT_DIR={OUTPUT_BASE}/point_stat
POINT_STAT_OUTPUT_FLAG_MPR = BOTH
FCST_POINT_STAT_INPUT_DIR = {INPUT_BASE}/model_applications/land_surface/PointStat_fcstUFS_obsISMN_SoilMoistureTemp/fcst
FCST_POINT_STAT_INPUT_TEMPLATE = {init?fmt=%Y%m%d%H}/{init?fmt=%Y%m%d%H}.f{lead?fmt=%3H}.grib2
OBS_POINT_STAT_INPUT_DIR = {ASCII2NC_OUTPUT_DIR}
OBS_POINT_STAT_INPUT_TEMPLATE = {ASCII2NC_OUTPUT_TEMPLATE}
OBS_POINT_STAT_WINDOW_BEGIN = -10800
OBS_POINT_STAT_WINDOW_END = 10800
POINT_STAT_OBS_SUMMARY = NEAREST
LOG_POINT_STAT_VERBOSITY = 3
############################################################
# Field Info
# https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#field-info
############################################################
FCST_POINT_STAT_VAR1_NAME = SOILW
FCST_POINT_STAT_VAR1_LEVELS = Z0-0.1
FCST_POINT_STAT_VAR2_NAME = TSOIL
FCST_POINT_STAT_VAR2_LEVELS = Z0-0.1
OBS_POINT_STAT_VAR1_NAME = SOILW
OBS_POINT_STAT_VAR1_LEVELS = Z0-0.1
OBS_POINT_STAT_VAR2_NAME = TSOIL
OBS_POINT_STAT_VAR2_LEVELS = Z0-0.1
MODEL = UFS_HR1
############################################################
# StatAnaylsis Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#statanalysis
############################################################
FCST_VAR1_LEVELS = Z0.1-0
OBS_VAR1_LEVELS = Z0.1-0
FCST_VAR2_LEVELS = Z0.1-0
OBS_VAR2_LEVELS = Z0.1-0
STAT_ANALYSIS_CONFIG_FILE = {PARM_BASE}/met_config/STATAnalysisConfig_wrapped
STAT_ANALYSIS_OUTPUT_DIR = {OUTPUT_BASE}/stat_analysis
STAT_ANALYSIS_OUTPUT_FILE = {valid?fmt=%Y%m%d%H}_{valid?fmt=%Y%m%d%H}_OBS_SID_CNT.stat
STAT_ANALYSIS_JOB1 = -job aggregate_stat -line_type MPR -obtype {ISMN_NETWORK_LIST} -out_line_type CNT -by OBS_SID,OBS_VAR -set_hdr VX_MASK OBS_SID -out_stat {STAT_ANALYSIS_OUTPUT_DIR}/{STAT_ANALYSIS_OUTPUT_FILE}
MODEL1 = UFS_HR1
MODEL1_STAT_ANALYSIS_LOOKIN_DIR = {OUTPUT_BASE}/point_stat
LOG_STAT_ANALYSIS_VERBOSITY = 3
############################################################
# UserScript Settings for plot_script
# This script is used to create a graphic from the matched
# forecast and observation pairs created with PointStat
############################################################
[plot_soilw_0_0.1_bias_map]
USER_SCRIPT_RUNTIME_FREQ = RUN_ONCE
USER_SCRIPT_OUTPUT_DIR = {OUTPUT_BASE}
USER_SCRIPT_POINT_STAT_OUTPUT_DIR = {OUTPUT_BASE}/point_stat
USER_SCRIPT_STAT_ANALYSIS_INPUT_FILE = {STAT_ANALYSIS_OUTPUT_DIR}/{STAT_ANALYSIS_OUTPUT_FILE}
USER_SCRIPT_COMMAND = python3 {PARM_BASE}/use_cases/model_applications/land_surface/PointStat_fcstUFS_obsISMN_SoilMoistureTemp/plot_point_stat_bias_map_ISMN.py {USER_SCRIPT_STAT_ANALYSIS_INPUT_FILE} {USER_SCRIPT_POINT_STAT_OUTPUT_DIR} SOILW 0-0.1 {USER_SCRIPT_OUTPUT_DIR}
[plot_tsoil_0_0.1_bias_map]
USER_SCRIPT_RUNTIME_FREQ = RUN_ONCE
USER_SCRIPT_OUTPUT_DIR = {OUTPUT_BASE}
USER_SCRIPT_POINT_STAT_OUTPUT_DIR = {OUTPUT_BASE}/point_stat
USER_SCRIPT_STAT_ANALYSIS_INPUT_FILE = {STAT_ANALYSIS_OUTPUT_DIR}/{STAT_ANALYSIS_OUTPUT_FILE}
USER_SCRIPT_COMMAND = python3 {PARM_BASE}/use_cases/model_applications/land_surface/PointStat_fcstUFS_obsISMN_SoilMoistureTemp/plot_point_stat_bias_map_ISMN.py {USER_SCRIPT_STAT_ANALYSIS_INPUT_FILE} {USER_SCRIPT_POINT_STAT_OUTPUT_DIR} TSOIL 0-0.1 {USER_SCRIPT_OUTPUT_DIR}
[plot_soilw_0_0.1_2d_hist]
USER_SCRIPT_RUNTIME_FREQ = RUN_ONCE
USER_SCRIPT_OUTPUT_DIR = {OUTPUT_BASE}
USER_SCRIPT_POINT_STAT_OUTPUT_DIR = {OUTPUT_BASE}/point_stat
USER_SCRIPT_COMMAND = python3 {PARM_BASE}/use_cases/model_applications/land_surface/PointStat_fcstUFS_obsISMN_SoilMoistureTemp/plot_point_stat_bias_2d_hist_ISMN.py {USER_SCRIPT_POINT_STAT_OUTPUT_DIR} SOILW 0-0.1 {USER_SCRIPT_OUTPUT_DIR}
[plot_tsoil_0_0.1_2d_hist]
USER_SCRIPT_RUNTIME_FREQ = RUN_ONCE
USER_SCRIPT_OUTPUT_DIR = {OUTPUT_BASE}
USER_SCRIPT_POINT_STAT_OUTPUT_DIR = {OUTPUT_BASE}/point_stat
USER_SCRIPT_COMMAND = python3 {PARM_BASE}/use_cases/model_applications/land_surface/PointStat_fcstUFS_obsISMN_SoilMoistureTemp/plot_point_stat_bias_2d_hist_ISMN.py {USER_SCRIPT_POINT_STAT_OUTPUT_DIR} TSOIL 0-0.1 {USER_SCRIPT_OUTPUT_DIR}
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
Ascii2NcConfig_wrapped
////////////////////////////////////////////////////////////////////////////////
//
// Default ascii2nc configuration file
//
////////////////////////////////////////////////////////////////////////////////
//
// The parameters listed below are used to summarize the ASCII data read in
//
//
// Time periods for the summarization
// obs_var (string array) is added and works like grib_code (int array)
// when the obs name is given instead of grib_code
//
${METPLUS_TIME_SUMMARY_DICT}
//
// Mapping of input little_r report types to output message types
//
message_type_map = [
{ key = "FM-12 SYNOP"; val = "ADPSFC"; },
{ key = "FM-13 SHIP"; val = "SFCSHP"; },
{ key = "FM-15 METAR"; val = "ADPSFC"; },
{ key = "FM-18 BUOY"; val = "SFCSHP"; },
{ key = "FM-281 QSCAT"; val = "ASCATW"; },
{ key = "FM-32 PILOT"; val = "ADPUPA"; },
{ key = "FM-35 TEMP"; val = "ADPUPA"; },
{ key = "FM-88 SATOB"; val = "SATWND"; },
{ key = "FM-97 ACARS"; val = "AIRCFT"; }
];
//
// Indicate a version number for the contents of this configuration file.
// The value should generally not be modified.
//
//version = "V10.0";
tmp_dir = "${MET_TMP_DIR}";
${METPLUS_MET_CONFIG_OVERRIDES}
PointStatConfig_wrapped
////////////////////////////////////////////////////////////////////////////////
//
// Default ascii2nc configuration file
//
////////////////////////////////////////////////////////////////////////////////
//
// The parameters listed below are used to summarize the ASCII data read in
//
//
// Time periods for the summarization
// obs_var (string array) is added and works like grib_code (int array)
// when the obs name is given instead of grib_code
//
${METPLUS_TIME_SUMMARY_DICT}
//
// Mapping of input little_r report types to output message types
//
message_type_map = [
{ key = "FM-12 SYNOP"; val = "ADPSFC"; },
{ key = "FM-13 SHIP"; val = "SFCSHP"; },
{ key = "FM-15 METAR"; val = "ADPSFC"; },
{ key = "FM-18 BUOY"; val = "SFCSHP"; },
{ key = "FM-281 QSCAT"; val = "ASCATW"; },
{ key = "FM-32 PILOT"; val = "ADPUPA"; },
{ key = "FM-35 TEMP"; val = "ADPUPA"; },
{ key = "FM-88 SATOB"; val = "SATWND"; },
{ key = "FM-97 ACARS"; val = "AIRCFT"; }
];
//
// Indicate a version number for the contents of this configuration file.
// The value should generally not be modified.
//
//version = "V10.0";
tmp_dir = "${MET_TMP_DIR}";
${METPLUS_MET_CONFIG_OVERRIDES}
STATAnalysisConfig_wrapped
////////////////////////////////////////////////////////////////////////////////
//
// STAT-Analysis configuration file.
//
// For additional information, see the MET_BASE/config/README file.
//
////////////////////////////////////////////////////////////////////////////////
//
// Filtering input STAT lines by the contents of each column
//
//model = [
${METPLUS_MODEL}
//desc = [
${METPLUS_DESC}
//fcst_lead = [
${METPLUS_FCST_LEAD}
//obs_lead = [
${METPLUS_OBS_LEAD}
//fcst_valid_beg =
${METPLUS_FCST_VALID_BEG}
//fcst_valid_end =
${METPLUS_FCST_VALID_END}
fcst_valid_inc = [];
fcst_valid_exc = [];
//fcst_valid_hour = [
${METPLUS_FCST_VALID_HOUR}
//obs_valid_beg =
${METPLUS_OBS_VALID_BEG}
//obs_valid_end =
${METPLUS_OBS_VALID_END}
obs_valid_inc = [];
obs_valid_exc = [];
//obs_valid_hour = [
${METPLUS_OBS_VALID_HOUR}
//fcst_init_beg =
${METPLUS_FCST_INIT_BEG}
//fcst_init_end =
${METPLUS_FCST_INIT_END}
fcst_init_inc = [];
fcst_init_exc = [];
//fcst_init_hour = [
${METPLUS_FCST_INIT_HOUR}
//obs_init_beg =
${METPLUS_OBS_INIT_BEG}
//obs_init_end =
${METPLUS_OBS_INIT_END}
obs_init_inc = [];
obs_init_exc = [];
//obs_init_hour = [
${METPLUS_OBS_INIT_HOUR}
//fcst_var = [
${METPLUS_FCST_VAR}
//obs_var = [
${METPLUS_OBS_VAR}
//fcst_units = [
${METPLUS_FCST_UNITS}
//obs_units = [
${METPLUS_OBS_UNITS}
//fcst_lev = [
${METPLUS_FCST_LEVEL}
//obs_lev = [
${METPLUS_OBS_LEVEL}
//obtype = [
${METPLUS_OBTYPE}
//vx_mask = [
${METPLUS_VX_MASK}
//interp_mthd = [
${METPLUS_INTERP_MTHD}
//interp_pnts = [
${METPLUS_INTERP_PNTS}
//fcst_thresh = [
${METPLUS_FCST_THRESH}
//obs_thresh = [
${METPLUS_OBS_THRESH}
//cov_thresh = [
${METPLUS_COV_THRESH}
//alpha = [
${METPLUS_ALPHA}
//line_type = [
${METPLUS_LINE_TYPE}
//column = [
${METPLUS_COLUMN}
//weight = [
${METPLUS_WEIGHT}
////////////////////////////////////////////////////////////////////////////////
//
// Array of STAT-Analysis jobs to be performed on the filtered data
//
//jobs = [
${METPLUS_JOBS}
////////////////////////////////////////////////////////////////////////////////
//
// Confidence interval settings
//
out_alpha = 0.05;
boot = {
interval = PCTILE;
rep_prop = 1.0;
n_rep = 0;
rng = "mt19937";
seed = "";
}
////////////////////////////////////////////////////////////////////////////////
//
// WMO mean computation logic
//
wmo_sqrt_stats = [ "CNT:FSTDEV", "CNT:OSTDEV", "CNT:ESTDEV",
"CNT:RMSE", "CNT:RMSFA", "CNT:RMSOA",
"VCNT:FS_RMS", "VCNT:OS_RMS", "VCNT:RMSVE",
"VCNT:FSTDEV", "VCNT:OSTDEV" ];
wmo_fisher_stats = [ "CNT:PR_CORR", "CNT:SP_CORR",
"CNT:KT_CORR", "CNT:ANOM_CORR" ];
////////////////////////////////////////////////////////////////////////////////
//
// Skill score index options
//
//ss_index_name =
${METPLUS_SS_INDEX_NAME}
//ss_index_vld_thresh =
${METPLUS_SS_INDEX_VLD_THRESH}
////////////////////////////////////////////////////////////////////////////////
//hss_ec_value =
${METPLUS_HSS_EC_VALUE}
rank_corr_flag = FALSE;
vif_flag = FALSE;
tmp_dir = "${MET_TMP_DIR}";
//version = "V10.0";
${METPLUS_MET_CONFIG_OVERRIDES}
Python Embedding
This use case does not use Python embedding.
User Scripting
This use case uses two Python UserScripts to perform plotting of the output. These scripts are called via the METplus configuration file, and produce a map plot of the mean error (ME, bias) at each observation location from PointStat for both variables analyzed, as well as a 2D histogram plotting the frequency of each forecast minus observation and forecast value combination.
parm/use_cases/model_applications/land_surface/PointStat_fcstUFS_obsISMN_SoilMoistureTemp/plot_point_stat_bias_map_ISMN.py
# Script to read Point-Stat MPR files and Stat-Analysis -out_stat files and plot
# ME (bias) at station locations on a map
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import glob
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import numpy as np
import os
import pandas as pd
import sys
stat_analysis_file = sys.argv[1]
point_stat_output_dir = sys.argv[2]
obs_var = sys.argv[3]
obs_lev = sys.argv[4]
out_dir = sys.argv[5]
# Read the stat_analysis_file
sa_df = pd.read_csv(stat_analysis_file,sep='\\s+')
# Subset the stat_analysis file to only those lines with the requested variable
sa_df = sa_df[sa_df['OBS_VAR'].astype('str')==obs_var].reset_index()
# Sort the stat_analysis file by VX_MASK (OBS_SID)
sa_df = sa_df.sort_values(by='VX_MASK').reset_index()
# Read all available _mpr.txt files from point_stat at the point_stat_output_dir
# The VX_MASK column in the stat_analysis_file contains the OBS_SID values where there
# are bias values. In order to get the OBS_LAT and OBS_LON values for plotting, we
# need to collect all stations from all MPR files from PointStat. This is because if we pick
# only a single MPR file from a single time, all stations may not be present. So we must
# create a superset from all MPR files to ensure we get all the OBS_LAT/OBS_LON values
# for all stations in the stat_analysis output
point_stat_mpr_files = glob.glob(os.path.join(point_stat_output_dir,'*_mpr.txt'))
point_stat_mpr_files.sort()
df_list = [pd.read_csv(sf,sep='\\s+') for sf in point_stat_mpr_files]
sid_list = [df[['OBS_SID','OBS_LAT','OBS_LON','OBS_VAR']] for df in df_list]
all_sid = pd.concat(sid_list)
all_sid = all_sid[all_sid['OBS_VAR'].astype('str')==obs_var]
unique_sid = all_sid.drop_duplicates(subset=['OBS_SID']).sort_values(by='OBS_SID').reset_index()
# Ensure the list of stations from stat_analysis matches the unique list of stations from all MPR files
# and that they are in the same order
if (unique_sid['OBS_SID'].astype('str')==sa_df['VX_MASK'].astype('str')).all():
sa_df['OBS_LAT'] = unique_sid['OBS_LAT']
sa_df['OBS_LON'] = unique_sid['OBS_LON']
else:
print("ERROR!")
print("LIST OF STATIONS FROM MPR FILES DOES NOT MATCH STAT_ANALYSIS OUTPUT.")
print("Error in UserScript: plot_point_stat_bias_map_ISMN.py")
exit()
# Obtain some metadata to adorn the plot
fcst_valid_beg = sa_df['FCST_VALID_BEG'].astype('str').iloc[0]
fcst_valid_end = sa_df['FCST_VALID_END'].astype('str').iloc[0]
var_level = sa_df['FCST_LEV'].astype('str').iloc[0]
title_string_left = f"F-O (Mean Error, ME) for all forecasts valid \nfrom: {fcst_valid_beg}\nto: {fcst_valid_end}"
title_string_right = f"VAR_NAME: {obs_var}\nVAR_LEVEL: {var_level}"
# Now the sa_df dataframe contains OBS_LAT, OBS_LON, VX_MASK (OBS_SID), and ME, which we can plot on a map
# Create a figure with a map
fig = plt.figure(1,figsize=(22,15))
proj = ccrs.LambertConformal(central_longitude=-97.5,central_latitude=38.5)
ax1 = plt.subplot(111,projection=proj)
ax1.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5)
ax1.add_feature(cfeature.STATES, linewidth=0.5)
ax1.add_feature(cfeature.BORDERS, linewidth=0.5)
ax1.set_extent([-125,-65,20,55])
ax1.set_title(title_string_left,loc='left',fontsize=18)
ax1.set_title(title_string_right,loc='right',fontsize=18)
# Create a color normalization to use for coloring the dots based on bias values
if obs_var=='SOILW':
levels = np.arange(-0.1,0.11,0.01)
units = 'fraction'
elif obs_var=='TSOIL':
levels = np.arange(-10,11,1)
units = 'K'
cmap = plt.get_cmap('RdBu_r')
norm = mcolors.BoundaryNorm(levels, ncolors=cmap.N)
scatter = ax1.scatter(sa_df['OBS_LON'],sa_df['OBS_LAT'],c=sa_df['ME'],marker='.',s=500,cmap=cmap,norm=norm,transform=ccrs.PlateCarree())
cb = plt.colorbar(scatter,ax=ax1,orientation='vertical',extend='both',ticks=levels,shrink=0.75,pad=0.01)
cb.set_label(units,fontsize=18)
cb.ax.tick_params(labelsize=18)
plt.savefig(f"{out_dir}/{fcst_valid_beg}_{fcst_valid_end}_{obs_var}_{obs_lev}_ME_map.png")
parm/use_cases/model_applications/land_surface/PointStat_fcstUFS_obsISMN_SoilMoistureTemp/plot_point_stat_bias_2d_hist_ISMN.py
# Script to read Point-Stat MPR files and Stat-Analysis -out_stat files and plot
# ME (bias) at station locations on a map
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import glob
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import numpy as np
import os
import pandas as pd
import sys
stat_analysis_file = sys.argv[1]
point_stat_output_dir = sys.argv[2]
obs_var = sys.argv[3]
obs_lev = sys.argv[4]
out_dir = sys.argv[5]
# Read the stat_analysis_file
sa_df = pd.read_csv(stat_analysis_file,sep='\\s+')
# Subset the stat_analysis file to only those lines with the requested variable
sa_df = sa_df[sa_df['OBS_VAR'].astype('str')==obs_var].reset_index()
# Sort the stat_analysis file by VX_MASK (OBS_SID)
sa_df = sa_df.sort_values(by='VX_MASK').reset_index()
# Read all available _mpr.txt files from point_stat at the point_stat_output_dir
# The VX_MASK column in the stat_analysis_file contains the OBS_SID values where there
# are bias values. In order to get the OBS_LAT and OBS_LON values for plotting, we
# need to collect all stations from all MPR files from PointStat. This is because if we pick
# only a single MPR file from a single time, all stations may not be present. So we must
# create a superset from all MPR files to ensure we get all the OBS_LAT/OBS_LON values
# for all stations in the stat_analysis output
point_stat_mpr_files = glob.glob(os.path.join(point_stat_output_dir,'*_mpr.txt'))
point_stat_mpr_files.sort()
df_list = [pd.read_csv(sf,sep='\\s+') for sf in point_stat_mpr_files]
sid_list = [df[['OBS_SID','OBS_LAT','OBS_LON','OBS_VAR']] for df in df_list]
all_sid = pd.concat(sid_list)
all_sid = all_sid[all_sid['OBS_VAR'].astype('str')==obs_var]
unique_sid = all_sid.drop_duplicates(subset=['OBS_SID']).sort_values(by='OBS_SID').reset_index()
# Ensure the list of stations from stat_analysis matches the unique list of stations from all MPR files
# and that they are in the same order
if (unique_sid['OBS_SID'].astype('str')==sa_df['VX_MASK'].astype('str')).all():
sa_df['OBS_LAT'] = unique_sid['OBS_LAT']
sa_df['OBS_LON'] = unique_sid['OBS_LON']
else:
print("ERROR!")
print("LIST OF STATIONS FROM MPR FILES DOES NOT MATCH STAT_ANALYSIS OUTPUT.")
print("Error in UserScript: plot_point_stat_bias_map_ISMN.py")
exit()
# Obtain some metadata to adorn the plot
fcst_valid_beg = sa_df['FCST_VALID_BEG'].astype('str').iloc[0]
fcst_valid_end = sa_df['FCST_VALID_END'].astype('str').iloc[0]
var_level = sa_df['FCST_LEV'].astype('str').iloc[0]
title_string_left = f"F-O (Mean Error, ME) for all forecasts valid \nfrom: {fcst_valid_beg}\nto: {fcst_valid_end}"
title_string_right = f"VAR_NAME: {obs_var}\nVAR_LEVEL: {var_level}"
# Now the sa_df dataframe contains OBS_LAT, OBS_LON, VX_MASK (OBS_SID), and ME, which we can plot on a map
# Create a figure with a map
fig = plt.figure(1,figsize=(22,15))
proj = ccrs.LambertConformal(central_longitude=-97.5,central_latitude=38.5)
ax1 = plt.subplot(111,projection=proj)
ax1.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5)
ax1.add_feature(cfeature.STATES, linewidth=0.5)
ax1.add_feature(cfeature.BORDERS, linewidth=0.5)
ax1.set_extent([-125,-65,20,55])
ax1.set_title(title_string_left,loc='left',fontsize=18)
ax1.set_title(title_string_right,loc='right',fontsize=18)
# Create a color normalization to use for coloring the dots based on bias values
if obs_var=='SOILW':
levels = np.arange(-0.1,0.11,0.01)
units = 'fraction'
elif obs_var=='TSOIL':
levels = np.arange(-10,11,1)
units = 'K'
cmap = plt.get_cmap('RdBu_r')
norm = mcolors.BoundaryNorm(levels, ncolors=cmap.N)
scatter = ax1.scatter(sa_df['OBS_LON'],sa_df['OBS_LAT'],c=sa_df['ME'],marker='.',s=500,cmap=cmap,norm=norm,transform=ccrs.PlateCarree())
cb = plt.colorbar(scatter,ax=ax1,orientation='vertical',extend='both',ticks=levels,shrink=0.75,pad=0.01)
cb.set_label(units,fontsize=18)
cb.ax.tick_params(labelsize=18)
plt.savefig(f"{out_dir}/{fcst_valid_beg}_{fcst_valid_end}_{obs_var}_{obs_lev}_ME_map.png")
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/land_surface/PointStat_fcstUFS_obsISMN_SoilMoistureTemp.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}/model_applications/land_surface/PointStat_fcstUFS_obsISMN_SoilMoistureTemp and will contain the following files:
* ascii2nc/2020080512_ismn.nc
* point_stat/point_stat_600000L_20200805_120000V_mpr.txt
* point_stat/point_stat_600000L_20200805_120000V.stat
* stat_analysis/2020080512_2020080512_OBS_SID_CNT.stat
* 20200805_120000_20200805_120000_SOILW_0-0.1_2D_hist.png
* 20200805_120000_20200805_120000_SOILW_0-0.1_ME_map.png
* 20200805_120000_20200805_120000_TSOIL_0-0.1_2D_hist.png
* 20200805_120000_20200805_120000_TSOIL_0-0.1_ME_map.png
Keywords
Note
ASCII2NCToolUseCase
PointStatToolUseCase
StatAnalysisToolUseCase
UserScriptUseCase
LandSurfaceAppUseCase
VxDataISMN
Navigate to the METplus Quick Search for Use Cases page to discover other similar use cases.
sphinx_gallery_thumbnail_path = ‘_static/land_surface-PointStat_fcstUFS_obsISMN_SoilMoistureTemp.png’