Enhance TC-Diag to actually compute and write diagnostics

[JohnHalleyGotway]

Describe the Enhancement

Review details of the initial development version in #2168.

The initial version of the TC-Diag tool is slated for inclusion in the MET-11.1.0 release. It transforms a configurable list of input gridded data fields into range/azimuth cylindrical coordinates and calls python script to compute diagnostics on the result. If the storm track (i.e. OFCL) differs from the gridded data (i.e. GFS), the center of the storm vortex and the track location may not coincide. This task is to incorporate NHC's vortex removal logic into TC-Diag prior to converting to cylindrical coordinates.

CSU/CIRA will provide the algorithm for doing so.

Time Estimate

Estimate the amount of work required here.
Issues should represent approximately 1 to 3 days of work.

Sub-Issues

Consider breaking the enhancement down into sub-issues.
None needed.

Relevant Deadlines

List relevant project deadlines here or state NONE.

Funding Source

Define the source of funding and account keys here or state NONE.

Define the Metadata

Assignee

• Select engineer(s) or no engineer required
• Select scientist(s) or no scientist required

Labels

• Select component(s)
• Select priority
• Select requestor(s)

Projects and Milestone

• Select Repository and/or Organization level Project(s) or add alert: NEED CYCLE ASSIGNMENT label
• Select Milestone as the next official version or Future Versions

Define Related Issue(s)

Consider the impact to the other METplus components.

• METplus, MET, METdataio, METviewer, METexpress, METcalcpy, METplotpy
  This may impact the METplus wrapper of TC-Diag.
• Enhance METbaseimage to install the YAML Python package METbaseimage#15 adds the new required yaml Python package to the base image testing environment.

Enhancement Checklist

See the METplus Workflow for details.

• Complete the issue definition above, including the Time Estimate and Funding Source.
• Fork this repository or create a branch of develop.
  Branch name: feature_<Issue Number>_<Description>
• Complete the development and test your changes.
• Add/update log messages for easier debugging.
• Add/update unit tests.
• Add/update documentation.
• Push local changes to GitHub.
• Submit a pull request to merge into develop.
  Pull request: feature <Issue Number> <Description>
• Define the pull request metadata, as permissions allow.
  Select: Reviewer(s) and Development issues
  Select: Repository level development cycle Project for the next official release
  Select: Milestone as the next official version
• Iterate until the reviewer(s) accept and merge your changes.
• Delete your fork or branch.
• Close this issue.

[JohnHalleyGotway]

As discussed at the TC-Diag project meeting on 9/12/23, we should ideally test with at least 3 input data sources:

1. NOAA generated GRIB2 inputs (with variable names like TMP, UGRD, and VGRD).
2. ECMWF generated GRIB2 inputs (with variable names like T, U, and V).
   ECMWF GRIB data are available at: mohawk:/d2/projects/TCGP/data/data_realtime/ecmwf/2023/09/12/06z
3. NetCDF wrfout files on which the GFDL vortex tracker can be run.

[JohnHalleyGotway]

From @robertdemariacira on 8/17/23:

I have a version of the code (v0.5.0) that can now read the provided temp netCDF files and perform diagnostic computations from them:
https://github.com/robertdemariacira/tc_diag_driver/tree/main

Instructions for how to use the code are included at the bottom of the README. Example configs and a manual test of the code are provided. Let me know if you have any questions.

[JohnHalleyGotway]

@robertdemariacira I've followed the instructions listed in the README file and am trying to figure out how to make a call to the driver code to have it compute diagnostics. I'd like to pass to it a NetCDF file for a single time step and see its output. Can you advise on that?

[JohnHalleyGotway]

Based on TC-Diag meeting on 9/20/23:

• Update TC-Diag to error out if handed a track to define diagnostics with an 'I' or a '2'.
• Update TC-Pairs to error out if the user requests to match diagnostics to track ending with an 'I' or a '2'.
  Note that using an 'I' or '2' will produce unexpected results.

[robertdemariacira]

@JohnHalleyGotway I've updated the tc_diag_driver repo to v0.6.0. It no longer requires the ATCF tech ID to be specified in the config file. You can find an example for how to call the driver here. The text distance to land LUT the code needs (for now) can be downloaded from here. You'll need to update the example config files described in the above readme section to point to the location of the LUT file.

[JohnHalleyGotway]

See the TC-Diag project-wide meeting notes from 10/4/23.

Relevant for the TC-Diag tool:

• @JohnHalleyGotway will define lists of standard storm data and sounding diagnostic names. Anything not in those lists goes in the CUSTOM DATA section at the bottom in the output file.
• @robertdemariacira will move the land_lut_file setting from the post_resample.yml to be an argument for the diag_calcs function.
• @robertdemariacira will update the python diagnostics code to compare the maximum radius requested in post_resample.yml to the actual cyl coord data handed to it. If the cyl coord radius < max config radius, it'll raise an exception.
• @JohnHalleyGotway will update compute_tc_diag.py to catch that exception.
• @JohnHalleyGotway will add 2 new config options to the TC-Diag domain_info section:

diag_config = "MET_BASE/python/tc_diag/config/post_resample.yml"; // for parent
diag_config = "MET_BASE/python/tc_diag/config/post_resample_nest.yml"; // for nest

Defines the python config file to be applied to that domain. Will be different for the parent and nest.

override_diags = [ "RMW", "SST" ]; // default setting only uses the nest RMW and SST

