Unexpectedly high RMSE from grid_stat for wind fields vs temperature fields

[maddenp-cu]

Verifying some GRIB model output with grid_stat, I am expecting RMSE close to zero for all fields at time zero, but am consistently seeing unexpectedly high figures for u and v wind fields, and am wondering if these need special treatment somehow. I can provide a reproducer, but thought I'd start with a simple description in case the cause is obvious.

First, the forecast fields are shape (530, 900) and the truth fields are shape (1059, 1799). I have verified that, as expected, the even-numbered lat/lon values of the truth grid correspond to the lat/lon values of the forecast grid:

>>> np.all(ds_fcst.latitude.values == ds_u.u.latitude.values[::2, ::2])
np.True_
>>> np.all(ds_fcst.longitude.values == ds_u.u.longitude.values[::2, ::2])
np.True_

So, I suppose that

regrid = {
  method = NEAREST;
  to_grid = FCST;
}

in my grid_stat config should produce little error via regridding, since the nearest neighbor is at the same lat/lon location.

I did a manual RMSE calculation in the Python REPL with xarray and cfgrib and found comparable low values for both the temperature field at 800mb and the u-component-of-wind field at the same height:

>>> np.sqrt(np.mean((ds_t.t.values[::2,::2] - ds_fcst.t.values[8]) ** 2))
np.float32(0.0027963952)
>>> np.sqrt(np.mean((ds_u.u.values[::2,::2] - ds_fcst.u.values[8]) ** 2))
np.float32(0.0030535057)

The t value here matches what grid_stat calculates:

$ python -c "import pandas as pd; print(pd.read_csv('grid_stat_hrrrcast_t_isobaricInhPa_0800_000000L_20260313_000000V_cnt.txt', sep=r'\s+', usecols=['RMSE']))"
       RMSE
0  0.002796

But the grid_stat value for u is much larger:

$ python -c "import pandas as pd; print(pd.read_csv('grid_stat_hrrrcast_u_isobaricInhPa_0800_000000L_20260313_000000V_cnt.txt', sep=r'\s+', usecols=['RMSE']))"
      RMSE
0  1.70536

I'm guessing there's something fundamental that I'm ignorant of here. Any ideas? Again, I'm happy to provide a reproducer if it would be helpful. Thanks in advance!

[j-opatz]

Thank you for your question,

And I think we should dive a bit deeper into the problem you're experiencing. In order to do that, we will need to be able to recreate the issue on our systems. Given that you have already benchmarked some data, please provide those same datasets, along with the configuration set up and any relevant commands to run the setup that is producing these unexpected RMSE values, and we can starting digging into where the problem might be.

Here is a link on how to submit data to us.

Additionally, please be aware that several of our troubleshooting engineers are on leave this week, which may result in a slower response to this issue. But rest assured we will look into this and get back to you with updates or additional questions if we have any.

[maddenp-cu]

Thanks @j-opatz. The file is too large for either of the suggested methods, but please try this link. I was able to download it in a Chrome Incognito window, so I think it should be freely downloadable.

After expanding the .tgz file, you can execute the run script. I think the rest should be self-explanatory for MET experts, but please let me know if you have questions.

[j-opatz]

Hello,

Thank you for your patience as I reviewed the data you sent over thoroughly.

I was able to replicate the 0.0030535057 direct RMSE calculation for the UGRD variable field. I'm not convinced MET is producing a bad RMSE value, but there's more digging to do as this is not as simple as "MET doesn't apply the RMSE equation correctly to UGRD data".

I'll walk you through what I've tested so far, as well as what we can do next.

First I made sure that all of the data you're testing is orientated correctly by Plot-Data-Plane; if MET incorrectly flips the data, sometimes a Python check (like the one you performed) will catch it, whereas just using a statistical tool like Grid-Stat will just output odd numbers. But all of the data appears in the correct orientation. I want to note here that when you review the forecast UGRD field at 800mb, compared to the observation UGRD field at 800mb, I do see quite a bit of difference, especially along and off the east coast of CONUS. A bigger difference than what I saw in the TMP fields for 800mb from both sources. That has me suspicious that an RMSE value of 1.7 is closer to what should be expected than the 0.003 value provided from Python, but that's just a visual comparison.

Here is a copy of the UGRD forecast field:

And here is a copy of the UGRD observation field:

After confirming that the orientation of the data was correct, I ran the two configuration files you provided and confirmed the RMSE values that you provided (0.002796 and 1.70536 for the TMP and UGRD variables, respectively) are produced by Grid-Stat.

I then began to focus on Python comparisons where I would systematically remove MET from the calculations and see if the RMSE values were the same or different. I started by creating netCDF files using the Regrid-Data-Plane tool, where I regridded the observation fields to the forecast field using NEAREST at a width of 1. This regridding step is already performed in Grid-Stat, but I wanted to read those netCDF files into a Python script and remove Grid-Stat's calculation of RMSE. I used the following RMSE calculation:
print('the direct RMSE is',np.sqrt(np.mean((obs_data - fcst_data) ** 2)))
and the result produced the same RMSE values for the TMP and UGRD fields as Grid-Stat. I was then going to convert both of the observation files and the forecast file into a netCDF format, with no MET tools involved or regridding, to see what RMSE values were calculated.

