dtcenter / MET

Model Evaluation Tools
https://dtcenter.org/community-code/model-evaluation-tools-met
Apache License 2.0
74 stars 22 forks source link

Bugfix: Fix the TC-Diag and TC-RMW tools to correctly handle the range and azimuth settings in range/azimuth grids #2833

Closed JohnHalleyGotway closed 3 months ago

JohnHalleyGotway commented 4 months ago

Describe the Problem

During development for issue #2729, I uncovered problems in how the range/azimuth grids are defined based on the configuration file settings. These grids are used by the by the TC-Diag and TC-RMW tools.

  1. There is an off-by-one error when defining the azimuths.

This issue is in the TcrmwGrid::set_from_data(const TcrmwData & d) function on this line of code:

# INCORRECT
RLLD.delta_rot_lon = 360.0/(Azimuth_n - 1);

The longitude rotation increment is defined as the number azimuths minus 1, and that - 1 should be removed. While doing development for #2729 and comparing the diagnostic output produced by TC-Diag to the output from the Python TC Diagnostics code, removing the - 1 made the previously inconsistent results much, much closer.

  1. There is an off-by-one error when defining the range increment.

Internally, MET computes and stores the maximum range value. Since the ranges always start with a range of 0 km, this line of code is correct. We should divide by n-1 when computing the range delta:

# CORRECT
RLLD.delta_rot_lat = range_max_deg/(Range_n - 1);

One bug is in the TC-RMW source code when defining that maximum range value on this line of code. We are multiplying by n_range when it should be n_range-1.

Another bug is in the TC-Diag source code when defining that maximum range value of this line of code. We are multiplying by n_range when it should be n_range - 1.

  1. TC-RMW always defines ranges relative to the radius of maximum wind and never as a fixed number of kilometers.

That same TC-RMW line of code from above contains another bug. Note that we're always multiplying by the radius of maximum winds to define the maximum range. The logic to differentiate between RMW ranges and fixed KM ranges is missing.

  1. The RMW grid configuration options are over-specified.

In the TC-RMW config file, specifying n_range, max_range_km, and delta_range_km (or rmw_scale) over-specifies the ranges. Only 2 of the 3 are actually required. After consulting with @KathrynNewman, we determined that max_range_km should be removed, and rmw_scale should be set to a default value of bad data (i.e. NA).

To be consistent with the existing documentation, if rmw_scale is set to something other than NA, its value should take precedence over the delta_range_km setting.

Note that the TC-RMW application code needs to be updated to implement this logic correctly.

Expected Behavior

The longitudinal rotation increment should be be computed simply as 360 divided by the number of azimuths (360.0/Azimuth_n). This differs from the latitudinal increment for which Range_n - 1 really should be used. Since the xml/unit_tc_rmw.xml tests use 180 azimuths, it's easy to overlook this very small difference in results since the 360/180 and 360/179 are very similar numbers.

To Reproduce

I tested this using the simplest setup I could imagine, using only 2 azimuths. With only 2 azimuths, the bug is obvious: RLLD.delta_rot_lon = 360.0/(Azimuth_n - 1) = 360 / (2 - 1) = 360 ... versus ... RLLD.delta_rot_lon = 360.0/Azimuth_n = 360 / 2 = 180

With this bug in place, rotating a fully 360 degrees should result in identical values being extracted for both azimuths. And indeed, that's exactly the behavior I found.

Describe the steps to reproduce the behavior: 1. Modify internal/test_unit/config/TCRMWConfig_pressure_lev_out by changing:

-n_azimuth      = 180;
+n_azimuth      = 2;

2. Run perl/unit.pl xml/unit_tc_rmw.xml 3. Observe the duplicate output values:

ncdump -v TMP test_output/tc_rmw/tc_rmw_pressure_lev_out.nc
 TMP =
  302.540002441406, 302.540002441406,
  302.615292718035, 302.615292718035,
  302.690582994665, 302.690582994665,
  302.2719353873, 302.2719353873,
  301.782571564577, 301.782571564577,
...

Setting n_azimuth = 2 should extract values from opposite sides of the storm but clearly the same values are being extracted here.

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

