Raw forecast values from rotated lat-lon grid in MET not matching VER values

[j-opatz]

Been playing with PointStat using the facility added with issue #1055 (Enhance the MET library code to read Rotated Lat/Lon data from CF-compliant NetCDF files. #1055)
Was playing with our 1.5km model data (which is on a rotated pole). Have run it through GridStat and the output seems to match up OK with my comparator data (VER). So moved onto PointStat (as we need this to do some work with the HiRA code).
METplus PointStat job set up and it runs through fine. Was then trying to compare raw data values extracted through the MET run and something from VER. Obs match up OK but the forecast values didn't seem to be anywhere close.
I have subsequently then run PlotPointObs using the netcdf file which I was using in the PointStat job. I have attached (hopefully) a screenshot of the output. The image on the left was created using a slightly different model which isn't rotated. I have circled the observation location which I was trying to look at (note the model domain is slightly larger than the rotated pole model). The middle plot is then the same job but run with the rotated pole netcdf file, and the right-hand image is the rotated pole netcdf file image. On the middle image I have again circled the ob in question, but it appears as though there is a shift to the left compared to where things should be. On the right-hand plot I have shown a circle where the ob should be coming from and an X where the forecast value is being picked up from. So you can see they are not the same. There is no re-gridding to take place, so I have not included any of that in my job, so my question is, is there a setting which you need to include in your script somewhere when using rotated lat/long model grid, or is something not quite working correctly and is shifting things off centre (as shown in the plotpointobs output

)

[j-opatz]

I'm posting a Discussions issue on behalf of @robdarvell. This could generate some METplus community discussion and if an issue in the coding is discovered, will provide better documentation to where it originated from.

I'm going to tag @JohnHalleyGotway on this issue in the hopes that he can provide some guidance on how we can troubleshoot this discrepancy in point locations. He's currently out of the office but will respond as soon as he's able. In the meantime, thank you for your patience.

[mpm-meto]

