Singularity HDF5 Point Obs Readin

[fearon10]

Discussed in #2391

Hi,

I'm following this previous issue closely as I'm attempting to do something very similar, to use python embedding to ingest a point text file (a raob text file) into point_stat. I've had similar success of getting the data into point_stat but receive no matched pairs.

DEBUG 2: Processing airtmp/pre_0850.0 versus t/P85000, for observation type RAOB, over region FULL, for interpolation method NEAREST(1), using 0 matched pairs.
DEBUG 2: Number of matched pairs = 0
DEBUG 2: Observations processed = 297447
DEBUG 2: Rejected: station id = 0
DEBUG 2: Rejected: obs var name = 213861
DEBUG 2: Rejected: valid time = 0
DEBUG 2: Rejected: bad obs value = 0
DEBUG 2: Rejected: off the grid = 0
DEBUG 2: Rejected: topography = 0
DEBUG 2: Rejected: level mismatch = 82998
DEBUG 2: Rejected: quality marker = 0
DEBUG 2: Rejected: message type = 588
DEBUG 2: Rejected: masking region = 0
DEBUG 2: Rejected: bad fcst value = 0
DEBUG 2: Rejected: bad climo mean = 0
DEBUG 2: Rejected: bad climo stdev = 0
DEBUG 2: Rejected: mpr filter = 0
DEBUG 2: Rejected: duplicates = 0

As a preliminary step, I have used ascii2nc with python embedding to create a netcdf file -- this works well. I then attempt to ingest the netcdf file into point_stat and receive the result above of no matched pairs. My netcdf and config files (within point.zip) are attached along with my log file. My point_stat call is as follows:

singularity exec met_v12.0.0-beta1.simg point_stat PYTHON_NUMPY test.nc
PointStatConfig_IODA_navgem_v2.1.1_2021030400_000000 -outdir outdir -v 5 -log log.point_stat

point.zip

log.point_stat.txt

Thanks,
Matt

 

^{Originally posted by justintsu October 25, 2023}
Hello,

I am trying to get point observations (radiosondes for now) in an HDF5 file read into point_stat using embedding. I've written the script test_h5_embed.py which reads in the hdf5 file and rearranges the incoming raobs in a dataframe whose indices match what I read on the MET documentation

df_master = pd.DataFrame(columns=['typ','sid','vld','lat','lon','elv','var',
'lvl','hgt','qc','obs'])

