Representation of NetCDF data

[sandorkertesz]

This discussion is about the representation of NetCDF data in earthkit-data.

Currently, a NetCDF file is handled in the following way when loaded with from_source("file", path, ...):

• the file is loaded it into an xarray dataset
• scanned for fields
• if fields are found the resulting object will be a NetCDFFieldList
• otherwise a NetCDFReader will be generated. Then we can call to_xarray() on it to convert it into a format we can actually work with.

There are various problems with the approach:

• not all NetCDF data can be/should be represented as a fieldlist. E.g. observations, point data etc.
• Improve initial NetCDF fieldlist scanning speed #198 demonstrates that the current field scanning can be very slow for non field data. The data in the issue is basically a time series but in the end (after > 20 s) it is identified as field data with fields containing a single point only!

The questions we need to address:

• What other representations we could use for non-field NetCDF data?
• How to avoid the (potentially costly) initial scan to figure out what type of data we have?

[EddyCMWF]

Hi @sandorkertesz , thanks for starting the dicussion thread. I want to say that this is could be a higher discussion than just NetCDF representation, I think the conversation is also relevant for formats such as coverageJSON which can contain gridded data, and therefore could also be represented as fieldlist. So, if the argument for converting things to fieldlist is consistency then we should also be trying to convert coverageJSON to fieldlists.

I do not think representing things as fieldlist automatically is a good idea, I think it should be on the user/applcation/downstrem-earthkit-method to decide what the representation should be. Also, as you point out, the behaviour in response to netDF can be different depending on whether it is a fieldlist or not, this already breaks our consistency.

If the consistency is along the lines of "gridded" vs "point/observations", I think this is a very messy (if-ey) game to code, and will not necessarily benefit users. I think we can expect a user to know if they are handling observation or gridded data, therefore can choose the visual representation/method of iteration. Similarly, the downstream code we write should also be able to make these decisions.

[JamesVarndell]

Thanks @sandorkertesz for opening up this discussion.

To play devil's advocate, do we need a default data representation for NetCDF - or any source type for that matter? Would we make things easier by simply not offering a default data representation at all in earthkit-data, or at least only offering something very basic at the top Data object level?

I don't think we can provide a one-size-fits-all data representation that will work for all of our users all of the time, but we can offer them the flexibility to choose the data representation that works for them. I know we already do this with all our to_<data_representation> methods, but by attempting to offer a default choice for the user, we increase the likelihood that we might make bad assumptions about the data, or about what the user actually wants to do with the data - while also creating a lot of work for ourselves to write the logic which works out what kind of data we're trying to represent (i.e. fields vs obs etc.).

Using NetCDF as an example, we know that the community tool of choice for working with NetCDF in Python is xarray, so a lot of users will want to call to_xarray() on their Data. We don't need to make any assumptions about the input data to provide the xarray representation the user has asked for, because we can offload all of that work to xarray. At the moment we're occasionally creating problems when reading NetCDF because we always try to scan the data for fields first, which as @sandorkertesz has pointed out can often be very slow, and isn't a required step unless a FieldList has been specifically requested.

With this in mind, as suggested by @EddyCMWF, I think it's worth asking this question at a higher level - i.e. what should be the representation of data in earthkit-data? Should the from_source() method do anything beyond identifying the source type? Would it be better to delay any attempt to find a representation of the data until the point at which the user or downstream software has asked for something specific?

[corentincarton]

Hi @sandorkertesz, @EddyCMWF and @JamesVarndell! Just want to add my point of view here but I don't know that well the internals of earthkit-data so correct me if I'm wrong.
I feel that the concept of Field we have in earthkit is too tied to a GRIB message, and that's why it completely breaks for NetCDF, because NetCDF is defined with arbitrary dimensions and coordinates, while GRIB follows a strict structure with a spatial grid definition.
My naive definition of a field follows more the NetCDF definition of a dataset than the definition of GRIB message. For instance, a 1D spatial array (slice or cloud of points), with a time dimension should be a Field, and it's only when you try to concatenate different instances of these Fields that they become a FieldList. It shouldn't be tied to spatial coordinates or a grid at that point.
Aren't we missing one level of abstraction here? Shouldn't we have a Field representing an arbitrary nd-array with attached metadata (no concept of grid) and a SpatialField something that recognises Lat/Lon or other spatial coordinates? When reading a NetCDF file, we could choose the right level of abstraction following the dimensions that we have in the dataset. And a list of NetCDF file would become a FieldList of those objects, but one field could contain one year of spatial data for instance, or a full ensemble! For GRIB, it would directly take the SpatialField (SpatialFieldList) direction.
I think it's important that we define well the top level of abstraction, which has to match all the requirements from the data formats we have (cloud of points, gridded data, one-point observations with time, etc.), and then we can create the proper child classes from there. Do we have that level of abstraction at the moment? Looks like we're missing if from what I see there (Base is an abstract class): https://github.com/ecmwf/earthkit-data/blob/develop/earthkit/data/core/fieldlist.py#L20

[pmaciel]

I agree with almost everything here, but I'll point out what I believe is significant but left out of discussino so far.

I don't see why GRIBs being stricter than netCDF is a problem. I can imagine we could extend the concept of "gridspec" for describing gridded data, coming from any format. Another significant "non-spatial" case is spectral data. Observations are generally (?) spacial data, unstructured in time and space -- in physics they are "events": describing a time and a place.

It might be useful to have a hierarchy starting from Field (base), then SpatialField, SpectralField, ..., (still abstract), and FieldLists can collect base Field's of any type (just like a GRIB file, a bit different on netCDF ass they share mostly the same descriptive metadata in general). Since metadata is easy to transport, it can be attached to every single field in a FieldList, hence all formats I know of are taken care of.

[corentincarton]

Fully agree with you @pmaciel, especially on the SpatialField, SpectralField, etc. and the collection of Field of any type.
I think the current efficiency issue with NetCDF is that we split a NetCDF dataset into a FieldList to match the definition of a Field, which is currently defined as an array on a grid (which could be 1D, I agree) but nothing else (no time, no ensemble, etc.). If I'm not mistaken, we duplicate all the NetCDF metadata per Field when we do that, which makes it very inefficient for long time series. I think a Field should be able to contain more than a grid, it should be able to host an arbitrary number of dimensions (and that links to the reshape/flatten discussion we had the other day when discussing the interface for xarray).
And to challenge the concept of a grid, how would we treat spatial average? A time series of global temperature on the globe for instance.