How to encode a 3D level

[Marston]

Where and how to encode the vertical levels of a 3D array into a series of grib messages in the same file:

for lev in range(nlev):   
    msg = grib2io.Grib2Message(discipline=<int>,idsect=<list>)
    # what section holds the vertical level information?
    msg.end() # Add Section 8 (End Section).  This terminates the GRIB2 message in it binary encoded form.
    g2file.write(msg)

[EricEngle-NOAA]

Have you tried this workflow with v2.0 beta of grib2io? Some of the API has changed with v2.0.

[Marston]

I also have this v2 demo notebook available: https://github.com/NOAA-MDL/grib2io/blob/master/grib2io-v2-demo.ipynb
This is a very good start, thanks for this. But there is, at least in my non-grib mind, a major difference: we are taking Xarray objects and writing them to grib2. This is truly from scratch, so to speak :-)
I saw that in section 4 there is a level encoding. I am going to try and encode the different levels there.

[EricEngle-NOAA]

So glad you mentioned xarray 😃 . grib2io v2.0 will introduce an xarray backend. We are working on documentation specifically for the backend.

[Marston]

Have you tried this workflow with v2.0 beta of grib2io? Some of the API has changed with v2.0.

@EricEngle-NOAA

Tested the new grib2io 2.0.0b2.post1 with the following:

import grib2io as gb
import xarray as xr

data = xr.open_dataset(ifile, engine='netcdf4', group='galwem')
data = data.assign(lon = [l if l>0 else l+360 for l in data.lon])
data = data.sortby(data.lon)
# remove 9999 for investigation but put it back b4 writing to grib2
# data = xr.where(data == 9999, np.nan, data)
varb = 15
av = data.variables3d.isel(time=0, features3d=varb)
name = data.names3d.isel(features3d=varb).item()
levels = data.msllevels.values
dlat = av.coords['lat'].values
dlon = av.coords['lon'].values
lat_first = dlat[0]
lat_last = dlat[-1]
lon_first = dlon[0]
lon_last = dlon[-1]
time = data.coords['time']

# Grib2 encoding
discipline = 0 
id_section = [center, subcenter, 4, 0, 3, 
                      time.dt.year.item(),
                      time.dt.month.item(),
                      time.dt.day.item(),
                      time.dt.hour.item(),
                      time.dt.minute.item(),
                      time.dt.second.item(),
                      0,
                      0, 
                      1]
gb.GribMessage(discipline=discipline, idsect=id_section)

This gave the following error:

Exception has occurred: ValueError
Must provide GRIB2 Grid Definition Template Number or section 3 array
  File "/home/mward/contrail_write_grib/write_grib.py", line 68, in <module>
    msg = gb.Grib2Message(discipline=discipline, idsect=id_section)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: Must provide GRIB2 Grid Definition Template Number or section 3 array 

This is a change I think you alluded to earlier but I had not gotten to this point as yet.
I looked for examples of the API call in in the code but this is different python style for classes than I am used to. I cannot find the documentation in the code on how the class is instantiated and what are the mandatory arguments. I checked to see if there is a test for writing grib but there doesn't seem to be any in the code so far.
This old way of calling the Grib2Message() has obviously changed. What is the new way of calling the class?

[EricEngle-NOAA]

The biggest change from grib2io v1.x to v2 is the way you create a Grib2Message object from scratch. Please see the grib2io v2 demo IPython Notebook, specifically the section "Creating Grib2Message Object from Scratch".

[Marston]

Thanks. I saw that but it was not obvious at the time. Too focused on getting the NCEPLIBS to work.
This way seems similar to the way ECMWF does theirs. Makes sense now.

[Marston]

@EricEngle-NOAA

I'm studying the JN mentioned above and adjusted my code to follow it's example. But I'm having a bit of trouble understanding and updating sections 4 and 5 with my info. I've taken a NCEP grib and extracted a message record in order to use it as a template to write my grib. For the other sections, i.e., 1 and 3, updating the parameters is simple.
I could update sections 1 and 3 in this manner, i.e., through direct assignment:

nsg.ny = nlat
nsg.nx = nlon
nsg.latitudeFirstGridpoint = lat_first
nsg.longitudeFirstGridpoint = lon_first
nsg.latitudeLastGridpoint = lat_last
nsg.longitudeLastGridpoint = lon_last
nsg.gridlengthXDirection = xincr
nsg.gridlengthYDirection = -yincr

But problems arose when I tried the same in sections 4 and 5. I got around this by assigning my new parameters:

type_of_process = 0 # https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table4-3.shtml
# Section 4 encoding
pdt = {'fullName': name, 'units': units, 'shortName': shortname, 'leadTime': leadtime, 'unitOfFirstFixedSurface': 102,
       'valueOfFirstFixedSurface': 1.0, 'unitOfSecondFixedSurface': None, 
       'validDate': datetime.strptime(str(time.values[0]).split('.')[0], '%Y-%m-%dT%H:%M:%S'),
       'duration': None, 'level': levels[0], 'parameterCategory': param_category, 'parameterNumber': param_number, 
       'typeOfGeneratingProcess': type_of_process, 'generatingProcess': 133, 'forecastTime': leadtime,
       'typeOfFirstFixedSurface': levelType, 'scaleFactorOfFirstFixedSurface': 0, 'scaledValueOfFirstFixedSurface': 1,
       'typeOfSecondFixedSurface': 255, 'scaleFactorOfSecondFixedSurface': 0, 'scaledValueOfSecondFixedSurface': 0}