That's where things got interesting.

I used cfgrib to create an output netCDF file for the HRRR UGRD variable data, and created a netCDF file from Regrid-Data-Plane where the HRRR UGRD data is "regridded" to itself (i.e. it should just convert the file type). As a quick check I did a print of the data fields in each file, just to make sure that they are more or less the same. But I saw some pretty significant differences between the two dataset's values!

Here's a small snippet of MET's netCDF data for the UGRD field (remember that this should be exactly what's in the grib2 file, since I'm using NEAREST method with a width of 1):

[[-4.628847 -4.682236 -4.7233887 ... -2.0220306 -2.020733 -2.0032055]
[-4.6282306 -4.69782 -4.7416234 ... -1.8964787 -1.878995 -1.8776833]
[-4.7112494 -4.6972013 -4.7571816 ... -1.8221624 -1.8208594 -1.8195564]
...
[ 9.271801 9.306747 9.286443 ... 19.654175 20.43404 21.310799 ]
[ 9.440822 9.471828 9.428335 ... 19.794317 20.425964 21.142052 ]
[ 9.518241 9.553062 9.675732 ... 20.03468 20.501877 20.982584 ]]

and here's the UGRD field in the netCDF created by cfgrib:

[[-5.48 -5.55 -5.61 ... -2.55 -2.55 -2.55]
[-5.48 -5.55 -5.61 ... -2.42 -2.42 -2.42]
[-5.55 -5.55 -5.61 ... -2.36 -2.36 -2.36]
...
[ 9.89 9.95 9.95 ... 10.95 11.52 12.14]
[10.02 10.08 10.08 ... 11.02 11.45 11.95]
[10.08 10.14 10.27 ... 11.14 11.45 11.77]]

That's a significant difference, I'd say! And just to be sure, I did the exact same netCDF creation process, but for the TMP field. And as you saw, those two netCDFs had field values that were identical. I did not test the observations, but I'd expect to see something similar in behavior.

Given that the cfgrib, Xarray method is producing a netCDF file with values that differ from MET's, my next step will be to discuss these results with the lead programmers and see what they suggest for additional testings. I appreciate your patience as we try to get to the bottom of this; I'll let you know when I know more.

[j-opatz]

Hi Paul,

I wanted to provide a second update on this, as the lead MET programmer has had some time to examine your test data and draw some conclusions.

Note that in his tests he used the P500 values, rather than the P800 values you and I had been examining. However, the conclusion he arrived at is valid for all U and V components.

He ran wgrib2 to extract UGRD and VGRD at P500, and then did the same but added -set_int 3 47 0 to reset the resolution and component flag to 0 to prevent the wind rotation from occurring.

When those outputs were complete, he ran both through Grid-Stat and looked at the difference fields in the netCDF output. What he saw seems to provide an explanation for why the RMSE from MET is different than the RMSE directly from the grib2 files using cfgrib and Xarray: when MET reads UGRD and VGRD from GRIB1 or GRIB2 files, it automatically checks the metadata and, if needed, rotates the winds from grid-relative to earth-relative. This rotation is NOT applied if the metadata does not indicate a grid-relative or earth-relative setting.

So while you would see this log message at -v 3 when reading in the HRRR UGRD variable:

DEBUG 3: Rotating U and V wind components from grid-relative to earth-relative.

You would NOT see that message for your observation fields. This difference between forecast rotation and observation non-rotation is the most likely culprit for why I was seeing differences between the Regrid-Data-Plane field values for the HRRR data and the HRRR netCDF output from cfgrib. And as you can assume, these field value differences will definitely affect the statistical output like RMSE.

Unfortunately, there's not many options at the moment to get around this wind component rotation issue in MET. The best method would be to utilize Python Embedding to pass in your forecast HRRR files to Grid-Stat, avoiding any rotations. This should provide the more realistic comparison (and smaller RMSE). However, I do have good news; there's an existing MET issue that seeks to expand on support for wind direction, including an ability to toggle the usage of wind rotation between grid- and earth-relative. I've made a note of this Discussion in that issue, so when work resumes on it we can ensure this situation can be addressed.

Please let me know if there are any additional questions or issues on this topic.

[maddenp-cu]

Thanks @j-opatz for this careful analysis!

I'm going to look again at the data, but am I correctly understanding that the forecast (not obs) GRIB data in this case is providing metadata indicating either grid-relative or earth-relative orientation and that, based on that, MET is performing a rotation? And is the issue that the metadata is incorrect, and that a solution might be to correct it? Or maybe a solution would be to remove this piece of metadata? I am in contact with the producers of the forecast GRIB data, so could work with them on correcting the GRIB data at generation time if that's the root cause of the issue here.

