Using IODA2NC output with POINT_STAT

[davis-wesley-12345]

To begin, I pushed the data I'm using to your server under DAVIS_DATA.

I was able to process a netCDF file using IODA2NC pretty easily, but I'm just confused at how to configure point_stat in order to read in the correct levels.

raob_all_v1_20201215T1200Z.nc4 = input to IODA2NC

raob_all_v1_20201215T1200Z.nc = output from IODA2NC, input to POINT_STAT

PointStatConfig_dev-nep1.3-L89-baseline-Dec_2020121512 = Sample POINT_STAT config

Example of Obs portion I need help with

obs = {
field = [
{
name = "eastward_wind";
cat_thresh = [ NA ];
level = ["8500"];
},
{
name = "eastward_wind";
cat_thresh = [ NA ];
level = ["2500"];
},
{
name = "northward_wind";
cat_thresh = [ NA ];
level = ["8500"];
},
{
name = "northward_wind";
cat_thresh = [ NA ];
level = ["2500"];
},
{
name = "air_temperature";
cat_thresh = [ NA ];
level = ["8500"];
convert(x) = K_to_C(x);
},
{
name = "air_temperature";
cat_thresh = [ NA ];
level = ["2500"];
convert(x) = K_to_C(x);
}
];
}

Example of ncdump-h of IODA2NC output

netcdf raob_all_v1_20201215T1200Z {
dimensions:
mxstr = 16 ;
mxstr2 = 40 ;
mxstr3 = 80 ;
nhdr = 72257 ;
nobs = 706320 ;
nhdr_typ = 1 ;
nhdr_sid = 1 ;
nhdr_vld = 4132 ;
nobs_qty = 2 ;
obs_var_num = 8 ;
variables:
int hdr_typ(nhdr) ;
hdr_typ:long_name = "index of message type" ;
hdr_typ:_FillValue = -9999 ;
int hdr_sid(nhdr) ;
hdr_sid:long_name = "index of station identification" ;
hdr_sid:_FillValue = -9999 ;
int hdr_vld(nhdr) ;
hdr_vld:long_name = "index of valid time" ;
hdr_vld:_FillValue = -9999 ;
float hdr_lat(nhdr) ;
hdr_lat:long_name = "latitude" ;
hdr_lat:_FillValue = -9999.f ;
hdr_lat:units = "degrees_north" ;
float hdr_lon(nhdr) ;
hdr_lon:long_name = "longitude" ;
hdr_lon:_FillValue = -9999.f ;
hdr_lon:units = "degrees_east" ;
float hdr_elv(nhdr) ;
hdr_elv:long_name = "elevation" ;
hdr_elv:_FillValue = -9999.f ;
hdr_elv:units = "meters above sea level (msl)" ;
int obs_qty(nobs) ;
obs_qty:long_name = "index of quality flag" ;
int obs_hid(nobs) ;
obs_hid:long_name = "index of matching header data" ;
obs_hid:_FillValue = -9999 ;
int obs_vid(nobs) ;
obs_vid:long_name = "index of variable names at var_name" ;
obs_vid:_FillValue = -9999 ;
float obs_lvl(nobs) ;
obs_lvl:long_name = "pressure level (hPa) or accumulation interval (sec)" ;
obs_lvl:_FillValue = -9999.f ;
float obs_hgt(nobs) ;
obs_hgt:long_name = "height in meters above sea level (msl)" ;
obs_hgt:_FillValue = -9999.f ;
float obs_val(nobs) ;
obs_val:long_name = "observation value" ;
obs_val:_FillValue = -9999.f ;
char hdr_typ_table(nhdr_typ, mxstr2) ;
hdr_typ_table:long_name = "message type" ;
char hdr_sid_table(nhdr_sid, mxstr2) ;
hdr_sid_table:long_name = "station identification" ;
char hdr_vld_table(nhdr_vld, mxstr) ;
hdr_vld_table:long_name = "valid time" ;
hdr_vld_table:units = "YYYYMMDD_HHMMSS UTC" ;
char obs_qty_table(nobs_qty, mxstr) ;
obs_qty_table:long_name = "quality flag" ;
char obs_var(obs_var_num, mxstr2) ;
obs_var:long_name = "variable names" ;
char obs_unit(obs_var_num, mxstr2) ;
obs_unit:long_name = "variable units" ;
char obs_desc(obs_var_num, mxstr3) ;
obs_desc:long_name = "variable descriptions" ;

// global attributes:
:FileOrigins = "File / generated 20220718_185944 UTC on host narwhal03 by the MET ioda2nc tool" ;
:MET_version = "V10.1.2" ;
:MET_tool = "ioda2nc" ;
:MET_Obs_version = "1.02" ;
:use_var_id = "true" ;
}