for attr,key in zip(msg.attrs_by_section(4), pdt.keys()):
    try:
        setattr(msg, attr, pdt[key])
    except:
        pass

But as you probably already know, some of the parameters cannot be set this way therefore a try:except was necessary. But I'm wondering if this is an efficient way to go about this.

Section 5 is particularly confusing:

dptt = {
    'numberOfPackedValues': nlat * nlon,
    'dataRepresentationTemplateNumber': 0, #- Grid Point Data - Simple Packing (see Template 5.0)
    'typeOfValues': 0, # - Floating Point
    'refValue': 0.0,
    'binScaleFactor': decimal_scale_factor,
    'decScaleFactor': bit_width,
    'nBitsPacking': 0}

#msg.section5 = np.zeros([nlat, nlon], dtype=np.float32)    
#msg.section5[0] = nlat * nlon
for attr,key in zip(msg.attrs_by_section(5), dptt.keys()):
    try:
        setattr(msg, attr, dptt[key])
    except:
        pass

Could I get some help in understanding the API in these two areas?
Whenever I try the above method, the packing fails.

[Marston]

One more thing, I'm not sure why some of the table parameters values are inside a list/array and some are not?
Take the https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_sect1.shtml, the first 2 values are not included in the list:

print(msg.attrs_by_section(1))
['originatingCenter', 'originatingSubCenter', 'masterTableInfo', 'localTableInfo', 'significanceOfReferenceTime', 'year', 'month', 'day', 'hour', 'minute', 'second', 'refDate', 'productionStatus', 'typeOfData']

What is the logic behind this setup? I think it is linked to the issues I'm having with the sections 4 and 5.

[EricEngle-NOAA]

Generally, the section length and number are not a part of the public API. Despite those being formal parts of each section, the user does not have to sorry about providing those values.

[EricEngle-NOAA]

OK. Lets take 1 section at a time. Lets start at section 4. What template do you need to use for section 4?

Keep in mind that not all metadata need to filled in by the user. For example, validDate does not need to be set because, depending on the section 4 template, it will be determine from other date/time metadata.

[Marston]

I think I'm getting the hang of this API. What is causing the confusion is the I'm using the tables on the website and trying to set/account for all of those, but as you pointed out, that is not necessary.
I've refactored section 1 to now be a bit more dynamic:

# Sect 1, develop and populate. Note the string used to create the reftime cannot encode microseconds
sect1 = [center, subcenter, 4, 1, 1, 
         time.dt.year.item(), time.dt.month.item(), time.dt.day.item(), time.dt.hour.item(), time.dt.minute.item(), time.dt.second.item(),
         datetime.strptime(str(time.values[0]).split('.')[0], '%Y-%m-%dT%H:%M:%S'), 0, 1]
names = ['originatingCenter', 'originatingSubCenter', 'masterTableInfo', 'localTableInfo', 'significanceOfReferenceTime', 'year', 'month', 'day', 
         'hour', 'minute', 'second', 'refDate', 'productionStatus', 'typeOfData']
index = 0
for name,value in zip(names, sect1):
    if name not in ['refDate']: 
        msg.section1[index] = value
    else:
        setattr(msg, name, value)
    index += 1

for k,v in msg.attrs_by_section(1,values=True).items():
    print(f'{k}: {v}')

For section 4, the I'm using is pdtn is 0: Analysis or forecast.

[EricEngle-NOAA]

I guess it would be beneficial to the masses if there was a mechanism to inform a user as to what attributes are required to be populated. Sections 0 and 1 are easy because they are static, but 3, 4, and 5 are dynamic according their template numbers. And even with section 5, some of the metadata, should not be edited by the user.

Let me think about an added feature re: this.

[Marston]

Exactly! This is what has taken me so long to try and figure out. The Python logic uses features that I'm not familiar with and I plan to take some extra classes and decorators and repr, for example. The nature of this API also causes a difference in its behavior if it is being run in a NB or in a python file.

[Marston]

@EricEngle-NOAA

Was thinking about this and I think you already have a similar method in place that can be used as an example.
for k,v in msg.attrs_by_section(3) returns a list, only this time create a method that collects the mutable attributes in the given section that can be updated using the setattr() method. Here's some pseudo code given there is only 1 mutable attribute in section 3, for example:

myvalues = {'key': value}
for attr, key in zip(msg.mutable_attrs_by_section(3), myvalues.keys()):
    setattr(msg, attr, myvalues[key])

[Marston]

@EricEngle-NOAA

Just want to follow up on #66 (reply in thread).
Should I promote this request to an issue so that it doesn't fall between the cracks?

[EricEngle-NOAA]

Please make this an issue. This will have to wait until after 2.0 release though.

[Marston]

OK. Will do. When will 2.0 be released?