If correcting the forecast GRIB data is not a solution, would extracting the grid from the forecast GRIB and writing it as netCDF and having MET work from the netCDF potentially also solve the problem?

Running grid_stat via Python would be an invasive change in my application, so I would like to find another solution.

[maddenp-cu]

I'm looking at the data with ecCodes' grib_dump in octet mode, and for both the HRRRCast forecast data and the HRRR obs data, I see, in section 3, byte 47

47        resolutionAndComponentFlags = 8 [00001000 (grib2/tables/2/3.3.table) ]

consulting GRIB docs, and if I'm reading the bits in the right direction (bit 1 on the left, bit 8 on the right), and seeing bit 5 as 1, it seems that both GRIB messages are specifying

Resolved u and v components of vector quantities relative to the defined grid in the direction of increasing x and y (or i and j) coordinates, respectively.

i.e. grid-relative. If so, it seems that -- if the metadata isn't lying -- both grids are oriented the same way, so is MET rotating only the forecast grid and not the obs grid? You mentioned the Rotating U and V wind components log message earlier and wrote that

You would NOT see that message for your observation fields.

Does that only mean that no message is logged for rotated observation fields, or that rotation is not performed for observation fields?

Also, in re:

This rotation is NOT applied if the metadata does not indicate a grid-relative or earth-relative setting.

The bit has to be 1 (grid-relative) or 0 (earth-relative) so one or the other is indicated -- unless this section is optional and simply omitted for some grids. I guess there's no difference between grid-relative and earth-relative for lat/lon grids. In any case, for this grid type, not indicating grid-relative vs earth-relative in order to prevent MET from doing rotation might not be an option.

[j-opatz]

Thank you for doing a bit of digging on your end,

I think there's a few things that we can clear up, and some that I will need to defer to @JohnHalleyGotway for interpretation of the logic MET is using.

I can confirm that when I ran the UGRD data through Grid-Stat, only the forecast has the rotational correction applied; the observation field does not. The log statements will only be printed when an action triggering them occurs, so I was seeing

DEBUG 3: Rotating U and V wind components from grid-relative to earth-relative.

for the forecast's UGRD (so it was applying a rotation), but that log message did NOT appear for the observation's UGRD, so it is not being rotated. So any resulting comparison is not logically sound.

The reason for this inconsistent rotation logic is due to the contents of your observation file: you only have the UGRD component, so MET does not apply the rotational logic. When both the U and V components are available in a GRIB1 or GRIB2 formatted file, it applies the rotational correction. John HG will need to provide an answer as to why both components are necessary inside the file.

All of this actually provides a few non-Python solutions for the time being: you could make sure that both the forecast and observation files have the same number of U and V components. So in the case where both the forecast and observation files have U and V components, MET should rotate both fields; if they only have the U components, MET should not apply any rotation logic. Again, run the Grid-Stat tool at -v 3 to see a log message detailing the rotation.

You could also use netCDF inputs rather than GRIB2 format; MET does not apply any rotational logic to wind components inside netCDF files due to their non-standard formatting.

[maddenp-cu]

Very interesting that the behavior depends on whether both u and v are present in the same file! Yes, I'd love to understand why that's the case. Maybe if I thought about vectors for a bit it'd be clear to me, but if @JohnHalleyGotway or anyone else can beat me to it, I'd be grateful.

make sure that both the forecast and observation files have the same number of U and V components

This seems like the best immediate fix for my use case. I'll start looking into that.

[maddenp-cu]

I think I see now why u and v have to be present together: The resulting wind direction and speed cannot change, but there are many u / v combinations that could produce it, and they can only be changed in relation to each other.

[j-opatz]

Thanks for looking into that, Paul. Keeping the U and V components rotated at the same time makes sense if multiple combinations of the components result in the same second-tier wind variables.

Did you find if keeping the number of U and V components consistent across your forecast and observation data resolved your issue in METplus? If so, we can go ahead and mark this Discussion as resolved.

You are welcome to track dtcenter/MET/issues/1518 for updates on a more refined solution to this problem.

[maddenp-cu]

@j-opatz Yes, I can confirm that making sure the availability of u and v matches for truth and forecast GRIB files resolves the issue. In my case, since I am retrieving individual truth GRIB messages from remote GRIB datasets and so e.g. have a truth GRIB file containing a single u message, I am creating a corresponding single-message u forecast GRIB file, and now MET does no rotation and all's well. This is unfortunate because extracting messages from the forecast GRIB file uses additional disk space, so I appreciate that you made the case in dtcenter/MET#1518 for placing the rotation under user control -- at least having an override available for the implicit behavior would be a great help. Much appreciated. I subscribed to that issue and hope to hear more on it. Thanks again for your kind help here, and please feel free to close the issue when you like.