Testing update #2893 in Pointstat (finally!!)

[robdarvell]

I’ve finally had a chance to look at update #2893 which was requested by me to allow the insertion of the Group Map value into the mpr (and orank) files rather than simply the key within the OBTYPE column of the output files.

I can confirm that setting the newly introduced variable name into the apps works as expected, ie it is using the ‘val’ rather than the ‘key’ in the output.

However, it is not putting in the correct value into the OBTYPE column in some instances.

I’ve been trawling through things like the pair_base and pair_data_point code trying to work out what is going on. The case which I have been looking at gave me the following output in my point_stat job: Please note that I have added the ‘typ:’ element to the code seen in the output below in trying to work out what is happening.

DEBUG 4: Computed NEAREST of 3 observations = 0 for [lat:lon:level:elevation] = [69.059:18.540:-9999.00:-9999.00]
DEBUG 4: [sid: 01023 typ: 10102 vld: 20250119_220000 obs: 102320]
DEBUG 4: [sid: 01023 typ: 10102 vld: 20250119_230000 obs: 102360]
DEBUG 4: [sid: 01023 typ: 10102 vld: 20250120_000000 obs: 102440]

Then the nearest ob to 00Z on the 20th is of course the bottom one in the list, and that is the one which appears in the output file, and the value in the OBTYPE column is 10102.

So that all looks good. Unfortunately, the ob at 12Z is of a different obs type to the other two (10101 rather than 10102).
Further up in the point_stat output (when I have changed the verbosity so I didn’t have to run too high) you can see this: (from around line 630 in pair_data_point.cc)

DEBUG 4: For air_pressure_at_sea_level(,) versus PMSL/Z0, skipping observation since it is a duplicate:
DEBUG 4: 10101 01023 20250120_000000 69.05900 -18.54000 76.00000 PMSL -9999.00000 -9999.00000 00016448 102440.00000

The way the code tests this, even though the header value is different, 10101 rather than 10102, the key it is being tested on is showing that it is the same as something previously (lat/lon/lvl/elev)??.
I think then the first occurrence of the key is then set within the station_values_t framework and then when things are set up to write out to the files then the value in the svt is used. (pair_base.cc approx. line 750), with the subsequent ob values for the relevant ob to be used.

I am wondering whether or not from this part of the code whether the typ_sa.add value needs to be taken from the ob.?? style structure, but that wouldn’t work at present as the type isn’t found within that (assuming it was ob_val_t). But then what other changes may be required I haven’t worked through. And then what that may or may not break elsewhere …

Of course I may be looking in totally the wrong places within the code !!

[JohnHalleyGotway]

@robdarvell, just saw your discussion come in. Thanks for the testing feedback!

I'm working on the first development release for the next cycle (MET-12.1.0-beta1) today. Wondering if you could share or point me to the sample data with which you're testing? I'd like to replicate this discrepancy and try to track it down. If it's "small" you could tar/zip it up and attach it to this discussion. If it's "large" you can use our Google Form Data Submission option.

[robdarvell]

Surface_20250120T0000Z.nc.gz
Fcast_20250119_0000_024_global.nc.gz

Attached a couple of files. Data is valid for 00Z on the 20th Jan, so a 24hr forecast.
I'm using the following FCST/OB settings

FCST_VAR1_NAME = air_pressure_at_sea_level
FCST_VAR1_LEVELS = "(,)"
FCST_VAR1_OPTIONS = set_attr_name="PMSL"; set_attr_level="Z0"
OBS_VAR1_NAME = PMSL
OBS_VAR1_OPTIONS = set_attr_level="PMSL"
OBS_VAR1_LEVELS = Z0

And the following group mapping:
POINT_STAT_MESSAGE_TYPE_GROUP_MAP = { key = "SURFACE"; val = "10101,10102"; }

Hopefully all the info you need is there (as you will have jobs set up to run these sorts of things)

As I said above, I think it maybe down to the duplicate key test not including the obtype part (which I can see why not), or when checking which ob to use (eg: when picking the nearest) the obtype isn't being transferred from the particular ob it selects, it is just using what it found when initially setting up the key structure.

[JohnHalleyGotway]

@robdarvell thanks for sending along this data. I narrowed it down a bit to just test the for the "01023" station which contains 4 observations for PMSL:

> cat Surface_20250120T0000Z_SID_01023_PMSL.txt
10101 01023    20250119_210000  69.0590   18.5400   76 PMSL             NA NA 00016448 102320.00000000
10102 01023    20250119_220000  69.0590   18.5400   76 PMSL             NA NA 00016448 102320.00000000
10102 01023    20250119_230000  69.0590   18.5400   76 PMSL             NA NA 00016448 102360.00000000
10101 01023    20250120_000000  69.0590   18.5400   76 PMSL             NA NA 00016448 102440.00000000