Defines which diagnostics computed for this domain should override previously defined diagnostics. The default value of an empty list means that all diagnostics returned from python should be used.

[robertdemariacira]

@JohnHalleyGotway The v0.7.0 TC Diag Driver code addresses the land_lut_file and radii validation issues described above.

The land_lut_file config parameter is now optional and can be removed from the config as long as the land_lut_override arg is provided to diag_calcs. If the config parameter is not removed, it will be ignored if land_lut_override is provided. See the README for an example.

I added a config parameter: radii_to_validate which is just a list of radii to check against the min/max radii found in the input file. A ValueError will be raised if any of the entries in radii_to_validate are found to be outside a valid range. It should be noted that the code could potentially raise ValueError, IndexError, KeyError, or IOError exceptions if other problems are found when the code is preparing to perform diagnostic calculations. These errors can't be suppressed when suppress_exceptions is set to True since they're raised when the code is figuring out information needed to pre-populate the returned data with missing values.

[JohnHalleyGotway]

@jvigh, @robertdemariacira, and @musgrave-kate, we have an issue with units that needs to be addressed.
Listed below are diagnostic names from a sample diag file (sal092022_avno_doper_2022092600_diag.dat). Each name is followed by the units used in the doper file and that's followed by the units reported by the python diagnostics code. I only listed the 1000mb sounding data for brevity.

You'll notice some simple differences in naming conventions (e.g. DEG vs DEGREES) along with more substantive differences like 10C versus K and 10KT vs M/S.

LAT (DEG) (DEGREES_NORTH)
LON (DEG) (DEGREES_EAST)
MAXWIND (KT) (KTS)
RMW (KM) (NA)
MIN_SLP (MB) (MILLIBARS)
SHR_MAG (KT) (M/S)
SHR_HDG (DEG) (DEGREES)
STM_SPD (KT) (KT)
STM_HDG (DEG) (DEGREES)
SST (10C) (K)
OHC (KJ/CM2) (KJ/CM2)
TPW (MM) (KG/M^2)
LAND (KM) (KM)
850TANG (10M/S) (M/S)
850VORT (/S) (/S)
200DVRG (/S) (/S)
T_SURF (10C) (K)
R_SURF (%) (%)
P_SURF (MB) (PA)
U_SURF (10KT) (M/S)
V_SURF (10KT) (M/S)
T_1000 (10C) (K)
R_1000 (%) (%)
Z_1000 (DM) (GPM)
U_1000 (10KT) (M/S)
V_1000 (10KT) (M/S)
TGRD (10^7C/M) (10^7C/M)

@robertdemariacira where are these unit discrepancies currently handled when the python code writes the doper output file?

@jvigh is it desirable for the NetCDF output to use precisely the same units as the existing ascii output or are you going to want them to differ?

How should we solve this? Options include:

• Update the call to post_resample_driver.diag_calcs() to return the exact diagnostic names (e.g. change current VMAX to MAXWIND) and units that should be reported in the output.
• Leave post_resample_driver.diag_calcs() as-is but update MET's python wrapper code to rename things and do unit conversion.
• Do the same as above, but in C++ rather than python (this is potentially messy and slow to be fixed).
• Abandon the approach of C++ writing the diagnostics ascii output and figure out how MET can leverage the existing python code to write it and, presumably, handle the discrepancies in naming conventions and unit conversion.

[robertdemariacira]

The diagnostic calculations are (mostly) performed using whatever units come from the input file. We then have a function that takes the Xarray Dataset objects returned by the diagnostic calculations and some configuration information and writes the doper file. The function uses the configuration information to map between the diagnostic variable name and the output name, the name of the output units, and performs any unit conversion needed.

[JohnHalleyGotway]

@robertdemariacira and @musgrave-kate, FYI, I was able to fix the unit problem for height with a simple change in the tc_diag_driver config files. Please see changes in this commit. Robert, I'd recommend that you make the same change in https://github.com/robertdemariacira/tc_diag_driver to keep the config files in sync.

[JohnHalleyGotway]

@robertdemariacira I'm wondering about scipy. It's used in the Python code here:

.//diag_lib/cylindrical_grid.py:from scipy import interpolate

And interpolate is used in:

def interpolate_to_cyl_grid(...)
...
class ScipyLinearNDInterpolator(CylindricalGridInterpolator)

I don't think that MET is actually using that SciPy interpolation functionality. I'm wondering if I should add SciPy as a new Python requirement even though MET is not actually using it. Or perhaps I should just comment out or delete the portion of the code that uses it.

Any suggestions?

[robertdemariacira]

@JohnHalleyGotway I don't think the code you're using should ever make a call to scipy. I could wrap the scipy import in a try/except clause and make it an optional dependency so that you don't need to include it as a dependency for MET.

[JohnHalleyGotway]

@robertdemariacira however you'd like to handle it is fine with me. But yes, I'd prefer not to add a dependency on scipy that we won't actually use. For now, I just modified this version of cyclindrical_grid.py by commenting out the scipy related code.

If you do implement a more elegant solution in your tc_diag_driver repo, please let me know, and I'll go grab the updated version.

Thanks!

[robertdemariacira]

@JohnHalleyGotway Sounds good. I'll let you know when I have that updated.

[robertdemariacira]

@JohnHalleyGotway I updated tc_diag_driver to v0.9.0. This version includes diag_lib v0.5.0 which now has Scipy as an optional dependency. I tested the driver in an environment where scipy wasn't installed and it appeared to run without any problems.