Singularity HDF5 Point Obs Readin

[justintsu]

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)

[justintsu]

Here is the observation file split into two files. You can extract the full data file by running (post unzip)

cat NAF_raob_2023033000.tar.gz.a* | tar xzvf -

NAF_raob_2023033000.tar.gz.ab.zip
NAF_raob_2023033000.tar.gz.aa.zip

[DanielAdriaansen]

Greetings @justintsu. I see that you are having trouble with Python embedding for point observations using Point Stat. Thank you for providing your script, data files, Point Stat config file, and command.

Prior to MET v11.0.0, Python Embedding for point observations was not well supported in tools like Point Stat, which require a special format of point observations typically provided by the PB2NC tool. In MET v11.0.0 and v11.0.1 (the version you are using), in order to use Python embedding for point observations with Point Stat you'll need to call a special reformatter on your point_data variable in your Python embedding script so that it conforms to the format Point Stat expects. This was a temporary bridge solution to the permanent fix which is included in MET v11.1.0+, where this is handled automatically by MET. So just be aware the suggestion here only applies to MET v11.0.0 and 11.0.1.

If you cannot upgrade to the latest version of MET (v11.1.0), I would recommend the following changes to your Python script test_h5_embed.py:

1. On line 10, add: from met_point_obs import convert_point_data
2. Comment out line 226
3. Uncomment line 227
4. On line 228, add: met_point_data = convert_point_data(point_data)

Could you please try that and let us know how it works?

[DanielAdriaansen]

Can you try unsetting MET_BASE, by either logging out and back in to your shell or manually unsetting the value, and try this command again:

singularity exec /app/projects/NEPTUNE/software/singularity/met_v11.0.1.simg point_stat 'PYTHON_NUMPY' 'PYTHON_NUMPY=/p/home/fnmocalf/user_dirs/tsu/uv/scripts/COAMPS/test_h5_embed.py /p/work1/fnmocalf/user_dirs/tsu/data/obs/NAF_raob_2023033000.h5 tmp' PointStat2023033000NEST1a0121x0121 -v 5 -outdir ./out/point_stat/radiosonde

Are you setting MET_BASE to try and get this to work? Or are you setting MET_BASE because you're using a MET installation has had it's shared files moved somewhere else?

[justintsu]

looks like unsetting MET_BASE is working, my observations seemed to have made it into point_stat however all of them are being rejected due to obs var name:
DEBUG 2: Searching 2893618 observations from 624 messages.
DEBUG 2:
DEBUG 2: --------------------------------------------------------------------------------
DEBUG 2:
WARNING:
WARNING: check_hdr_str() -> output header column OBS_VAR is empty!
WARNING:
DEBUG 2: Processing airtmp/pre_000100 versus /P100, for observation type radiosonde, over region FULL, for inter polation method NEAREST(1), using 0 matched pairs.
DEBUG 2: Number of matched pairs = 0
DEBUG 2: Observations processed = 2893618
DEBUG 2: Rejected: station id = 0
DEBUG 2: Rejected: obs var name = 2893618
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 = 0
DEBUG 2: Rejected: quality marker = 0
DEBUG 2: Rejected: message type = 0
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

Yes, I was setting MET_BASE because it's shared filed were moved elsewhere (somewhere separate from the singularity package)

[DanielAdriaansen]

OK. If all of the shared files at MET_BASE were moved, then setting your MET_BASE should not produce this error. Are you sure you are using the correct location for MET_BASE?

Please let me know if you need guidance setting MET_BASE, but if not I'd recommend leaving it unset.

Can you modify the obs entry in this file PointStat2023033000NEST1a0121x0121 to look like this:

obs = {
    field = [
{name="";desc ="ARCTIC_Nest1a0121x0121";}
    ];
}

The level entry doesn't need to be in there, as it should be provided in your Python script in the lvl column of your DataFrame, which it looks like you are already doing:

 lvl = data[typ]['pres'][:]

Run the command again after and see if matches are made.

[justintsu]

I am not exactly sure why MET_BASE was set to what it was, I was provided this wrapper script with it pre defined and decided not to touch it. Removing MET_BASE seems to not break anything though so I will keep it out.

I've modified the config file and reran and got this error:

DEBUG 4: VarInfoFactory::new_var_info() -> created new VarInfo object of type "FileType_Python_Numpy".
WARNING:
WARNING: RegridInfo::validate() -> Resetting the regridding width from 4 to 2 for regridding method "BILIN".
WARNING:
WARNING:
WARNING: RegridInfo::validate() -> Resetting the regridding width from 4 to 2 for regridding method "BILIN".
WARNING:
ERROR :
ERROR : VarInfo::set_level_info_grib() - either level or GRIB_lvl_typ, GRIB_lvl_val1 and GRIB_lvl_val2 (if necessary) must be specified in field information
ERROR :