I spent a while testing the behavior for this station, varying the settings for:

• the time window defined by obs_window
• the duplicate logic defined by duplicate_flag
• the observation summary logic defined by obs_summary

I confirmed that varying the time window really does include/exclude obs based on their valid time as I'd expect.

For the duplicate_flag logic, I definitely see that the obs_key defined on this line of pair_base.cc does NOT currently include the message type. Point-Stat was originally designed to compute stats separately for each message type rather than grouping them together. So this hadn't come up previously - where we're combining 2 message types into a group.

Currently, when I run with:

obs_window = { beg = -6*3600; end = 6*3600; }
duplicate_flag = UNIQUE;
obs_summary = NEAREST;

I get a single matched pair:

DEBUG 4: Computed NEAREST of 4 observations = 6.952780267e-310 for [lat:lon:level:elevation] = [69.059:18.540:-9999.00:-9999.00]
DEBUG 4:   [sid: 01023 vld: 20250119_210000 obs: 102320]
DEBUG 4:   [sid: 01023 vld: 20250119_220000 obs: 102320]
DEBUG 4:   [sid: 01023 vld: 20250119_230000 obs: 102360]
DEBUG 4:   [sid: 01023 vld: 20250120_000000 obs: 102440]
DEBUG 2: 
DEBUG 2: --------------------------------------------------------------------------------
DEBUG 2: 
DEBUG 2: Processing air_pressure_at_sea_level(*,*) versus PMSL/Z0, for observation type SURFACE, over region 01023, for interpolation method NEAREST(1), using 1 matched pairs.
DEBUG 3: Number of matched pairs   = 1

What do you WANT to happen in this case? Should the obs for message type 10101 and 10102 be considered unique? Meaning, would you want basically want 2 matched pairs in this case:

• the 10101 obs closest to the forecast valid time of 20250120_000000.
• the 10201 obs closest to the forecast valid time of 20250120_000000.

Just trying to figure out what the goal is in this case.

I tried adding the message type to the obs_key in pair_base.cc:

   string obs_key = str_format("%s:%.3f:%.3f:%.2f:%.2f",
                       typ,         //  message type
                       lat,         //  lat
                       lon,         //  lon
                       lvl,         //  level
                       elv).text(); //  elevation

And reran point_stat to get 2 matched pairs instead of 1:

DEBUG 4: Computed NEAREST of 2 observations = 6.952798539e-310 for [message type:lat:lon:level:elevation] = [10101:69.059:18.540:-9999.00:-9999.00]
DEBUG 4:   [sid: 01023 vld: 20250119_210000 obs: 102320]
DEBUG 4:   [sid: 01023 vld: 20250120_000000 obs: 102440]
DEBUG 4: Computed NEAREST of 2 observations = 6.952798539e-310 for [messge type:lat:lon:level:elevation] = [10102:69.059:18.540:-9999.00:-9999.00]
DEBUG 4:   [sid: 01023 vld: 20250119_220000 obs: 102320]
DEBUG 4:   [sid: 01023 vld: 20250119_230000 obs: 102360]
DEBUG 2: 
DEBUG 2: --------------------------------------------------------------------------------
DEBUG 2: 
DEBUG 2: Processing air_pressure_at_sea_level(*,*) versus PMSL/Z0, for observation type SURFACE, over region 01023, for interpolation method NEAREST(1), using 2 matched pairs.
DEBUG 3: Number of matched pairs   = 2

Is that the behavior you expect in this case?

[JohnHalleyGotway]

In addition, I'll note that we were doing a bad job reporting the computed summary value! The code was updating a local copy of an object rather than the one stored in the underlying map object.

Here's the BAD LOG OUTPUT:

DEBUG 4: Computed NEAREST of 2 observations = 6.952798539e-310 for [message type:lat:lon:level:elevation] = [10101:69.059:18.540:-9999.00:-9999.00]
DEBUG 4:   [sid: 01023 vld: 20250119_210000 obs: 102320]
DEBUG 4:   [sid: 01023 vld: 20250120_000000 obs: 102440]
DEBUG 4: Computed NEAREST of 2 observations = 6.952798539e-310 for [messge type:lat:lon:level:elevation] = [10102:69.059:18.540:-9999.00:-9999.00]
DEBUG 4:   [sid: 01023 vld: 20250119_220000 obs: 102320]
DEBUG 4:   [sid: 01023 vld: 20250119_230000 obs: 102360]

And here's the FIXED VERSION in my local branch:

DEBUG 4: Computed NEAREST of 2 observations = 102440 for [message type:lat:lon:level:elevation] = [10101:69.059:18.540:-9999.00:-9999.00]
DEBUG 4:   [sid: 01023 vld: 20250119_210000 obs: 102320]
DEBUG 4:   [sid: 01023 vld: 20250120_000000 obs: 102440]
DEBUG 4: Computed NEAREST of 2 observations = 102360 for [message type:lat:lon:level:elevation] = [10102:69.059:18.540:-9999.00:-9999.00]
DEBUG 4:   [sid: 01023 vld: 20250119_220000 obs: 102320]
DEBUG 4:   [sid: 01023 vld: 20250119_230000 obs: 102360]

Where the bogus 6.952798539e-310 value is replaced by real values of 102440 and 102360. This is just an artifact of logging and does not actually impact the output values.

[JohnHalleyGotway]

Here's a patched MET-12.0.0 version of pair_base.cc.gz if you'd like to test with that updated logic to see if it that provides the functionality you'd expect.

[robdarvell]

Thanks John for the investigative work. I have downloaded the pair_base code you attached above and will have a play when I have chance.

In respect of what I am looking for (from within one of your comments above): As you noted there are 4 obs for station 01023, so using that as the example:

I set a time window around which obs we wish to include (in what we operate with in this case it is 2.5hrs before the VT and 30min afterwards. So for the 4 obs we just pick up the 22Z, 23Z and 00Z as the 21Z fails this test.
From those 3 successful obs we want to only include the observation which is nearest to the validity time, so in this case we just want the 00Z observation, which is the ob value which is being included in the mpr file (102440).
Then, from that 00Z report we also want to include the obtype (message type value rather than message type key) information code (in this case 10101) in the mpr.

We don't want to run the ob types as 2 separate setups as then as you have shown you would pick up the 23Z ob for 10102 and the 00Z for 10101. We only want the 00Z ob. We only wish to use the station Id once in what we store.

As the code already has the ability to find the observation which is NEAREST to a time I don't feel as though we need to do this extra processing outside of the code base. It does all I need, just without that last little bit. The correct ob is picked up by the code as it is, it is just not containing the correct obtype part of the observation required which is what I would like reporting.

There is then a follow-on query, along the same lines (but not to be tackled in this discussion) which is along the same sort of lines. And again this maybe different to anything which you have covered previously.
It is with the case of ships (or buoys or aircraft (but aircraft maybe different)) which move within the time window you request. Again, there is just the wish to include each Id once within the window (primarily the NEAREST). At present, each time the ship changes lat/long information within the window it becomes a different ob. We would again only like the ob which is closest to the valid time for the Id (not worrying about the lat/long changes directly). So we have a tendency to look at the Id rather than the lat/long (which I note doesn't always fit) when it comes to our obs processing.
But that is for another time.

[JohnHalleyGotway]

@robdarvell, OK, great. I understand the situation better now.

So the MET-12.0.0 released version does already pick out the single ob you'd like for the station named "102440". As for how the OBTYPE column is reported in the output, that is exactly what the new obtype_as_group_val_flag configuration option was added to control.

With obtype_as_group_val_flag = FALSE; I see in the MPR output line that the OBTYPE column is set to SURFACE. Re-running with obtype_as_group_val_flag = TRUE; I see in the MPR output line that the OBTYPE column is set to 10101. And this is exactly what you'd like, as I understand it. Please confirm that this really does work the way you expect.

I'll note that my proposed changes to pair_base.cc would BREAK this behavior. Including the message type in the "case information" would give you 2 separate pairs for the "102440" station: one for message type "10101" and another for message type "10102", and I understand that is NOT the behavior you want. You just want to use the ob whose time is closest to forecast valid time, and that can be either for "10101" or "10102". So there's really no need for you to compile/test with that updated file since it'd make things worse. Although, part of the changes in that file really do fix a legitimate bug in the log messages about the value being selected... but I'll fix that separately for MET version 12.1.0.

The issue about observations from stations that MOVE through time (like ships, buoys, flights, or even drifting radiosondes is really interesting), and I'm surprised this hasn't come up in the past. I suspect we'll need more control over how this "uniqueness key" is constructed on line 452 of pair_base.cc:

   //  build a uniqueness test key
   string obs_key = str_format("%.3f:%.3f:%.2f:%.2f",
                       lat,         //  lat
                       lon,         //  lon
                       lvl,         //  level
                       elv).text(); //  elevation

For "moving" observations, we'll likely want to replace latitude, longitude, and elevation with the station ID name instead. So making that uniqueness key definition configurable is probably what's needed. But yes, we should do that using a separate GitHub issue.