[CPKalb]

Hi:

With MET's version netCDF for point data (like the output from IODA2NC), levels need to be values that either exist in the obs_lvl variable, or formatted using P for pressure levels, Z for heights, A for accumulations, etc.. Do you want the 850 and 250mb levels? If so, my suspicion is that you want "P85000" and "P25000" for the level values

I wasn't able to test this since I don't have a copy of the model data. But, does this work and produce suspect output for you? More information about the letter codes can be found here in this section.

Christina

[davis-wesley-12345]

That seemed to move me in the right direction, but now I'm getting a different error:

ERROR : Dictionary::lookup_string_array() -> array "message_type" does not contain StringType.

Looking in my ioda2nc output and searching for "message_type" I get the this:
char hdr_typ_table(nhdr_typ, mxstr2) ; hdr_typ_table:long_name = "message type" ;
Searching for "hdr_typ_table" I find this:
hdr_typ_table = "" ;
So it looks like point stat is expecting the hdr_typ_table to be populated with at least one message type string, but there's nothing there.

These are all RAOB obs that I'm using but it looks like they're only marked as such in the filename and since that part of the IODA2NC output is empty, it looks like it's breaking POINT_STAT.

I tried manually setting it in the config file with an empty string, but that didn't work. POINT_STAT errored out saying that message_type couldn't be empty.

I tried giving it a dummy name like RAOB, but that didn't either. POINT_STAT ran, but the output was empty and the debug showed the following:

DEBUG 2: Processing air_temperature/250.0 versus air_temperature/P25000, for observation type RAOB, over region EASTERN_PACIFIC, for interpolation method NEAREST(1), using 0 matched pairs.
DEBUG 2: Number of matched pairs   = 0
DEBUG 2: Observations processed    = 706320
DEBUG 2: Rejected: station id      = 0
DEBUG 2: Rejected: obs var name    = 618030
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  = 87668
DEBUG 2: Rejected: quality marker  = 0
DEBUG 2: Rejected: message type    = 622
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

[CPKalb]

Hi:

What is the message_type variable set to in your point_stat configuration file? Point-stat will filter out observations based on the values given in the message_type variable. So, if the values listed in message_type don't match any of the strings in your header_type_table, then it will filter out all of the observations. I see, from your output, that there is 622 points filtered based on message type. Does changing the following in your point stat configuration file fix the empty data?
message_type = [ "RAOB" ];

Christina

[davis-wesley-12345]

So, if the values listed in message_type don't match any of the strings in your header_type_table

And therein lies the rub. There's only one entry in my header_type_table and it's an empty string.

From my last entry:

Looking in my ioda2nc output and searching for "message_type" I get the this:

char hdr_typ_table(nhdr_typ, mxstr2) ; hdr_typ_table:long_name = "message type" ;

Searching for "hdr_typ_table" I find this:

hdr_typ_table = "" ;

It seems like the only option is to put in message_type=[ "" ]; since that's the only string point_stat is going to find in the hdr_type_table, but when I do that, I still get obs filtered out because of message type and end up with zero matched pairs.

[davis-wesley-12345]

Next, I tried to add it in to the ioda2nc output and then use that to overwrite the original file
(base) daviswj@narwhal09:/p/work1/daviswj/UV-123/uv/scripts/NEPTUNE/dev> ncap2 -s 'hdr_typ_table="RAOB"' raob_all_v1_20201215T1200Z.nc raob_all_v1_20201215T1200Z.nc4 ncap2: raob_all_v1_20201215T1200Z.nc4 exists---e'xit, o'verwrite (i.e., clobber existing file), or a'ppend (i.e., replace duplicate variables in, and add metadata and new variables to, existing file) (e/o/a)? a Segmentation fault (core dumped)

