Problem getting Point-Stat to produce matched pairs for ADPUPA point obs

[jwraby]

I am having problems getting Point-Stat to produce matched pairs for ADPUPA point observations. My observations are as follows:
ADPSFC MSATWR10 20220705_120000 32.58 -106.67 1349.39 11 869.2578 1351.39 NA 298.97753906
ADPSFC MSATWR10 20220705_120000 32.58 -106.67 1349.39 52 869.2578 1351.39 NA 43.18756866
ADPUPA MSATWR10 20220705_120000 32.58 -106.67 1349.39 11 868.3243 1359.39 NA 298.93899536
ADPUPA MSATWR10 20220705_120000 32.58 -106.67 1349.39 33 869.2578 1351.39 NA 1.54999995
ADPUPA MSATWR10 20220705_120000 32.58 -106.67 1349.39 34 869.2578 1351.39 NA 0.39000000
ADPSFC MSATWR10 20220705_120000 32.58 -106.67 1349.39 33 868.3243 1359.39 NA 1.90000000
ADPSFC MSATWR10 20220705_120000 32.58 -106.67 1349.39 34 868.3243 1359.39 NA 1.05999994

The UPP post-processed GRIB WRF forecast file contains the following variables to match with the above observations:
ADPUPA variables for vertical level verification are:
3574:134172502:d=22070503:UGRD:2 m above gnd:9hr fcst:NAve=0
3575:134235316:d=22070503:VGRD:2 m above gnd:9hr fcst:NAve=0
3576:134298130:d=22070503:TMP:10 m above gnd:9hr fcst:NAve=0

The ADPSFC variables for surface verification are:
3505:132378496:d=22070503:RH:2 m above gnd:9hr fcst:NAve=0
3501:132178200:d=22070503:TMP:2 m above gnd:9hr fcst:NAve=0
3507:132484524:d=22070503:UGRD:10 m above gnd:9hr fcst:NAve=0
3508:132527738:d=22070503:VGRD:10 m above gnd:9hr fcst:NAve=0

The Point-Stat MPR file only has the ADPSFC pairs, but no ADPUPA pairs. I'm trying to figure out why no ADPUPA pairs.

I will provide my Point-Stat config file (PointStatConfig_v11.1.0_WRF_LES_twr10_obs_test11.txt), the forecast file (WRFPRS_d04.09_00), the Point-Stat log file (WRF_twr10_test11) and the NetCDF point observation file (d04.2022-07-05_12:00:00.0000_09FEB2024_test9.nc), and the MPR file (point_stat_WRF_twr10_test11_090000L_20220705_120000V_mpr.txt) to Julie Prestopnik via the Army file server (SAFE).

This test was conducted to see if the UPP GRIB forecast file containing the above output can be used to verify ARL tower data using surface verification and vertical level verification. This relates to GitHub issue:#2487

Thanks.

[DanielAdriaansen]

Hi @jwraby if you're using Python embedding, can you also share your Python script? Thanks!

[jwraby]

Hi @DanielAdriaansen for this testing I haven't used Python embedding.

[DanielAdriaansen]

OK thanks @jwraby, I got a little crossed up. I'll check with @jprestop and have a look at your files and get back to you. Thanks for all the details!

[jprestop]

@DanielAdriaansen and @jwraby I have grabbed the data and put it on seneca at /d1/projects/METplus/discussions/2498.

[JohnHalleyGotway]

@jwraby I see you have some questions about how Point-Stat constructs matched pairs. You are seeing MPR output when verifying forecasts of 2-meter relative humidity, 2-meter temperature, and 10-meter winds against observations with message type ADPSFC. However, you're not seeing MPR output when verifying 2-meter winds or 10-meter temperature against observations with message type ADPUPA. Hopefully I can provide some information to help clarify what's going on.

First, the Point-Stat tool provides options for handling the "level" information when matching forecast data to point observations. And those options are what need to be clarified here. Apologies, if I restate anything already obvious to you.

