Include Grib2 format capability

[fossell]

Describe the New Feature

UPP output is in grib2 format. Need to assess the feasibility of reading in grib2 files for plotting.

Acceptance Testing

Example UPP grib2 output data can be found on Cheyenne here. Tests should include verifying grib2 data can be read and used in downstream plotting tools.

Time Estimate

1-3 days

Sub-Issues

Consider breaking the new feature down into sub-issues.

• [ ] Add a checkbox for each sub-issue here.

Relevant Deadlines

Sep 30, 2022

Define the Metadata

Assignee

• [x] Select engineer(s) or no engineer required
• [ ] Select scientist(s) or no scientist required

Labels

• [x] Select component(s)
• [x] Select priority

Projects and Milestone

• [x] Select Project
• [x] Select Milestone as the next official version or Backlog of Development Ideas

New Feature Checklist

• [x] Complete the issue definition above, including the Time Estimate and Funding source.
• [ ] Fork this repository or create a branch of develop. Branch name: feature_<Issue Number>/<Description>
• [ ] Complete the development and test your changes.
• [ ] Add/update log messages for easier debugging.
• [ ] Add/update tests.
• [ ] Add/update documentation.
• [ ] Push local changes to GitHub.
• [ ] Submit a pull request to merge into develop. Pull request: feature <Issue Number> <Description>
• [ ] Define the pull request metadata, as permissions allow. Select: Reviewer(s), Project, and Development issue Select: Milestone as the next official version
• [ ] Iterate until the reviewer(s) accept and merge your changes.
• [ ] Delete your fork or branch.
• [ ] Close this issue.

[hahnd]

@michelleharrold @fossell I got the GRIB2-to-GeoJson mostly working today. Couple of questions:

• When I read the variable data from the GRIB2 file, I get a 1-D array of floats back for the 2-D variables. For example, I get an array of length 348,801, rather than a 2-D array with dimensions 499x699. I have code to restructure this as a 2-D float array so that the contouring code still works. Where can I find the original grid dimensions? I currently have this hard-coded to 499x699, but that only works for our test domain.
• How to the variables in the output products table map to the variable names in the GRIB2 file?

[michelleharrold]

@hahnd - wow, that is great news! For your questions:

• In the WRF namelist.input, e_we = 700 and e_sn = 500, so it looks like e_we - 1 and e_sn - 1, but I will defer to @fossell if there is a better way to get around hard-coding.
• First thing tomorrow, I will add a column in the output products that tells what the UPP output variable name is. I should have it done by our meeting tomorrow AM.

[fossell]

@JohnHalleyGotway - recommends using an existing package that already handles grib2 data instead of reinventing the wheel. Grib2 packing can be tricky. John will check with MET folks and get some advice on doing this.

[JohnHalleyGotway]

@bikegeek and @georgemccabe, I'm writing to request some advice about reading GRIB2 data via python for this WRF Cloud project. For METplus, I see that we use pygrib, as shown in this list of testing environments.

I warned @hahnd that reading GRIB2 data can become a rat's nest and recommended that we use pygrib to do so. But he ran into headaches with pygrib's dependency on eccodes... concerns that were echoed by @michelleharrold.

I'm wondering about your perspectives about pygrib. Is working through the headaches of eccodes worth the effort? If you had it to do over again, would you make any different decisions related to this?

[georgemccabe]

I have used pygrib in an IFI script to concatenate records to create 3D radar files and fix inconsistent timing information. There is also at least 1 METplus use case that uses it. I vaguely remember some issues with eccodes but I don't remember the details. I haven't experienced any issues lately, so I can't speak to working through the headaches.

It looks like we install pygrib using conda for the use case testing which appears to handle the eccodes dependency automatically. If you can use conda and have conflicting eccodes version dependencies, you could create a separate conda environment that contains pygrib and isolate it from the other tools that require eccodes.

It just checked and pygrib is installed on RAL machines under /usr/local/python3, but when I imported it, I did get an eccodes warning.

/usr/local/python3/bin/python3

Python 3.8.3 | packaged by conda-forge | (default, Jun 1 2020, 17:43:00) [GCC 7.5.0] on linux Type "help", "copyright", "credits" or "license" for more information.

>>> import pygrib

:219: UserWarning: Warning: ecCodes 2.19.1 or higher is recommended. You are running 2.17.0.

There are conda environments on seneca (RAL machine) that have pygrib installed and do not report any warning when you import it.

[hahnd]

I am currently using xarray+cfgrib to read the grib2 file. It is working fine except I cannot figure out how to get the original grid dimensions. All of the variables have a 1-D shape with a length of 348,801, rather than a 2-D shape with dimensions of 699x499. Do either of you have any advice on how to get the dimensions via xarray+cfgrib, or know that pygrib+eccodes is easier to work with?

[michelleharrold]

@hahnd - I am not sure I will be more help than @georgemccabe or @bikegeek, but I can at least offer examples of what I have used with post-created grib2 output (take or leave!):

For xarray+cfgrib, here are two example scripts I have used (this was a couple years ago, so I might be rusty at remembering the exact details!).

For pygrib+eccodes, we installed pygrib with conda, so as @georgemccabe mentioned, the dependencies were handled automatocally (dockerfile for reference). Example of plotting 2-m T is here.

[fossell]

An alternative, although I don't think most efficient, option we talked about yesterday was using wgrib2 utility to convert the grib2 files to netcdf, and then use the netcdf files directly. I created a test file and David was able to use that and generate plots. But it's an extra conversion step and will inflate the disk usage too, so probably not the greatest long term solution.

[georgemccabe]

@hahnd, the script that is used in the METplus use case uses cfgrib to read data and it preserves the shape of the data. I tried running the script line-by-line in the Python interpreter and confirmed that it is read as a 3D dataset.

>>> ds['t'].shape
(21, 361, 720)

Here is the script: https://github.com/dtcenter/METplus/blob/develop/parm/use_cases/model_applications/medium_range/TCStat_SeriesAnalysis_fcstGFS_obsGFS_FeatureRelative_SeriesByLead_PyEmbed_Multiple_Diagnostics/gfs_pv_analysis.py#L33

I didn't write this script, but if you'd like to send me the file you are using I can try to take a look if you can't borrow logic from the script to accomplish what you need.

[hahnd]

Thanks for the examples! @michelleharrold - it was not clear from the code how the plotting library is reading the dimensions. @georgemccabe - when I look at the shape of the variable, it is a 1D array.

Yesterday afternoon, I switched over to using pygrib and got the code working fine. If I have to switch back to xarray+cfgrib for any reason, then I will provide a sample file for @georgemccabe and see if you can get the grid dimensions. Thanks!

[fossell]

@hahnd - If pygrib is working, that seems like maybe the best solution moving forward. I had put adding wgrib2 to the AMI on my list of things to do, but if you don't think there's a need at this point, I'll table that.

[fossell]

Closed with #93