[CPKalb]

It looks like you are missing the specification of the output file in the ncap2 command. It should be added on to the end as follows:
ncap2 -s 'hdr_typ_table="RAOB"' raob_all_v1_20201215T1200Z.nc output_file_name_here

You can use the same filename as the input file and it will ask if you want to overwrite, or select a different name.

Christina

[davis-wesley-12345]

Then I tried amending the original file, the one that got processed by ioda2nc. That didn't work either:

(base) daviswj@narwhal09:/p/work1/daviswj/UV-123/uv/scripts/NEPTUNE/dev> ncap2 -s 'hdr_typ_table="RAOB"' raob_all_v1_20201215T1200Z.nc4 raob_all_v1_20201215T1200Z_postprod.nc Segmentation fault (core dumped)

[davis-wesley-12345]

`

It looks like you are missing the specification of the output file in the ncap2 command. It should be added on to the end as follows: ncap2 -s 'hdr_typ_table="RAOB"' raob_all_v1_20201215T1200Z.nc output_file_name_here
`

I don't think this will work because it's telling me that the dataset is a netCDF4 and that command will only work on netCDF3

When no output file is specified, ncap2 tries to write directly to the input file but this can only work with netCDF3 not netCDF4 input files.

[CPKalb]

It looks like you are pointing to the original ioda file in the ncap2 command. There is no hdr_typ_table variable in the original file so it's crashing. Does it work if you point to the file that is output from IODA2NC, raob_all_v1_20201215T1200Z.nc as follows?
ncap2 -s 'hdr_typ_table="RAOB"' raob_all_v1_20201215T1200Z.nc raob_all_v1_20201215T1200Z_postprod.nc

I'm checking in the one of our software engineers on whether it is possible to get this output directly when running ioda2nc, and I'll let you know when I have that answer.

Christina

[CPKalb]

Hi:
To get a hdr_typ_table other than an empty string, I had to run IODA2NC with the configuration file IODA2NCConfig_hdr_typ.txt. This has the following changes from the default configuration file:
Change the metadata map line { key = "message_type"; val = "msg_typ"; } to { key = "message_type"; val = "station_ob"; }

Then run:
ioda2nc raob_all_v1_20201215T1200Z.nc4 raob_all_v1_20201215T1200Z.nc -config name_of_config_file

This should produce multiple values in the hdr_typ_table. When I run it, my hdr_typ_table looks like this:
hdr_typ_table =
"n_rbraob",
"n_rlraob",
"n_pbpibal",
"n_plpibalB",
"n_plpibalD",
"n_plpibalA",
"n_plpibalC"

To retain all of these in point_stat, a key will need to be added to the message_type_group_map. I added the last key to the dictionary as shown below:
message_type_group_map = [
{ key = "SURFACE"; val = "ADPSFC,SFCSHP,MSONET"; },
{ key = "ANYAIR"; val = "AIRCAR,AIRCFT"; },
{ key = "ANYSFC"; val = "ADPSFC,SFCSHP,ADPUPA,PROFLR,MSONET"; },
{ key = "ONLYSF"; val = "ADPSFC,SFCSHP"; },
{ key = "LANDSF"; val = "ADPSFC,MSONET"; },
{ key = "WATERSF"; val = "SFCSHP"; },
{ key = "RAOB"; val = "n_rbraob,n_rlraob,n_pbpibal,n_plpibalB,n_plpibalD,n_plpibalA,n_plpibalC";}
];

Then, set
message_type = ["RAOB"];

"RAOB" could really be any name that you want to set it to, as it's just an identifier to retain all the hdr_typ_table values that are listed in the message_type_group_map key. However, it needs to the same name in the "key" as what is listed in message_type.

Christina

[davis-wesley-12345]

JUST saw this. Let me give it a whirl and get back with you

[davis-wesley-12345]

That got the job done, thank you!

[CPKalb]

Great. I'm going to mark this discussion as answered. If you have other questions or problems, feel free to start up a new discussion.

Christina