The Point-Stat configuration file includes separate dictionaries to specify information about the gridded fcst data and the point obs data. The fcst.field entry is an array specifying a list of gridded forecast fields to process. And the obs.field entry is an array describing the point observations to be used in that comparison. While we often set obs = fcst; to use the same settings for the forecast and observation data, they can be specified separately. The number of fcst.field entries must be identical to the number of obs.field entries, and they are matched up. So the first fcst.field gridded data is verified against the first set of obs.field point observations, and so on.

Each field entry is defined primarily by the name and level settings. For example, you're using name="TMP"; level="Z10"; to extract 10-meter temperature from your forecast file. If you used obs = fcst;, that same setting would be used for the point observation data as well. That setting tells Point-Stat to only use point observations of temperature when a height value exactly equal to 10.

Here's the usage statement from ascii2nc which describes the 11-column point observation format:

	The "met_point" ASCII format consists of 11 columns:
		Message_Type Station_ID Valid_Time(YYYYMMDD_HHMMSS)
		Lat(Deg North) Lon(Deg East) Elevation(msl)
		Var_Name(or GRIB_Code) Level Height(msl or agl)
		QC_String Observation_Value

		where	"Level" is the pressure level (hPa) or accumulation interval (HH[MMSS]).
			"Height" is meters above sea level or above ground level for the observation (msl or agl).

In the sample data you sent, the 9th column (i.e. Height) has a value between 1351 and 1359. That's not a value of 10, so you get 0 matches. None of the observations match the level value you requested. If you were to set the obs dictionary as name="TMP"; level="Z1300-1400"; then you would get matches since the height of the observation data would fall inside the range you requested.

But wait, if Point-Stat only uses point observations that match the requested levels, why are you getting MPR output for the ADPSFC message type? The answer is the message_type_group_map configuration option. By default, that setting includes this entry:

{ key = "SURFACE"; val = "ADPSFC,SFCSHP,MSONET"; }

And those three message types are assumed to indicate that these are point observations at the SURFACE. To make surface verification easier, Point-Stat includes a rule that when verifying a vertical-level-type forecast (i.e. name="TMP"; level="Z2";... with the Z in the level type) against a surface message type (i.e. message_type = "ADPSFC";), then it should ignore the vertical level matching requirement. And that's what's resulting in MPR's for the ADPSFC message type.

If you were to add ADPUPA to the list of surface message types, then you'd also see matched pairs for ADPUPA in the output. But I would not recommend that. ADPUPA is the PREPBUFR mnemonic indicating upper-air observations. Generally, they contain sounding data from weather balloons... and that's definitely not surface data.

It sounds like you'd like to know if/how you can verify 2-meter winds and 10-meter temperature against ADPUPA point observations. I'm guessing the real problem here is how the Height of those observations is defined. The forecast 2-meter wind and 10-meter temperature are defined as meters above-ground-level (i.e. AGL). I'm guessing the point observations are defined as meters above-mean-sea-level (i.e. MSL). And MET does not currently include any direct conversion between MSL and AGL for point observations. I certainly could see it could be enhanced to do so, using topography data provided by the topo_mask dictionary, but it definitely does not right now.

Big picture, I haven't really encountered folks verifying 2-meter winds and 10-meter temperature against ADPUPA point observations in the past. Instead, they verify forecast at pressure levels against ADPUPA point observations. For example, people often verify an 850 mb temperature forecast (name="TMP"; level="P850";) against a range of pressure level observations (`name="TMP"; level="P825-875").

Hope that helps clarify. What other questions do you have?

[JohnHalleyGotway]

@jwraby, thanks for confirming that that works as expected. Let me caution you about the pitfalls of this approach. While this gives you the 7 matched pairs you expect, if you run it with real data, it'll probably give you WAY TOO MANY matched pairs. ADPUPA observations usually contain an entire column of data with pressure levels for around 850 mb near the surface plus obs for at least all the mandatory pressure levels as the radiosonde rises through the atmosphere. Telling Point-Stat to count all of those as surface obs will make all of them appear in the list of matched pairs! And those will not be meaningful matches.

In practice, you'll need to some reasonable method for dealing with the vertical levels.

Options include:

1. If the topography is relatively constant in your domain, you could hard-code the observation MSL heights, like this:

