Fix the wind-based diagnostics computations in TC-Diag

[JohnHalleyGotway]

Describe the Enhancement

TC-Diag development history:

1. MET-12.0.0-beta2: Issue Enhance TC-Diag to actually compute and write diagnostics #2550 and PR Feature #2550 compute_diagnostics #2728 enhanced the TC-Diag tool to actually compute tropical cyclone diagnostics and write output files.
2. MET-12.0.0-beta3: Issue Enhance TC-Diag to use tc_diag_driver version 0.11.0 #2769 and PR Feature #2769 tc_diag_driver_v0.11.0 #2812 updated the Python driver code, addressed unit issues, and fixed the diagnostic file formatting.
3. Note that issue Enhance TC-DLand to use updated NHC land definition data #2751 will also impact the TC-Diag tool but is documented as a separate issue.

The following tasks were not completed for #2769 and should be addressed during the beta4 development cycle:

• Investigate why the computed diagnostics DO NOT MATCH @robertdemariacira. As discussed on 11/1/23:
  • It is OK if diag values differ by 1 due to rounding.
  • There are clearly problems with the winds.
  • @JohnHalleyGotway and @robertdemariacira will coordinate on testing/debugging.
• Update unit_tc_diag.xml to run a realistic HWRF parent + nest example. Be sure to test with with and without MET_PYTHON_EXE set. Use Hurricane Ian example for which we have diagnostics output and @KathrynNewman is providing gridded HWRF test data. Note that d01 should be used as the "parent" and d03 as the "nest" data source.
• MET-12.0.0-beta3 replaced writing temporary Python files in NetCDF format with JSON + NumPy serialization. TC-Diag still writes NetCDF cylindrical coordinate files for the Python tc_diag_driver code to read. Investigate whether or not this approach works on NOAA machines (like jet) where there was a NetCDF conflict. If not, consider always making a system call to Python (with MET_PYTHON_EXE set) rather than actually embedding it.

Additional task list (last updated based on Nov 29, 2023 project meeting):

• Keep NetCDF variable long names consistent with the ones listed in this Google Doc. See @jvigh.
• Enhance the NetCDF diagnostics output to report data in more typical units (e.g. use C rather than 10C). However, recommend keeping knots rather than switching to m/s.
• Run the TC-Diag NetCDF output through the https://ioos.noaa.gov CF-convention compliance checker. It may not be possible to achieve full compliance, but recommend getting as close as reasonably possible.
• Parallelize the processing of valid times. Simple OpenMP does not work. Global python instance is a problem.
• Update documentation @JohnHalleyGotway:
  • Describe CIRA and NetCDF diag outputs.
  • Update Appendix F with TC-Diag python embedding info.
• Review existing configuration options and documentation for details about vortex removal. Make any necessary updates since it will not be included in MET version 12.0. See issue Enhance TC-Diag to support vortex removal #2752.

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.
See related issues:

• Enhance TC-Diag to use tc_diag_driver version 0.11.0 #2769
• Enhance TC-DLand to use updated NHC land definition data #2751
• Enhance TC-Diag to support vortex removal #2752

Relevant Deadlines

List relevant project deadlines here or state NONE.

Funding Source

2770043

Define the Metadata

Assignee

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

Labels

• Review default alert labels
• Select component(s)
• Select priority
• Select requestor(s)

Milestone and Projects

• Select Milestone as the next official version or Backlog of Development Ideas
• For the next official version, select the MET-X.Y.Z Development project

Define Related Issue(s)

Consider the impact to the other METplus components.

• METplus, MET, METdataio, METviewer, METexpress, METcalcpy, METplotpy
  No impacts currently.

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 issue
  Select: Milestone as the next official version
  Select: MET-X.Y.Z Development project for development toward the next official release
• Iterate until the reviewer(s) accept and merge your changes.
• Delete your fork or branch.
• Close this issue.

[JohnHalleyGotway]

Here's a CF-compliance checker report provided by @jvigh:
IOOS Compliance Checker.pdf

[JohnHalleyGotway]

@robertdemariacira and @musgrave-kate please find the sample data we discussed at the meeting on 11/29/23 in the following tar file on the RAL public ftp server:
ftp://ftp.rap.ucar.edu/incoming/irap/met_help/for_robert/tc_diag_test_case.tar.gz

This includes the inputs, outputs, and temp files I get when running TC-Diag. Please let me know if you have any questions.

