Working with GSI Diagnostics files

[esheehan-gsl]

@ian-noaa I’ve been looking into these diagnostic files produced by GSI. I haven’t had much luck parsing them directly with Python. It’s possible—although I haven’t been able to confirm this—that this binary file format is actually a Fortran-specific format for serializing Fortran data structures; similar to Python’s pickle format. I’m beginning to think our best bet for working with these files is to use Fortran to read them.

I have two ideas about how to build a pipeline that can convert these files into something we can work with:

1. Use utilities from GSI to build the pipeline
2. Write our own parser in Fortran and wrap it for use in Python

GSI Utilities

I found a python_utilities/ folder in the GSI repo, which has a lib_GSI.py module. lib_GSI.py appears to have utilities for reading from the diag files. It calls into a read_diag.f90 module, which must be wrapped for use in Python.

If we can figure out how to get read_diag and lib_GSI installed, we can just write a little CLI app in Python that imports lib_GSI and do whatever we need to do with the data from the diag files.

What I like about this solution is that it helps us keep our parser in sync with GSI. What I don’t like is that I suspect we’d have to build all of GSI just to get to these utilities, which seems like a lot of extra work that we’re not even really going to use.

@guoqing-noaa do you know anything about installing and using these Python utilities?

Custom Fortran Parser

Alternatively, we could implement our own parser based on what we see in read_diag.f90 and read_diag_conv.f90 and we can use NumPy’s f2py module to create a wrapper that enables us to call into the Fortran module from Python.

What I like about this solution is that we can build the parser so that it does exactly what we need and not have to worry about the build process for GSI. What I don’t like about it is that we’d have to keep our eye on GSI for any changes to the format and update our parser.

---

What do you think?

[JacobCarley-NOAA]

Speaking from experience, read_diag_conv.f90 might be the way to go. I've used this extensively in the past, but just made an intermediate ascii file and then read that with Pandas instead of having a python wrapper. I doubt you will need to be too concerned with GSI changes here, as development with GSI, while still active, is waning.

However - we may well switch to the netcdf-based diag file format from GSI in the semi-near future. And that would make this all OBE - perhaps we ought to just make that transition sooner rather than later. With netcdf there are plenty of nice tools out there (like PyGSI). For what it's worth, I cannot see a reason why we should wait to make the switch over. Perhaps that's a discussion for an upcoming 3DRTMA and RRFS devs meeting.

[JacobCarley-NOAA]

tagging a few folks here who can also weigh in @ManuelPondeca-NOAA @ShunLiu-NOAA @hu5970

[ian-noaa]

Agreed. NetCDF would be great. In that case, we wouldn't need to deal with GSI or with a custom parser. Any idea how "semi-near" that future could be?

If GSI is relatively stable, I could see just wrapping read_diag.f90 and read_diag_conv.f90. NumPy's f2py looks relatively easy to use and integrate into the CI process if needed. I'd bet that this would be easier to get going (and understand) than trying to figure out how to get GSI installed in a container.

[guoqing-noaa]

I agree that NetCDF diag would be a better option. I can try to turn on generating NetCDF diag files.

[esheehan-gsl]

@guoqing-noaa if you produce some NetCDF diagnostic files for us to test, you can upload them to the Sample Data folder in Google Drive. 🙏🏻

[guoqing-noaa]

@esheehan-gsl I've upload sample NetCDF diag files. Please check the ncdiag.tar file in the google folder.

[ManuelPondeca-NOAA]

I agree, NetCDF diag files seems to be the way to go.