fcst = { field = [ { name = "TMP"; level="Z10"; } ]; }
obs  = { field = [ { name = "TMP"; level="Z1300-1400"; } ]; }

2. Pick a pressure-level proxy for those heights instead (i.e. verify against pressure levels 825-875):

fcst = { field = [ { name = "TMP"; level="Z10"; } ]; }
obs  = { field = [ { name = "TMP"; level="P825-875"; } ]; }

3. The best solution would be modifying the ASCII point obs to convert the observation height value in column 9 from MSL to AGL. I think all you'd need to do is reset the value in 9-th column as {column 9 MSL obs height} - {column 6 MSL station elevation} = {AGL obs height}. And one option for that would be using a Python embedding script.

fcst = { field = [ { name = "TMP"; level="Z10"; } ]; }
obs  = { field = [ { name = "TMP"; level="Z9-11"; } ]; }

Question for you @jwraby. Do you think Point-Stat itself should be updated to be able to perform that operation? We could potentially add a config option convert_obs_height_msl_to_agl = TRUE; and then have the tool do that conversion when the data is read?

[jwraby]

@JohnHalleyGotway. Thanks for this input on providing these very practical ways of dealing with observations such as those I'm working with from the ARL tower array. The range in station elevations (MSL) for all the towers is 1313m to 1754m. The MSL obs heights are fixed for all towers at 2m and 10m AGL for my current dataset. For your best solution, since all the obs heights are the same (Z2 or Z10) couldn't you specify a single level of either Z10 or Z2 (depending on the variable) in the obs field setting rather than a range like in your example (Z9-11)? Just making sure I'm understanding how this works. To understand better how this solution works, I will edit the input obs file to change the values in col 9 from MSL heights to AGL heights and then run ASCII2NC and the Point-Stat to see how that works.

I like your suggestion of possibly modifying Point-Stat to do the conversion when the data is read. It might be a good idea to allow for the possibility of having obs at other levels on the towers besides 2m and 10m, but for the example you give, having obs at other AGL levels on a given tower would be transparent since you derive the AGL level by subtracting the MSL levels.

All of this discussion directly relates to aspects of the issue #2487

[jwraby]

Hi @JohnHalleyGotway. Just want to mention a word of caution regarding the model output on the non-standard levels. (wind components at 2m AGL and TMP at 10m AGL). Although I was able to get the UPP to put those values in the post-processed output by modifying the code, I need to verify that the values are accurate. I would hold off on modifying Point-Stat per the suggestion you made until I can confirm that those values are valid. IF they are not valid, then my fallback position is to verify the wind components and TMP at the non-standard levels on pressure levels per your option 2. above. In the issue #2487, I provided test case data where I ran Point-Stat by verifying against a pressure level using the tower data as ground truth. In that case I had to compute a pressure at the obs level 10m AGL because there is no pressure sensor at the 10m level on the towers. See the examples of the input point tower observations and the documentation I provided for #2487 .

[jwraby]

Does anyone know of a person who is familiar with the UPP WRF post-processor code? I'm running an older version V3.1.1. I would like to ask some questions about the code. Thanks,

[jwraby]

Due to the concerns within our modeling team about the validity of the post-processed WRF output for winds at the 2m AGL level and TMP at teh 10m AGL level I have reached out to the to GitHub Link: https://github.com/orgs/ufs-community/discussions/categories/q-a mentioned on the DTC website for a contact who might be able to answer questions about the UPP code. Until I can connect with a UPP expert who could answer our questions, I will consider the matter at hand in this discussion as closed. The information provided in this discussion is very helpful as I now have the answer to my original question about why I didn't have the ADPUPA matched pairs. In addition, the info provide above will be helpful if I can resolve the issue with the validity of the UPP as modified to output the above variables. Thanks to @JohnHalleyGotway, @DanielAdriaansen , and @jprestop for all your assistance with answering my question.

[jwraby]

@DanielAdriaansen. Thanks for checking. I may have contributed to the confusion by creating a new discussion, but if the problem can be resolved, then it paves the way for a better way to handle the tower data for the use case.

[jwraby]

Thank you, Julie.