[DanielAdriaansen]

You are attempting to regrid your forecast data to an observation grid, but your observations are not gridded:

////////////////////////////////////////////////////////////////////////////////

//
// Verification grid
// May be set separately in each "field" entry
//
regrid = {
   to_grid    = OBS;
   method     = BILIN;
   width      = 4;
   vld_thresh = 0;
   shape      = SQUARE;
}

Try setting to_grid = NONE; instead. This may not be the only problem though. Could you also share this script:
/p/home/fnmocalf/user_dirs/tsu/uv/nrl/veriftools/model/coamps/read_COAMPS_general.py

[justintsu]

read_COAMPS_general.zip
Here is that file. I'll make the to_grid change and retry

[justintsu]

it seems like I should be able to modify my dataframe to have all the variable types in column 'var', and just run test_h5_embed.py once via point_stat and the obs block lines will access this dataframe according to what is in the config file:
obs->field->name

is my thinking correct?

[DanielAdriaansen]

When running point_stat with observation python embedding and multiple observations lines in the obs block in the config file, will point_stat attempt to execute the python embedding once or as many times as there are lines in the obs block?

I believe it does it one time for each time. So if you change to a new valid time, the Python embedding (HDF5 read) will occur again for that new time, but does not run each time a new line is encountered in the obs block for the current time.

Since the observation embedding is performed in the invocation of point_stat at the command line and since my observation processing python routine is only able to read in one variable into memory at a time, does this mean that I will need to run point_stat separately for each variable and hence point_stat will be reading in the same raob file once per variable?

If you are memory limited from a computing hardware perspective, then yes. But if not, I would suggest modifying your Python routine to read in all observations relevant for all lines in the obs block in the Configuration file for Point Stat, especially if you have to read the entire HDF5 file to memory.

it seems like I should be able to modify my dataframe to have all the variable types in column 'var', and just run test_h5_embed.py once via point_stat and the obs block lines will access this dataframe according to what is in the config file:
obs->field->name

That's right. Just give the "var" column in your DataFrame the name in each name= entry in the obs block and it should be able to find matches that way.

[justintsu]

thanks, I modified my script after my last post and deployed point_stat for all valid times/variables/pressure levels overnight. It looks like it is taking 6-7 hours to run each point_stat for a valid time, which is a little less than ideal. I've posted the cnt and config file for one of the valid times just because I am noticing some very large biases in the temperature, particularly at high and low elevations. I suspect perhaps a unit error so I am checking my code
tsu.zip

[DanielAdriaansen]

Keep in mind that while the observations are only read one time with Python for each time, your COAMPS Python script is called for each entry in the fcst block. So you may consider doing some profiling of how long that script takes to run. Is it global data? What is the nx/ny of the grid?

I'd suggest reducing your config file to a single time, single subset of variables at a single level. It's possible some variables are taking longer than others for some reason, but without some systematic profiling it's hard to know where the bottlenecks are, if there are any.

One question I'd ask is whether it is necessary for you to be setting MET_PYTHON_EXE? If you set this, MET will write out and read back in any data that is obtained by Python. It is only recommended to use this if there are Python packages not available (or not possible to be installed) into the Python version that was used to compile the MET software against.

If you are able to use the Python MET was compiled with to run your Python scripts, data is passed in memory from Python to MET and no write/read is performed.

For read_COAMPS.py I see:

import os
import sys
import re
import numpy as np
import datetime as dt
from IPython import embed as shell

I only see shell() used in one place, and it could probably be re-worked to not use that. All of the other imports here should be available in the Python3 standard library.

For your obs script, I see:

import h5py
from collections import OrderedDict
import pandas as pd

from IPython import embed as shell
import sys

import numpy as np
import math

#FOR PROGRESS BAR, COMMENT OUT IF YOU DONT HAVE MODULE
from tqdm import tqdm

So seems like you could get by without tqdm, and I don't see shell used anywhere here so that line can be removed. The only one that's special here is h5py. So, is it possible for you to check if the Python used when compiling MET has h5py installed?

If not, then you will need to set MET_PYTHON_EXE. Unfortunately there's no way to tell MET right now to use the compile time Python for FCST and the user's Python (at MET_PYTHON_EXE) for OBS so it will use MET_PYTHON_EXE for both which will force the writing/reading of the FCST data too.

Thus, it's probably advantageous given the scale of your task to see if you can switch to the compile time Python. Let me know if you have additional questions.

The mpr file should show the obs/fcst pair values to assist with units debugging, but I see that's off. You could turn it on and then run a single time for only 1 level for each variable if that would help.

[DanielAdriaansen]

Hi @justintsu, just following up on this. If you feel that we've resolved your issue, could you mark a response as the answer? Otherwise let us know if you need anything else. Thanks!