Here's the contents of the run_tc_diag.sh script:

#!/bin/sh

TC_DIAG=/d1/projects/MET/MET_regression/develop/NB20231130/MET-develop/bin/tc_diag

echo "input/gfs/gfs.0p25.2022092400.f000.grib2 \
      input/gfs/gfs.0p25.2022092400.f006.grib2 \
      input/gfs/gfs.0p25.2022092400.f012.grib2 \
      input/gfs/gfs.0p25.2022092400.f018.grib2 \
      input/gfs/gfs.0p25.2022092400.f024.grib2" > input/gfs_2022092400_file_list;

${TC_DIAG} \
-deck input/adeck/aal092022_OFCL_SHIP_AVNO.dat \
-data parent GFSO input/gfs_2022092400_file_list \
-config config/TCDiagConfig_ian \
-outdir output \
-log run_tc_diag.log \
-v 4

And here's a full directory listing of the tar file contents:

ls -lah tc_diag_test_case/*
-rw-r--r-- 1 met_test rap 1.9M Nov 30 08:58 tc_diag_test_case/run_tc_diag.log
-rwxr-xr-x 1 met_test rap  568 Nov 30 08:54 tc_diag_test_case/run_tc_diag.sh

tc_diag_test_case/config:
total 12K
drwxr-sr-x 2 met_test rap 4.0K Nov 30 08:55 .
drwxr-sr-x 6 met_test rap 4.0K Nov 30 08:57 ..
-rw-r--r-- 1 met_test rap 3.5K Nov 30 08:55 TCDiagConfig_ian

tc_diag_test_case/input:
total 20K
drwxr-sr-x 4 met_test rap 4.0K Nov 30 08:57 .
drwxr-sr-x 6 met_test rap 4.0K Nov 30 08:57 ..
drwxr-sr-x 2 met_test rap 4.0K Nov 30 08:57 adeck
drwxr-sr-x 2 met_test rap 4.0K Nov 30 08:55 gfs
-rw-r--r-- 1 met_test rap  229 Nov 30 08:57 gfs_2022092400_file_list

tc_diag_test_case/output:
total 13M
drwxr-sr-x 2 met_test rap 4.0K Nov 30 08:58 .
drwxr-sr-x 6 met_test rap 4.0K Nov 30 08:57 ..
-rw-r--r-- 1 met_test rap  13M Nov 30 08:58 sal092022_gfso_doper_2022092400_cyl_grid_parent.nc
-rw-r--r-- 1 met_test rap  13K Nov 30 08:58 sal092022_gfso_doper_2022092400_diag.dat
-rw-r--r-- 1 met_test rap  48K Nov 30 08:58 sal092022_gfso_doper_2022092400_diag.nc

tc_diag_test_case/tmp:
total 14M
drwxr-sr-x 2 met_test rap 4.0K Nov 30 08:57 .
drwxr-sr-x 6 met_test rap 4.0K Nov 30 08:57 ..
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f0_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f102_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f108_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f114_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f120_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f126_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f12_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f18_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f24_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f30_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f36_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f42_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f48_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f54_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f60_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f66_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f6_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f72_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f78_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f84_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f90_parent_32767_0_.nc
-rw-r--r-- 1 met_test rap 642K Nov 30 08:58 tmp_tc_diag_AL092022_GFSO_2022092400_f96_parent_32767_0_.nc

[robertdemariacira]

@JohnHalleyGotway @musgrave-kate Got it! Thanks! Also, I believe we gave you some incorrect information about the way that resampling to the cylindrical grid works yesterday (and very possibly before that). The Fortran version of the code (as well as the other Python driver) actually use bilinear interpolation from the 4 closest samples not nearest neighbor interpolation. This may explain some of the differences in output values we've been getting. In any case, I'll be diving into this data while I do some debugging.

[robertdemariacira]

@JohnHalleyGotway @musgrave-kate I have a new version (v0.11.0) of tc_diag_driver that addresses the following issues:

1. Missing SHR_MAG, SHR_HDG, TGRD values:
   These computations reference the computed U, V diagnostic variables so when I changed their names from u, v to U, V these computations were not longer able to find the variables they needed. This is an easy fix. In the config files the shear and TGRD entries now reference U, V.

2. Missing units:
   The code now initializes the returned Dataset variables with units defined in the config files so that even if a variable can not be computed, the returned Dataset will have all of the correct units. The div_vort_850 and div_vort_200 entries in the config files now have explicitly defined units of /S to ensure that these will appear in the returned results.

Check tests/post_resample_test/post_resample.yml and tests/post_resample_test/post_resample_nest.yml for relevant changes.

Remaining Issues:

1. Missing 850TANG, 850VORT, and 200DVRG values:
   These values are not missing for me when I run the code locally on the provided input .nc files, so I'm not entirely sure why they're missing in the provided output.

2. Differences in RMW and other values:
   I believe two different methods are being used to resample the input grib variables to cylindrical coordinates. The resampled variables in the provided temp .nc files differ from their counterparts in the Fortran/Python versions. In many variables, such as temperature, these differences are fairly minor. Temperature was usually only different by less than a degree, but in some cases was 5 degrees different. However, the UGRD, VGRD, UGRD_Z10, and VGRD_Z10 variables differ more significantly. I suspect that the wind vectors are simply more noisy so differences in resampling methods may make a larger difference.

[JohnHalleyGotway]

On 12/13/23 @robertdemariacira noted that we should switch from nearest neighbor interpolation to bilinear in the regrid dictionary of the TCDiagConfig_default file for consistency with the existing Fortran implementation.

[JohnHalleyGotway]

During the TC-Diag project meeting on Feb 7, 2024, we made the following decisions:

1. The first two columns of the diagnostics (name, units), including the space between those entries, should have a total width of 16 characters.

• Recommend that we remove the parenthesis surrounding the units string to have more space.
• Revise the existing logic for checking the diagnostic name and units length, and instead enforce their combined length.

2. Recommend running Python tc_diag_driver code with GRIB2 inputs and use this option to dump the cylindrical coordinate data.

During this meeting, we worked up a nice example that can be found in seneca:/d1/projects/TCDiag/al092022_test. The run_met.sh, run_python.sh, and run_python_dump.sh scripts run the MET TC-Diag tool and tc_diag_driver Python code directly to compute diagnostics using the exact same inputs. Dumping the cylindrical coordinate data, it looks like there's a pretty clear shift going on. Here's a screenshot to demonstrate. Across the top is the tc_diag_driver version of the cylindrical coordinates data for 850mb U wind. Down the left side is the MET output. Note that the values for (0,0) differ and there's a shift in the pattern of the red and blue blobs.

[JohnHalleyGotway]

Tested TC-Diag and confirmed that azimuths are currently computed in degrees clockwise from north, as seen by lat/lon's computed for 4 azimuths. Relative to the center point, these points are north - east - south - west:

DEBUG 1: Range 5 = 500, Azimuth 0 = 0, (lat, lon) = (18.6916, 71.9)
DEBUG 1: Range 5 = 500, Azimuth 1 = 90, (lat, lon) = (14.1555, 76.5325)
DEBUG 1: Range 5 = 500, Azimuth 2 = 180, (lat, lon) = (9.70843, 71.9)
DEBUG 1: Range 5 = 500, Azimuth 3 = 270, (lat, lon) = (14.1555, 67.2675)

[JohnHalleyGotway]

Moving this issue along to the MET-12.0.0-rc1 cycle and treating it as bugfix. Do careful testing and resolve any issues with unexpected output values or behavior. Do not add any new features, just fix any problems with the existing functionality.

[JohnHalleyGotway]

Revisiting these diffs by looking in seneca: /d1/projects/TCDiag/al092022_test and comparing met/output/sal092022_gfso_doper_2022092400_diag.dat to python/output/sal092022_avno_doper_2022092400_diag.dat, the unacceptable differences remain in the 850TANG, 850VORT, and 200DVRG diagnostics:

While many other small differences are present, these are certainly the big ones.

Log output from run_python.sh for computing 850VORT and 200DVRG for the first time:

2024-03-07 16:42:37,545;INFO;tc_diag_driver.diag_vars;550;Started computing divergence and vorticity u_name:u, v_name:v, level_hPa:850, radius_km:1000.000000
2024-03-07 16:42:37,737;INFO;tc_diag_driver.diag_vars;563;Finished computing divergence:0.000001 and vorticity:0.000002 u_name:u, v_name:v, level_hPa:850, radius_km:1000.000000
2024-03-07 16:42:37,738;INFO;tc_diag_driver.diag_vars;550;Started computing divergence and vorticity u_name:u, v_name:v, level_hPa:200, radius_km:1000.000000
2024-03-07 16:42:37,933;INFO;tc_diag_driver.diag_vars;563;Finished computing divergence:0.000005 and vorticity:-0.000020 u_name:u, v_name:v, level_hPa:200, radius_km:1000.000000

'rad_tan_area_avg': {
  'callable': 'tc_diag_driver.diag_vars.radial_and_tangential_area_average',
  'output_vars': ['850radial', '850tang'],
  'kwargs': {'u_name': 'u', 'v_name': 'v', 'level_hPa': 850, 'min_radius_km': 0, 'max_radius_km': 600}
},
'div_vort_850': {
  'callable': 'tc_diag_driver.diag_vars.divergence_vorticity',
  'output_vars': ['850dvrg', '850vort'],
  'kwargs': {'u_name': 'u', 'v_name': 'v', 'level_hPa': 850, 'radius_km': 1000}
},
'div_vort_200': {
  'callable': 'tc_diag_driver.diag_vars.divergence_vorticity',
  'output_vars': ['200dvrg', '200vort'],
  'kwargs': {'u_name': 'u', 'v_name': 'v', 'level_hPa': 200, 'radius_km': 1000}
}

This is similar to the configuration used by MET, but there is a small difference in adding unit conversion attributes:

'rad_tan_area_avg': {
  'callable': 'tc_diag_driver.met_diag_vars.radial_and_tangential_area_average',
  'output_vars': ['850RADIAL', '850TANG'],
  'kwargs': {'u_name': 'UGRD', 'v_name': 'VGRD', 'level_hPa': 850, 'min_radius_km': 0, 'max_radius_km': 600},
  'unit_converters': ['pass', 'tc_diag_driver.met_post_process.scale10'],
  'units': ['RADIANS', '10M/S']
}

Compared the range/azimuth grid data for 850 U-wind between MET and Python and the values are just different:
MET (vertical): /d1/projects/TCDiag/al092022_test/met/dump/tmp_tc_diag_AL092022_GFSO_2022092400_f0_parent_1481401_0_.nc
Python (horizontal): /d1/projects/TCDiag/al092022_test/python/dump/u_000h_0850_dump.nc

It looks like the data is off by one, the leftmost MET column appears in the bottom Python row instead of the top one as you'd expect. But since these are degrees, that wrap point should be arbitrary.

The progression of values don't match. In both, the lower-left corner is (0, 0).
Here's a screenshot of the UGRD/P850 field:

MET:

• From left to right, the rotation angles (in degrees) increase from 0, 45, 90, 135...
• UGRD/P850 values go from positive to negative.
• This is consistent with a clockwise rotation starting from due east, as described for the TC-RMW tool.

Python:

• From bottom to top, the rotation angles (in radians) increase from 0, .78, 2.35...
• UGRD/P850 values go from negative to positive.
• This is consistent with a counter-clockwise rotation starting from due east.

Perhaps this difference in rotation direction is the underlying source of the 850TANG differences?

[JohnHalleyGotway]

To test, I temporarily modified TC-Diag to switch from clockwise rotation to counter-clockwise. In TmpFileInfo::write_nc_data(...):

   DataPlane dp_rev = dp_out;
   int nx = dp_rev.nx();
   int ny = dp_rev.ny();
   for(int y=0; y<ny; y++) {
      for(int x=0; x<nx; x++) {
         int x_from = (nx - x) % nx; 
         dp_rev.set(dp_out(x_from, y), x, y);  
      }
   }
   dp_out = dp_rev;

And this did have the desired effect of reversing the data value written to the NetCDF output files:

However the resulting 850TAG winds are not fixed:

[JohnHalleyGotway]

Success.

Update the TC-Diag Python driver code to:

• Switch from clockwise to counter-clockwise
• Convert the MET azimuths from degrees to radians in scripts/python/tc_diag/tc_diag_driver/post_resample_driver.py:

    # MET #2729: Modify the azimuth values reported by MET
    #   - Subtract from 360 degrees to convert rotation angle
    #     from clockwise to counter-clockwise
    #   - Convert from degrees to radians
    azimuth_1d_deg = (360 - input_data[config.in_azimuth_name]) % 360
    azimuth_1d = np.deg2rad(azimuth_1d_deg)

To produce the desired result with an acceptable level of diffs: