Note
Go to the end to download the full example code
PointStat: CESM and FLUXNET2015 Terrestrial Coupling Index (TCI)
model_applications/land_surface/PointStat_fcstCESM_obsFLUXNET2015_TCI.conf
Scientific Objective
This use case ingests two CESM (CAM and CLM) files and a new FLUX2015 dataset (NETCDF) that Paul Dirmeyer, GMU prepared from FLUXNET2015 ASCII dataset. The use case will calculate Terrestrial Coupling Index (TCI) from CESM datasets. Utilizing Python embedding, this use case taps into a new vital observation dataset and compares it to CESM simulations TCI forecast. Finally, it will generate plots of model TCI and observed TCI.
Datasets
METplus Components
This use case utilizes the METplus PyEmbedIngest to read the CESM files and calculate TCI using python embedding and a NETCDF file of the TCI is generated. The METplus PointStat processes the output of PyEmbedIngest and FLUXNET2015 dataset (using python embedding), and outputs the requested line types. Then METplus PlotPointObs tool reads the output of PyEmbedIngest and FLUXNET2015 dataset and produce plots of TCI from CESM and point observations. A custom loop runs through all the pre-defined seasons (DJF, MAM, JJA, SON) and runs PyEmbedIngest, PointStat, and PlotPointObs.
METplus Workflow
The PyEmbedIngest tool reads 2 CESM files containing Soil Moisture and Sensible Heat Flux, each composed of daily forecasts from 1979 to 1983 and calculates TCI and generates a NETCDF file of the TCI. One FLUXNET2015 NETCDF file containing station observations of several variables including Coupling Index of Soil Moisture and Sensible Heat Flux is read by Python Embedding.
PointStat is used to compare the two new fields (TCI calculated from CESM dataset and FLUXNET2015). Finally, PlotPointObs is run to plot the CESM TCI overlaying the FLUXNET2015 point observations.
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. -c parm/use_cases/model_applications/land_surface/PointStat_fcstCESM_obsFLUXNET2015_TCI.conf
[config]
# Documentation for this use case can be found at
# https://metplus.readthedocs.io/en/latest/generated/model_applications/land_surface/PointStat_fcstCESM_obsFLUXNET2015_TCI.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 = PyEmbedIngest,PointStat,PlotPointObs
############################################################
# 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 = 1979060100
VALID_END = 1979060100
VALID_INCREMENT = 24H
LEAD_SEQ = 0
############################################################
# Pre-determined seasons for the use case
############################################################
CUSTOM_LOOP_LIST = DJF,MAM,JJA,SON
LOG_LEVEL=DEBUG
############################################################
# PyEmbedIngest Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#pyembedingest
############################################################
PY_EMBED_INGEST_1_OUTPUT_DIR =
PY_EMBED_INGEST_1_OUTPUT_TEMPLATE = {OUTPUT_BASE}/regrid_data_plane_{custom}.nc
############################################################
# The last argument is the Latent Heat Flux.
# User can change it to use any other variable present in CESM files.
############################################################
PY_EMBED_INGEST_1_SCRIPT = {PARM_BASE}/use_cases/model_applications/land_surface/PointStat_fcstCESM_obsFLUXNET2015_TCI/cesm_tci.py {INPUT_BASE}/model_applications/land_surface/PointStat_fcstCESM_obsFLUXNET2015_TCI/f.e21.FHIST.f09_f09_mg17.CESM2-CLM45physics.002.clm2.h1.1979-83_SoilWater10cm.nc {INPUT_BASE}/model_applications/land_surface/PointStat_fcstCESM_obsFLUXNET2015_TCI/f.e21.FHIST.f09_f09_mg17.CESM2-CLM45physics.002.cam.h1.1979-83_CIvars.nc {custom} LHFLX
PY_EMBED_INGEST_1_TYPE = NUMPY
#PY_EMBED_INGEST_1_OUTPUT_GRID = G129
PY_EMBED_INGEST_1_OUTPUT_GRID = "latlon 288 192 -90.0 -0.0 0.9375 1.25"
############################################################
# PointStat Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#pointstat
############################################################
FCST_POINT_STAT_INPUT_TEMPLATE = PYTHON_NUMPY
############################################################
# The default variables for the FLUXNET2015 observations id 10CM Soil Moisture
# and Latent Heat Flux. User can change them in the python Embedding script
# fluxnet_tci.py in the PointStat_fcstCESM_obsFLUXNET2015_TCI
# directory.
############################################################
OBS_POINT_STAT_INPUT_TEMPLATE = PYTHON_NUMPY={PARM_BASE}/use_cases/model_applications/land_surface/PointStat_fcstCESM_obsFLUXNET2015_TCI/fluxnet2015_tci.py {INPUT_BASE}/model_applications/land_surface/PointStat_fcstCESM_obsFLUXNET2015_TCI/F2015_LoCo_AllChains_F2015.nc4 {custom}
POINT_STAT_OUTPUT_DIR = {OUTPUT_BASE}/PointStat
############################################################
# Field Info
# https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#field-info
############################################################
POINT_STAT_ONCE_PER_FIELD = False
FCST_POINT_STAT_INPUT_DIR = {OUTPUT_BASE}/
FCST_POINT_STAT_INPUT_TEMPLATE = regrid_data_plane_{custom}.nc
FCST_POINT_STAT_VAR1_LEVELS =
OBS_POINT_STAT_VAR1_NAME = TCI
OBS_POINT_STAT_VAR1_LEVELS = L0
FCST_POINT_STAT_VAR1_NAME = TCI_10cm_soil_depth
FCST_POINT_STAT_VAR1_LEVELS = Z10
BOTH_POINT_STAT_VAR1_THRESH =
############################################################
# PointStat Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#pointstat
############################################################
LOG_POINT_STAT_VERBOSITY = 2
POINT_STAT_CONFIG_FILE ={PARM_BASE}/met_config/PointStatConfig_wrapped
POINT_STAT_INTERP_TYPE_METHOD = NEAREST
POINT_STAT_INTERP_TYPE_WIDTH = 1
POINT_STAT_OUTPUT_FLAG_CTC = BOTH
POINT_STAT_OUTPUT_FLAG_MPR = BOTH
POINT_STAT_OUTPUT_FLAG_CNT = BOTH
OBS_POINT_STAT_WINDOW_BEGIN = -1000000000
OBS_POINT_STAT_WINDOW_END = 2000000000
POINT_STAT_OFFSETS = 0
MODEL = CESM
POINT_STAT_DESC = TCI
OBTYPE =
POINT_STAT_REGRID_TO_GRID = NONE
POINT_STAT_REGRID_METHOD = BILIN
POINT_STAT_REGRID_WIDTH = 2
POINT_STAT_OBS_VALID_BEG = 19790101_000000
POINT_STAT_OBS_VALID_END = 20130101_000000
POINT_STAT_MASK_GRID = FULL
POINT_STAT_MASK_POLY =
POINT_STAT_MASK_SID =
POINT_STAT_MESSAGE_TYPE = ADPSFC
POINT_STAT_OUTPUT_PREFIX = {custom}
############################################################
# PlotPointObs Settings
# https://metplus.readthedocs.io/en/latest/Users_Guide/wrappers.html#plotpointobs
############################################################
PLOT_POINT_OBS_INPUT_TEMPLATE =PYTHON_NUMPY={PARM_BASE}/use_cases/model_applications/land_surface/PointStat_fcstCESM_obsFLUXNET2015_TCI/fluxnet2015_tci.py {INPUT_BASE}/model_applications/land_surface/PointStat_fcstCESM_obsFLUXNET2015_TCI/F2015_LoCo_AllChains_F2015.nc4 {custom}
PLOT_POINT_OBS_GRID_INPUT_DIR = {OUTPUT_BASE}/
PLOT_POINT_OBS_GRID_INPUT_TEMPLATE = regrid_data_plane_{custom}.nc
PLOT_POINT_OBS_OUTPUT_DIR = {OUTPUT_BASE}/plot_point_obs
PLOT_POINT_OBS_OUTPUT_TEMPLATE = cesm_fluxnet2015_{custom}.ps
PLOT_POINT_OBS_TITLE = {custom} CESM vs FLUXNET2015 TCI
LOG_PLOT_POINT_OBS_VERBOSITY = 2
PLOT_POINT_OBS_GRID_DATA_FIELD = { name = "TCI_10cm_soil_depth"; level = "10cm_soil_depth"; }
PLOT_POINT_OBS_GRID_DATA_GRID_PLOT_INFO_PLOT_MIN = -40.0
PLOT_POINT_OBS_GRID_DATA_GRID_PLOT_INFO_PLOT_MAX = 55.0
PLOT_POINT_OBS_POINT_DATA =
{
msg_typ = "ADPSFC";
obs_thresh = >-9999;
obs_var = "TCI";
dotsize(x) = 3.5;
line_color = [0, 0, 0];
fill_plot_info = {
flag = TRUE;
color_table = "MET_BASE/colortables/met_default.ctable";
plot_min = -40.0;
plot_max = 55.0;
colorbar_flag = FALSE;
}
}
MET Configuration
METplus sets environment variables based on the values in the METplus configuration file. These variables are referenced in the MET configuration file. 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 not controlled by an environment variable, you can add additional environment variables to be set only within the METplus environment using the [user_env_vars] section of the METplus configuration files. See the ‘User Defined Config’ section on the ‘System Configuration’ page of the METplus User’s Guide for more information.
////////////////////////////////////////////////////////////////////////////////
//
// 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}
${METPLUS_FCST_FIELD}
}
obs = {
${METPLUS_OBS_FILE_TYPE}
${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 = NONE;
obs_summary = NONE;
obs_perc_value = 50;
//
// 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}
////////////////////////////////////////////////////////////////////////////////
//
// 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}
////////////////////////////////////////////////////////////////////////////////
//
// Statistical output types
//
//output_flag = {
${METPLUS_OUTPUT_FLAG_DICT}
////////////////////////////////////////////////////////////////////////////////
// Threshold for SEEPS p1 (Probability of being dry)
//seeps_p1_thresh =
${METPLUS_SEEPS_P1_THRESH}
////////////////////////////////////////////////////////////////////////////////
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/land_surface/PointStat_fcstCESM_obsFLUXNET2015_TCI/cesm_tci.py
"""
Code adapted from Meg Fowler,CGD,NCAR
The code computes the Terrestrial Coupling Index (TCI)
from latent Heat Flux and 10 cm Soil Moisture
Designed to read Latent Heat Flux (from CAM) and Soil Temp (from CLM)
from two CESM files
User needs to provide the season (DJF, MAM, JJA, or SON) in the METplus conf file
User can change the variables to compute TCI
"""
import numpy as np
import xarray as xr
import pandas as pd
import datetime
import time
import sys
import os
import netCDF4 as nc
if len(sys.argv) < 4:
print("Must specify the following elements: sfc_flux_file soil_file season (DJF, MAM, JJA, or SON)")
sys.exit(1)
fileCLM = os.path.expandvars(sys.argv[1])
fileCAM = os.path.expandvars(sys.argv[2])
season = sys.argv[3]
var_y = sys.argv[4]
print('Starting Terrestrial Coupling Index Calculation for: ',season)
camDS_CLM45 = xr.open_dataset(fileCAM, decode_times=False)
print('Finished reading in CAM file')
clmDS_CLM45 = xr.open_dataset(fileCLM, decode_times=False)
print('Finished reading in CLM file')
if season=="DJF":
ss = 0
elif season=="MAM":
ss = 1
elif season=="JJA":
ss = 2
elif season=="SON":
ss = 3
else:
sys.exit('ERROR : URECOGNIZED SEASON, PLEASE USE DJF, MAM, JJA, OR SON')
units, reference_date = camDS_CLM45.time.attrs['units'].split('since')
camDS_CLM45['time'] = pd.date_range(start=reference_date, periods=camDS_CLM45.sizes['time'], freq='D')
camDS_time=camDS_CLM45.time[0]
dt_object = datetime.datetime.utcfromtimestamp(camDS_time.values.tolist()/1e9)
vDate=dt_object.strftime("%Y%m%d")
units, reference_date = clmDS_CLM45.time.attrs['units'].split('since')
clmDS_CLM45['time'] = pd.date_range(start=reference_date, periods=clmDS_CLM45.sizes['time'], freq='D')
ds = camDS_CLM45
ds['SOILWATER_10CM'] = (('time','lat','lon'), clmDS_CLM45.SOILWATER_10CM.values)
xname = 'SOILWATER_10CM' # Controlling variable
#yname = 'LHFLX' # Responding variable
yname=str(var_y)
xday = ds[xname].groupby('time.season')
yday = ds[yname].groupby('time.season')
# Get the covariance of the two (numerator in coupling index)
covarTerm = ((xday - xday.mean()) * (yday - yday.mean())).groupby('time.season').sum() / xday.count()
# Now compute the coupling index
couplingIndex = covarTerm/xday.std()
ci_season=couplingIndex[ss,:,:]
ci = np.where(np.isnan(ci_season), -9999, ci_season)
met_data = ci[:,:]
met_data = met_data[::-1].copy()
#trim the lat/lon grids so they match the data fields
lat_met = camDS_CLM45.lat
lon_met = camDS_CLM45.lon
print(" Model Data shape: "+repr(met_data.shape))
v_str = vDate
v_str = v_str + '_000000'
#print(" Valid date: "+v_str)
lat_ll = float(lat_met.min())
lon_ll = float(lon_met.min())
n_lat = lat_met.shape[0]
n_lon = lon_met.shape[0]
delta_lat = (float(lat_met.max()) - float(lat_met.min()))/float(n_lat)
delta_lon = (float(lon_met.max()) - float(lon_met.min()))/float(n_lon)
print(f"variables:"
f"lat_ll: {lat_ll} lon_ll: {lon_ll} n_lat: {n_lat} n_lon: {n_lon} delta_lat: {delta_lat} delta_lon: {delta_lon}")
attrs = {
'valid': v_str,
'init': v_str,
'lead': "000000",
'accum': "000000",
'name': 'TCI',
'standard_name': 'terrestrial_coupling_index',
'long_name': 'terrestrial_coupling_index',
'level': "10cm_soil_depth",
'units': "W/m2",
'grid': {
'type': "LatLon",
'name': "CESM Grid",
'lat_ll': lat_ll,
'lon_ll': lon_ll,
'delta_lat': delta_lat,
'delta_lon': delta_lon,
'Nlat': n_lat,
'Nlon': n_lon,
}
}
parm/use_cases/model_applications/land_surface/PointStat_fcstCESM_obsFLUXNET2015_TCI/fluxnet2015_tci.py
"""
This python embedding code reads in the FLUXNET2015 NETCDF data
and passed to point_data. User need to pass on the season.
User can change the verifying variable
"""
import numpy
import sys
import os
import glob
import numpy as np
import xarray as xr
import datetime
from datetime import date, timedelta
import pandas as pd
if len(sys.argv) < 2:
print("Must specify the following elements: FLUXNET2015_file season")
sys.exit(1)
obsfile = os.path.expandvars(sys.argv[1])
season = sys.argv[2]
f2015data = xr.open_dataset(obsfile, decode_times=False)
if season=="DJF":
ss = 0
elif season=="MAM":
ss = 1
elif season=="JJA":
ss = 2
elif season=="SON":
ss = 3
else:
print('Unrecognized season, please use DJF, MAM, JJA, SON')
exit()
start_year=f2015data['Start year'].values.tolist()
end_year=f2015data['End year'].values.tolist()
vld=pd.to_datetime(start_year,format='%Y')
vld=vld.strftime("%Y%m%d")
vld=vld+ '_000000'
sid=f2015data['Station'].values.tolist()
print("Length :",len(sid))
print("SID :",sid)
lat=f2015data['Latitude'].values.tolist()
lon=f2015data['Longitude'].values.tolist()
# User can change the name of the variable below
obs=f2015data['CI Sfc_SM Latent_Heat'].values.tolist()
obs=np.array(obs)
obs=obs[:,ss]
#create dummy lists for the message type, elevation, variable name, level, height, and qc string
#numpy is more efficient at creating the lists, but need to convert to Pythonic lists
typ = np.full(len(sid),'ADPSFC').tolist()
elv = np.full(len(sid),10).tolist()
var = np.full(len(sid),'TCI').tolist()
lvl = np.full(len(sid),10).tolist()
hgt = np.zeros(len(sid),dtype=int).tolist()
qc = np.full(len(sid),'NA').tolist()
obs = np.where(np.isnan(obs), -9999, obs)
obs = np.full(len(sid),obs).tolist()
vld = np.full(len(sid),vld).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]
#print("Data Length:\t" + repr(len(point_data)))
#print("Data Type:\t" + repr(type(point_data)))
#print(" Point Data Shape: ",np.shape(point_data))
print(" Point Data: ",point_data)
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_fcstCESM_obsFLUXNET2015_TCI.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 the use case will be found in 3 folders(relative to OUTPUT_BASE). Those folders are:
PyEmbedIngest
The OUTPUT_BASE folder contains all of the TCI output calculated using CESM files in NETCDF format:
regrid_data_plane_DJF.nc
regrid_data_plane_JJA.nc
regrid_data_plane_MAM.nc
regrid_data_plane_SON.nc
PointStat
The final folder, PointStat, contains all of the following output from the PointStat call:
point_stat_DJF_000000L_19790101_000000V_cnt.txt
point_stat_DJF_000000L_19790101_000000V_ctc.txt
point_stat_DJF_000000L_19790101_000000V_mpr.txt
point_stat_DJF_000000L_19790101_000000V.stat
point_stat_JJA_000000L_19790101_000000V_cnt.txt
point_stat_JJA_000000L_19790101_000000V_ctc.txt
point_stat_JJA_000000L_19790101_000000V_mpr.txt
point_stat_JJA_000000L_19790101_000000V.stat
point_stat_MAM_000000L_19790101_000000V_cnt.txt
point_stat_MAM_000000L_19790101_000000V_ctc.txt
point_stat_MAM_000000L_19790101_000000V_mpr.txt
point_stat_MAM_000000L_19790101_000000V.stat
point_stat_SON_000000L_19790101_000000V_cnt.txt
point_stat_SON_000000L_19790101_000000V_ctc.txt
point_stat_SON_000000L_19790101_000000V_mpr.txt
point_stat_SON_000000L_19790101_000000V.stat
PlotPointObs
The final folder plot_point_obs, contains all of the plots from the PlotPointObs call:
cesm_fluxnet2015_DJF.ps
cesm_fluxnet2015_JJA.ps
cesm_fluxnet2015_MAM.ps
cesm_fluxnet2015_SON.ps
Keywords
Note
PyEmbedIngestToolUseCase
PointStatToolUseCase
PlotPointObsToolUseCase
PythonEmbeddingFileUseCase
NETCDFFileUseCase
LandSurfaceAppUseCase
Navigate to the METplus Quick Search for Use Cases page to discover other similar use cases.
sphinx_gallery_thumbnail_path = ‘_static/land_surface-PointStat_fcstCESM_obsFLUXNET2015_TCI.png’
Total running time of the script: (0 minutes 0.000 seconds)