Further to this discussion @JohnHalleyGotway what happens when you're in gridded space? I ran grid_diag to get some joint distributions, and I wasn't entirely sure it's correct. A little stumped as to how to test this but I don't think that what I've done should work!?
The plot on the left is using model fields only, both on rotated grids. I know there is 1;1 line on this plot (and it shouldn't be there!) but you can sort of see the two variables are related. The plot on the right is a comparison between the same two variables but based on observations. The y-axis data are not rotated and the x-axis data are on rotated model grid. I specified that I wanted to use the coarser non-rotated grid for regridding. As you can see the observations do not appear to be correlated at all, so I'm not sure this is working properly. I'm not sure that MET will correctly unrotate and regrid as this isn't a trivial thing to do (each grid point needs to be converted separately so that the unrotated grid isn't regular in and of itself, though the coarser unrotated grid IS regular. [Added: I suspect that going the other way (i.e. specifying the rotated grid as the reference) may well be a better bet. I shall try that and report what I find.]

Testing which direction to do the regrid yields largely the same result. Now, this may be genuine (i.e. that there isn't a relationship in nature whilst the model thinks there is one but, we need to be sure!) Left is regridding to the higher resolution rotated grid and right is regridding to the unrotated grid.

[JohnHalleyGotway]

I've begun looking into this issue. I don't have an explanation yet, but thought I'd share my progress. I used a GRIB2 rotated lat/lon file for testing and ran plot_point_obs to plot fields of latitude and longitude values. Then I overlaid the location of some large cities in northern Europe, plotting the values as the latitude or longitude. The resulting images look fine.

There's no obvious problem with the location of the dots... all the cities are in the right spots. And the latitude or longitude "value" inside each dot blends nicely into the background color.

The plot_point_obs tool calls the latlon_to_xy() function for the rotated grid, and this confirms that there's no obvious problem in that functionality. I'll keep digging.

[robdarvell]

@JohnHalleyGotway - I have put the two files you have requested onto the ftp server. Hope these are of use.

[robdarvell]

Have run plot_data_plane and have got the following out:

[robdarvell]

Not sure whether or not this is of any help to the issue, but I have spent a couple of hours today repeating the process for other rotated pole models which we run. My observation for those is that what is being seen above is not being seen on the others. I will try and paste other images in later on (not currently possible at present due to where I am writing this).
The only difference which I have spotted between the ones which appear to be ok and the one above seem to be in the north_pole_longitude values. The one above has a value < 180 (ie 177.5) whereas the others which I have processed are all above 180.

[robdarvell]

[robdarvell]

Example from our Euro model: Rotation 193 long, 41 lat

[mpm-meto]

I've run regrid_data_plane separately and plot_data_plane and a screenshot of the outcomes is attached. An interesting issue is that neither ncview nor Iris will display or read the resulting netcdf file when going from unrotated to rotated. Iris is really unhappy. See here.

cube = iris.load_cube("unrotated_to_rotated.nc")
Traceback (most recent call last):
File "/net/project/ukmo/scitools/opt_scitools/environments/default/2021_03_18-1/lib/python3.6/site-packages/iris/init.py", line 337, in load_cube
cube = cubes.merge_cube()
File "/net/project/ukmo/scitools/opt_scitools/environments/default/2021_03_18-1/lib/python3.6/site-packages/iris/cube.py", line 402, in merge_cube
proto_cube.register(cube, error_on_mismatch=True)
File "/net/project/ukmo/scitools/opt_scitools/environments/default/2021_03_18-1/lib/python3.6/site-packages/iris/_merge.py", line 1358, in register
match = cube_signature.match(other, error_on_mismatch)
File "/net/project/ukmo/scitools/opt_scitools/environments/default/2021_03_18-1/lib/python3.6/site-packages/iris/_merge.py", line 465, in match
raise iris.exceptions.MergeError(msgs)
iris.exceptions.MergeError: failed to merge into a single cube.
cube.standard_name differs: None != 'longitude'
cube.long_name differs: 'Precipitation Amount' != 'longitude'
cube.var_name differs: 'precipitation_amount' != 'lon'
cube.units differs: Unit('mm') != Unit('degrees')
cube.attributes keys differ: 'accum_time', 'accum_time_sec', 'valid_time', 'init_time', 'level', 'init_time_ut', 'valid_time_ut', 'name'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "", line 1, in
File "/net/project/ukmo/scitools/opt_scitools/environments/default/2021_03_18-1/lib/python3.6/site-packages/iris/init.py", line 339, in load_cube
raise iris.exceptions.ConstraintMismatchError(str(e))
iris.exceptions.ConstraintMismatchError: failed to merge into a single cube.
cube.standard_name differs: None != 'longitude'
cube.long_name differs: 'Precipitation Amount' != 'longitude'
cube.var_name differs: 'precipitation_amount' != 'lon'
cube.units differs: Unit('mm') != Unit('degrees')
cube.attributes keys differ: 'accum_time', 'accum_time_sec', 'valid_time', 'init_time', 'level', 'init_time_ut', 'valid_time_ut', 'name'

Also, the rotated_to_unrotated domain seems truncated, which is a bit odd. I would expect the unrotated to give me back the full domain?

Though hard to see it does look as though values are more or less in the right place in their interpolated and rotated/unrotated state. That's why I wanted to read them into Iris to check.

The India netcdf file I uploaded before has both grids (rotated and unroated) stored in the header and will give the extent that should come out from unrotation. The coords were created using iris like this.

import iris.analysis.cartography as iac

rot_lat_coord = cube[0].coord('grid_latitude').points
rot_lon_coord = cube[0].coord('grid_longitude').points
rotated_lons, rotated_lats = np.meshgrid(rot_lon_coord, rot_lat_coord)

pole_lon = cube[0].coord('grid_longitude').coord_system.grid_north_pole_longitude
pole_lat = cube[0].coord('grid_latitude').coord_system.grid_north_pole_latitude

true_lons, true_lats = iac.unrotate_pole(rotated_lons, rotated_lats, pole_lon, pole_lat)

co_x = iris.coords.AuxCoord(true_lons, standard_name="longitude", units="degrees")
co_y = iris.coords.AuxCoord(true_lats, standard_name="latitude", units="degrees")

[JohnHalleyGotway]

@mpm-meto and @robdarvell, wondering about the status of this discussion. Are you waiting for some additional debugging/testing from the METplus team? Or are you still planning to test out recent changes to MET?

Note that as of PR dtcenter/MET#2161 for issue dtcenter/MET#1920, the structure of the MET repository has changed a lot. And it should now be much easier to compile MET directly from the repo. Here's how:

git clone https://github.com/dtcenter/MET
cd MET
git checkout develop
./configure --prefix=`pwd` # ADD YOUR CONFIG OPTIONS HERE
make install

[robdarvell]

We are still awaiting the person who is going to test this out through the HiRA code to carry out the work. I have done a very preliminary examination of the new output, and it does seem to match up things now as expected, but I haven't tested it in any great detail.

[JohnHalleyGotway]

Hi @robdarvell and @mpm-meto, I'm moving some summary information over from from Slack and email to this discussion. Please see below.

Marion had asked me in Slack what HiRA logic exists in Ensemble-Stat, as of the MET-11.0.0-beta3 release on 9/21/22. I said nothing but figured it'd be wise to take a closer look and clarify. Here are the HiRA issues related to ensembles along with a brief description:

MET#1583 is DONE and was included in the MET version 10.1.0 release. This enables interp.type.method to be set to "HIRA", as seen on this line of code. I've attached an "_ecnt.txt" file to this email to illustrate the effect. This test runs ensemble_stat with 6 ensemble members. You see that in "N_ENS = 6" in that output file. However, the lines with "INTERP_MTHD = HIRA_SQUARE" and "INTERP_PNTS = 4" have "N_ENS = 24" instead. The ensemble is defined as size 24 for HIRA since it's 6 ensemble members * 4 points for each member.

If you have not already tested out that functionality, please do so.

MET#1602 is NOT DONE but is assigned to the current, beta4, development cycle. Now that ensemble post-processing has been removed from Ensemble-Stat (in beta3), this issue needs more definition. Perhaps we should actually enhance Gen-Ens-Prod to compute HIRA-based probabilities and write them to the NetCDF output file? Or since Ensemble-Stat can compute relative frequencies on the fly and verify them now, do we want to enable HIRA-based ensemble probabilities and verification all in one step in Ensemble-Stat?

Please add comments to that issue to clarify what is needed.

MET#1583 is NOT DONE but is assigned to the current, beta4, development cycle. Can you confirm that we should add this new statistic to the ECNT line type for MET version 11.0.0?

[JohnHalleyGotway]

Since MET#1583 was included in the MET-11.0.0-beta4 development release, I'll mark this discussion as being answered and lock it.