Finally, the dataframe values are extracted and converted to a list and stored in a variable named met_point_obs (although I've also tried storing it in varialbes named point_obs_data and point_data, neither worked). The script is invoked with arguments specifying what particular observation to store in 'obs'. I've placed a temperature file along with this script attached to this discussion. The script is invoked like:

test_h5_embed.py NAF_raob_2023033000.h5 tmp

With the correct environmental variables set, I run:

singularity exec /app/projects/NEPTUNE/software/singularity/met_v11.0.1.simg point_stat 'PYTHON_NUMPY' 'PYTHON_NUMPY' PointStat2023033000NEST1a0121x0121 -v 5 -outdir ./out/point_stat/radiosonde

and get this non-descript error:

ERROR :
ERROR : process_obs_file() -> trouble parsing the python command PYTHON_NUMPY.
ERROR :

Code_Config.zip

(See below comment for the observation file)

[jprestop]

Hi Matt (@fearon10). I converted the METplus Issue to a METplus Discussion. We use Discussions to interact with our users about any questions they have and when work is identified for enhancements, bugfixes, new use cases, etc. we typically create Issues.

Thank you for this information. You mentioned that you are not receiving any matched pairs, which I also see. Based on the information in your log file and in your post, I see that 297447 observations were processed:

DEBUG 2: Observations processed = 297447

However, all 297447 observations are being rejected for one of the three reasons below:

DEBUG 2: Rejected: obs var name = 213861
DEBUG 2: Rejected: level mismatch = 82998
DEBUG 2: Rejected: message type = 588

I can see that you're looking for the obs name of t:

Processing airtmp/pre_0850.0 versus t/P85000, for observation type RAOB, over region FULL, for interpolation method NEAREST(1), using 0 matched pairs.

Running the following command, I can see that t is the first value for obs_var, with an index of 0:

...
data:

 obs_var =
  "t",  // obs_var(0,39)
    "mixing_ratio",  // obs_var(1,39)
    "u",  // obs_var(2,39)
    "v",  // obs_var(3,39)
    "ght";  // obs_var(4,39)

Looking at the obs_vid variable (index of variable names at var_name), I want to determine how many instances there are with a value equal to 0, which shows how many values there are for t. So, I created a file called out.txt which will output all of the values in the obs_vid array using the following command:

ncdump -fc -v obs_vid test.nc >& out.txt

To find a total of only the values with 0, and therefore the total number of t observations, I ran the following command:

grep "0," out.txt | wc -l

which gave me 83586. Going back to the log output, we see there were 297447 observations processed and 213861 were rejected due to the obs var name. 297447 observations processed minus 83586 t values gives us 213861 observations that do not have the name t just as we would expect. This explains why 213861 observations were rejected due to the obs var name not matching.

Regarding the 82998 observations rejected due to level mismatch:

DEBUG 2: Rejected: level mismatch = 82998

I see the following, along with other instances, in your log file:

DEBUG 4: For GRIB record number 26, the requested level values (850 and 850) do not fall within the current range of level values (500 and 500).
DEBUG 4: For GRIB record number 27, the requested level values (850 and 850) do not fall within the current range of level values (250 and 250).

It looks like the level you are requesting doesn't line up with current range of level values. Please revisit the level values that you are requesting vs. what is available.

Regarding the 588 rejected for the message type:

DEBUG 2: Rejected: message type = 588

I see in you are looking for a message type of RAOB in your configuration file, however, I see INNOVATION_RAOB in the test.nc file. Running:

 ncdump -v hdr_typ_table test.nc

produces the following:

...
data:

 hdr_typ_table =
  "INNOVATION_RAOB" ;

Please try modifying your configuration file and see if that helps.

[DanielAdriaansen]

Hi Matt,

Sorry for not being more clear. Inside your Python embedding script read_NRL_FLATFILE_CONVERT_CtoK.py, you have to set the level attribute of the dataset attributes. I am wondering what it is being set to inside the Python script. It would be the level setting in the attributes listed in this table:
https://met.readthedocs.io/en/latest/Users_Guide/appendixF.html#attributes-for-2d-gridded-dataplanes

[fearon10]

Hi Daniel,

I read a NAVGEM model output field in flat binary format with the read_NRL_FLATFILE_CONVERT_CtoK.py as follows:

./read_NRL_FLATFILE_CONVERT_CtoK.py /p/cwfs/projects/NEPTUNE/navgem_v2.1.1/2021030400/fields/airtmp_pre_0850.0_0000.0_glo1440x721_2021030400_00000000_fcstfld

In the reader script, I set the attributes including the level. I printed the settings shown below.
Data Shape: (721, 1440)
Data Type: dtype('float64')
Data Range: 234.5175498962402 to 302.00023498535154 K
Attributes: {'valid': '20210304_000000', 'init': '20210304_000000', 'lead': '0', 'accum': '00', 'name': 'airtmp', 'long_name': 'Air Temperature', 'level': 'pre_0850.0', 'units': 'K', 'grid': {'name': 'Global 0.25 Degree', 'type': 'LatLon', 'lat_ll': -90.0, 'lon_ll': 0.0, 'delta_lat': 0.25, 'delta_lon': 0.25, 'Nlat': 721, 'Nlon': 1440}}

The level is set to 'pre_0850.0'.

[DanielAdriaansen]

Thanks @fearon10!

The level value must be the same as the level of your observations. If you observations you want to match have a level of 85,000 Pa (850 hPa), then the level of your NAVGEM data must also be exactly 85,000. In this case, please change the level attribute from pre_0850.0 to 85000 (but still a string).

The name value must also be the same as your observations. Based on your point_stat config file, it looks like your observations have a name of t. Thus, please also change the name attribute from airtmp to t in your Python script.

I will also note that the documentation here:
https://met.readthedocs.io/en/latest/Users_Guide/appendixF.html#attributes-for-2d-gridded-dataplanes

Recommends that the accum and lead have a string type with format HHMMSS. You may also need to change these two attributes to be exactly 000000, but it may work the way you have it I'm not sure.

After making the changes to name and level, and possibly accum and lead if needed, run again and see if you get any matched pairs.

[fearon10]

Thanks Daniel.

I've used this format in the past and had success. However, I changed the level, variable, accum, and lead in the attribute dictionary just to check as you suggested. I get the same result of no matched pairs. Seems like there is something else happening?

Attributes: {'valid': '20210304_000000', 'init': '20210304_000000', 'lead': '000000', 'accum': '000000', 'name': 't', 'long_name': 'Air Temperature', 'level': '85000.0', 'units': 'K', 'grid': {'name': 'Global 0.25 Degree', 'type': 'LatLon', 'lat_ll': -90.0, 'lon_ll': 0.0, 'delta_lat': 0.25, 'delta_lon': 0.25, 'Nlat': 721, 'Nlon': 1440}}

[DanielAdriaansen]

Hi @fearon10, thanks for trying that suggestion. While there are no matched pairs, did the number of rejections for various reasons change? I would be curious to see if the number of name or level rejections differs between your tests, but perhaps the rejections are for another reason now? Could you share an updated log file?

If possible, could you also share your Python embedding script for the NAVGEM data and a sample NAVGEM file? Without those we're unable to completely reproduce your workflow. If you can't share those that's OK, but it may take several exchanges with suggestions for you to try to arrive at a solution, and for each suggestion and try on your end we'll need to updated log file from that run. Thanks!

[fearon10]

Hi Daniel,

Apologies for the slow response. Unfortunately, I can't share the NAVGEM reader script and data without going through public release. I did receive the same result, of no matched pairs, when changing the name and level for the NAVGEM data.

I attached my script for converting the raob file to netcdf, maybe there's something in there? I also attached my point-stat
log and config file -- all in point-stat.zip. I'm using MET V12.0.0, but decided to try MET V11 -- same result.

point-stat.zip

Thanks
Matt

[fearon10]

I also attached my python code for the observation portion (MET_embed_ascii.py) since it appears that the issue is coming from how this code is passing the data to point-stat. At the bottom of the .py file, I've tried both the following options:

#met_point_data = convert_point_data(point_data)
met_point_data = point_data

as this code is a modified version that posted in Justin Tsu's post (#2391) on the same subject. Both yield the same result/error.

pycode.zip

[DanielAdriaansen]

OK yes, my apologies. In MET-11.0.1, you'll need to use convert_point_data like Justin's code.

When you tried this line:

met_point_data = convert_point_data(point_data)

Did you also uncomment this line at the top:

#from met_point_obs import convert_point_data

The older version of PointStat required observations in the more complex NetCDF-like format that MET uses, however newer PointStat versions can take either the 11-column ASCII format or the NetCDF-like format and if it's the 11-column format it will convert it automatically for the user.

Either way, it does appear going to an older version eliminates the SEGFAULT because you're getting farther than before, so I suspect switching to MET-12.0.0-beta2+ would also help, and would avoid needing the special convert_point_data() function that v11.0.1 requires.

[fearon10]

HI Daniel,
I've tried a few more times with version V11.0.1. No dice. It will be a few weeks before our group here at NRL can switch to V12.0.0-beta2 or newer. In the meantime, I will use ascii2nc as a temporary step prior to point-stat. I really do appreciate all your assistance.
best wishes,
Matt

[jprestop]

@fearon10 Dan will likely follow up as well, but I wanted to mention that if you are going to use a beta release, I would suggested grabbing the most recent beta5 release. Also, here is a link to where you can view our release schedule for the upcoming beta6, rc1, and official releases, which I'll also paste in below for your reference, in case this impacts what your team would want to install:

• Beta6 releases are tentatively scheduled for 2024-10-02.
• Release Candidate 1 releases are tentatively scheduled for 2024-11-13.
• Official Release releases are tentatively scheduled for 2024-12-11.

[fearon10]

Thanks @jprestop for the information and suggestions on releases. This is really helpful.