Note
Go to the end to download the full example code
PointStat: read in satellite data and verify wind speeds or wave heights
model_applications/marine_and_cryosphere/PointStat_fcstGFS_obsJASON3_satelliteAltimetry.conf
Scientific Objective
Satellite data provides a wealth of information, especially over vast water bodies (eg. oceans) where traditional observation methods are sparse or unavailable. This use case shows how a satellite dataset can be used as observations to verify against a model forecast. While the use case is set up to verify using JASON-3 data, the Python script called on via Python Embedding is capabile of processing SARAL and Sentinel-6a datasets as well.
Datasets
METplus Components
This use case calls Python Embedding during PointStat, which is the only tool used.
METplus Workflow
PointStat kicks off a Python script execution, which reads in the file name, variable field of interest, and type of file (JASON, SARAL, or SENTINEL). After these points are passed back to PointStat as the point observation dataset, they are compared to gridded forecast data. CTC and CTS line types are output, which can be adjusted for additional wind speeds/ wave heights. The use case processes the following run 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/marine_and_cryosphere/PointStat_fcstGFS_obsJASON3_satelliteAltimetry.conf
[config]
# Documentation for this use case can be found at
# https://metplus.readthedocs.io/en/latest/generated/model_applications/marine_and_cryosphere/PointStat_fcstGFS_obsJASON3_satelliteAltimetry.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, PointStat(wind)
###
# 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_TIME_FMT = %Y%m%d_%H
VALID_BEG = 20240102_12
VALID_END = 20240102_12
VALID_INCREMENT = 1d
LEAD_SEQ = 12
###
# File I/O
# https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#directory-and-filename-template-info
###
CONFIG_DIR = {PARM_BASE}/use_cases/model_applications/marine_and_cryosphere/PointStat_fcstGFS_obsJASON3_satelliteAltimetry
FCST_POINT_STAT_INPUT_DIR = {INPUT_BASE}/model_applications/marine_and_cryosphere/PointStat_fcstGFS_obsJASON3_satelliteAltimetry
FCST_POINT_STAT_INPUT_TEMPLATE = gfswave.t00z.global.0p25.f{lead?fmt=%3H}.grib2
OBS_POINT_STAT_INPUT_DIR =
OBS_POINT_STAT_INPUT_TEMPLATE = PYTHON_NUMPY= {CONFIG_DIR}/read_satData.py {INPUT_BASE}/model_applications/marine_and_cryosphere/PointStat_fcstGFS_obsJASON3_satelliteAltimetry/JA3_GPSOPR_2PfS362_215_{valid?fmt=%Y%m%d}_102102_{valid?fmt=%Y%m%d}_121958.nc:swh_ocean:JASON
POINT_STAT_OUTPUT_DIR = {OUTPUT_BASE}/model_applications/marine_and_cryosphere/PointStat_fcstGFS_obsJASON3_satelliteAltimetry
POINT_STAT_OUTPUT_PREFIX = swh
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 = WVHGT
FCST_VAR1_LEVELS = Z0
OBS_VAR1_NAME = swh_ocean
OBS_VAR1_LEVELS = Z0
###
# 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_CLIMO_MEAN_TIME_INTERP_METHOD = NEAREST
POINT_STAT_INTERP_TYPE_METHOD = BILIN
POINT_STAT_INTERP_TYPE_WIDTH = 2
POINT_STAT_OUTPUT_FLAG_CNT = STAT
POINT_STAT_OUTPUT_FLAG_SL1L2 =
OBS_POINT_STAT_WINDOW_BEGIN = -1800
OBS_POINT_STAT_WINDOW_END = 1800
POINT_STAT_OFFSETS = 0
MODEL = GFS
POINT_STAT_DESC = NA
OBTYPE =
POINT_STAT_REGRID_TO_GRID = NONE
POINT_STAT_REGRID_METHOD = BILIN
POINT_STAT_REGRID_WIDTH = 2
POINT_STAT_MESSAGE_TYPE = WDSATR
[wind]
OBS_POINT_STAT_INPUT_TEMPLATE = PYTHON_NUMPY= {CONFIG_DIR}/read_satData.py {INPUT_BASE}/model_applications/marine_and_cryosphere/PointStat_fcstGFS_obsJASON3_satelliteAltimetry/JA3_GPSOPR_2PfS362_215_{valid?fmt=%Y%m%d}_102102_{valid?fmt=%Y%m%d}_121958.nc:wind_speed_alt:JASON
FCST_VAR1_NAME = WIND
FCST_VAR1_LEVELS = Z0
OBS_VAR1_NAME = wind_speed_alt
OBS_VAR1_LEVELS = Z0
POINT_STAT_OUTPUT_PREFIX = wind
POINT_STAT_OUTPUT_FLAG_CNT =
POINT_STAT_OUTPUT_FLAG_SL1L2 = STAT
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
Note
See the PointStat MET Configuration section of the User’s Guide for more information on the environment variables used in the file below:
////////////////////////////////////////////////////////////////////////////////
//
// 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}
}
obs = {
${METPLUS_OBS_FILE_TYPE}
//field = [
${METPLUS_OBS_FIELD}
}
////////////////////////////////////////////////////////////////////////////////
//
// 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}
////////////////////////////////////////////////////////////////////////////////
//
// 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}
////////////////////////////////////////////////////////////////////////////////
tmp_dir = "${MET_TMP_DIR}";
// output_prefix =
${METPLUS_OUTPUT_PREFIX}
//version = "V10.0.0";
////////////////////////////////////////////////////////////////////////////////
${METPLUS_MET_CONFIG_OVERRIDES}
Python Embedding
This use case calls the read_satData.py script to read and pass to PointStat the user-requested variable. The script needs 3 inputs in the following order: an input file, a variable field to extract, and where the data came from, passed as JASON (JASON-3), SARAL, or SENTINEL (Sentinel-6a). The location of the code is parm/use_cases/model_applications/marine_and_cryosphere/PointStat_fcstGFS_obsJASON3_satelliteAltimetry/read_satData.py
#This script is designed to ingest one of three types of satellite netCDF files: JASON-3, SARAL, and Sentinel-6a.
#It also grabs the lat, lon, and time information and puts it into a list of lists, the accepted format of MET.
#Currently the script can accept any variable in the netCDFs; however, it's important to note that
#JASON-3 and Sentinel-6a data use groups and non-groups for data organization, which could cause an issue if the requested
#variable is outside of the hardcoded group. The original intent of the use case was for wind speeds and significant wave heights.
#FOR FUTURE REFERENCE
#JASON-3 is swh_ocean, SARAL is swh, and Sentinel-6a is swh_ocean
from netCDF4 import Dataset
import sys
import numpy as np
import datetime as dt
import xarray as xr
import pandas as pd
#from met.point import convert_point_data
#Users are responsible for passing the following arguements at runtime:
##input file
##Varible field to read in
##type of file (JASON, SARAL, or SENTINEL)
if len(sys.argv[1].split(':')) == 3:
try:
input_file,field_name,file_type = sys.argv[1].split(':')
except:
print("input directory may not be set correctly, dates may not be in correct format. Please recheck")
sys.exit()
if file_type == 'JASON' or file_type == 'SENTINEL':
#need to check if the variable is in the data_01 group or data_01/ku group
try:
ds = xr.open_dataset(input_file, group="/data_01/ku")
du = xr.open_dataset(input_file, group="/data_01")
obs_hold = ds[field_name]
except KeyError:
ds = xr.open_dataset(input_file, group="/data_01")
du = xr.open_dataset(input_file, group="/data_01")
obs_hold = ds[field_name]
obs = obs_hold.values
latitude = np.array(du.latitude.values)
longitude = np.array(du.longitude.values)
time = np.array(du.time.values)
#convert times to MET readable
new_time = []
for i in range(len(time)):
new_time.append(pd.to_datetime(str(time[i])).strftime("%Y%m%d_%H%M%S"))
elif file_type == 'SARAL':
f_in = Dataset(input_file,'r')
obs = np.array(f_in[field_name][:])
latitude = np.array(f_in['lat'][:])
longitude = np.array(f_in['lon'][:])
time = np.array(f_in['time'][:])
#adjust times to MET time
new_time = []
for i in range(len(time)):
new_time.append(dt.datetime(2000,1,1) + dt.timedelta(seconds=int(time[i])))
new_time[i] = new_time[i].strftime("%Y%m%d_%H%M%S")
else:
print('file type '+file_type+' not supported. Please use JASON, SARAL, or SENTINEL')
sys.exit()
#Currently, all station IDs are assigned the same value. If this is not the desired behavior it will require
#additional coding
sid = np.full(len(latitude),"1")
#get arrays into lists, then create a list of lists
#this also requires creating the last arrays of typ, elv, var, lvl, hgt, and qc
typ = np.full(len(latitude), 'WDSATR').tolist()
elv = np.zeros(len(latitude),dtype=int).tolist()
var = np.full(len(latitude), field_name).tolist()
lvl = np.full(len(latitude),1013.25).tolist()
#adding additional check; if 'wind' appears in the variable name, it's assumed
#to be a wind speed and gets a height of 10m; otherwise its a height of 0
if field_name.rfind('wind') != -1:
hgt = np.full(len(latitude),0,dtype=int).tolist()
else:
hgt = np.full(len(latitude),0,dtype=int).tolist()
qc = np.full(len(latitude),'NA').tolist()
#if the data has nans, this replaces them with -9999 (default bad value in MET)
#obs = np.nan_to_num(obs, nan=-9999)
sid = sid.tolist()
vld = new_time
lat = latitude.tolist()
lon = longitude.tolist()
obs = obs.tolist()
l_tuple = list(zip(typ,sid,vld,lat,lon,elv,var,lvl,hgt,qc,obs))
point_data = [list(ele) for ele in l_tuple]
#met_point_data = convert_point_data(point_data)
print("Data Length:\t" + repr(len(point_data)))
print("Data Type:\t" + repr(type(point_data)))
#if the incorrect number of args are passed, the system will print out the usage statement and end
else:
print("Run Command:\n\n read_satData.py /path/to/input/input_file:variable_field_name:file_type\n\nCommands notes:\nIf only certain variable fields return errors, the field may not be supported or outisde of the expected netCDF group. Additional coding changes may be required.\nfile names currently supported: JASON SARAL SENTINEL\nCurrent Message_type is hard-coded to WDSATR\n")
sys.exit()
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/marine_and_cryosphere/PointStat_fcstGFS_obsJASON3_satelliteAltimetry.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 model_applications/marine_and_cryosphere/PointStat_fcstGFS_obsJASON3_satelliteAltimetry (relative to OUTPUT_BASE) and will contain the following files:
point_stat_swh_120000L_20240102_120000V.stat
point_stat_wind_120000L_20240102_120000V.stat
Keywords
Note
PointStatToolUseCase
PythonEmbeddingFileUseCase
GRIB2FileUseCase
MarineAndCryosphereAppUseCase
Navigate to the METplus Quick Search for Use Cases page to discover other similar use cases.
sphinx_gallery_thumbnail_path = ‘_static/marine_and_cryosphere-PointStat_fcstGFS_obsJASON3_satelliteAltimetry.png’
Total running time of the script: (0 minutes 0.000 seconds)