Labels

Milestone and Projects

Define Related Issue(s)

Consider the impact to the other METplus components.

Bugfix Checklist

See the METplus Workflow for details.

JohnHalleyGotway commented 4 months ago

@KathrynNewman, I am confident with this simple fix in the handling of the azimuths:

// MET #2833 divide by n rather than n-1 for the azimuth increment
RLLD.delta_rot_lon = 360.0/Azimuth_n;

However, I have questions about the handling of the ranges in the range-azimuth grid. Here's how the grid is defined in...

  1. The TC-RMW config file:
    n_range        = 100;
    n_azimuth      = 180;
    max_range_km   = 1000.0;
    delta_range_km = 10.0;
    rmw_scale      = 0.2;
  2. The TC-Diag config file:
    n_range        = 150;
    n_azimuth      = 8;
    delta_range_km = 10.0;

So TC-RMW supposedly has 3 options to specify the ranges... spacing in km, spacing in fraction of rmw, or the maximum range. But checking the actual code, while delta_range_km is parsed from the TC-RMW config file, I do NOT think its actually used at all! So the range can be defined using max range / n_range or by the rmw scale... but not specifying the delta_range_km directly.

Whereas TC-Diag ONLY supports defining the range by specifying delta_range_km.

If you think about it, specifying n_range, max_range_km, and delta_range_km in TC-RMW is an over-specification. And setting the latter doesn't actual even have an impact.

For simplicity, I'd be in favor of REMOVING the max_range_km config option from TC-RMW and updating the tool to support delta_range_km or rmw_scale.

The existing logic gets messy in regards to dividing by n or n-1 when handling the range rings. There's likely an off-by-one problem because the first range ring is always 0 km. But I think answering this config file question first is important and then we can discuss the handling of the n vs n-1 for the range rings.

Also, I'd recommend that we limit the fix for MET version 11.1.1 to only fixing the azimuths and make the larger change in the handling of the range rings for MET version 12.0.0.

JohnHalleyGotway commented 4 months ago

I note that testing TC-RMW before/after this fix reveals that the range and azimuth values reported in the output remain unchanged. However the differences can be seen in lat and lon NetCDF variables.

Tested in:

seneca:/d1/projects/MET/MET_pull_requests/met-12.0.0/beta4/MET-bugfix_2833_develop_azimuth/PR_2834/run_tc_rmw_tests.sh

This script does the following:

  1. Run tc_rmw for a single track point with a radius of maximum winds of 28 nm:

    AL, 14, 2016092900, 03,  fv3, 141, 233N,  760W,  73,  974, XX,  34, NEQ, 0133, 0147, 0082, 0116, 1006,  292,  28, ...

    And 28 nm = 51.856 km.

  2. Run tc_rmw using delta_range_km = 51.856; rmw_scale = NA;.

  3. Run tc_rmw using delta_range-km = NA; rmw_scale = 1;.

  4. Run ncdump on the output of both and prints the diffs:

    1c1
    < netcdf tc_rmw_KM {
    ---
    > netcdf tc_rmw_RMW {
    35c35
    <       range:units = "km" ;
    ---
    >       range:units = "fraction of RMW" ;
    107c107
    <  range = 0, 51.856, 103.712 ;
    ---
    >  range = 0, 1, 2 ;

    The good news is that all of the numbers are identical... other than the range units and values, as expected.

  5. Run tc_rmw with delta_range_km = 111.1; rmw_scale = NA;. This delta range in km is approximately the same as 1 degree. This should produce lat/lon points spaced roughly 1 degree apart, and that's what I do see in the output:

    range = 0, 111.1, 222.2 ;
    
    azimuth = 0, 90, 180, 270 ;
    
    lat =
    23.3, 23.3, 23.3, 23.3,
    24.2980278112276, 23.2962566691624, 22.3019721887723, 23.2962566691624,
    25.2960556224553, 23.2850284440261, 21.3039443775447, 23.2850284440261 ;
    
    lon =
    -76, -76, -76, -76,
    -76, -77.0866274956359, -76, -74.9133725043641,
    -76, -78.1731327440292, -76, -73.8268672559708 ;