Implement a data model

[pulkkins]

Pysteps operates on 2d grids that contain metadata (such as the grid definition, units, applied transformations and temporal resolution). Currently this information is carried in a separate metadata dictionary or as separate arguments (see e.g. pysteps.nowcasts.steps.forecast) in a very inconsistent way. It would be a good idea to store all this information into the same object.

Possible solution: use xarray (http://xarray.pydata.org/en/stable) that allows storing metadata into attributes, as well as using named columns and integration with dask.

[dnerini]

This is definitely needed and an excellent idea!

The use of xarray can indeed represent an ideal solution to this issue, at the cost of introducing an additional external dependency. xarray also offers a natural interface to netCDF files, pysteps output format.

However, this is not an easy task and will require considerable work. @pulkkins how do you plan to introduce such data model?

[aperezhortal]

May I suggest another alternative? Instead of using an xarray object, a class can be used to handle a collection of radar composite images. Internally, this class can store the radar 2D data, metadata, the array with the observations quality and the grid information. The benefit over an xarray object is that more method can be added to the class (for example the import, export, or maybe data preprocessing methods, etc.). In case that xarray need to be used, the class can include a method to export part of the data to a DataArray object.

This is how the class could look like (with an execution example at the end):

import numpy

from pysteps.io import (import_fmi_pgm, import_bom_rf3,
                        import_mch_gif, import_odim_hdf5)


class RadarComposites(object):
    """
    Class to store a collection of 2D precipitation fields.
    """

    def __init__(self):
        self.times = None  # Not used for now, but this could be useful to retrieve 2D fields at a given times.
        self.data = None  # precip data
        self.quality = None
        self.latitude = None  # Not used for now
        self.longitude = None  # Not used for now
        self.metadata = None

    def add(self, file_paths, file_format='fmi', **kwargs):
        """ Add radar 2D composites to the collection."""

        file_format = file_format.lower()

        available_importers = dict(bom_rf3=import_bom_rf3,
                                   bom=import_bom_rf3,
                                   fmi_pgm=import_fmi_pgm,
                                   fmi=import_fmi_pgm,
                                   mch_gif=import_mch_gif,
                                   odim_hdf5=import_odim_hdf5)

        if file_format not in available_importers:
            raise ValueError("Format not supported.\n"
                             + "The supported formats are: "
                             + str(list(available_importers.keys())))

        importer = available_importers[file_format]
        data = list()
        quality = list()
        metadata = list()

        for file in file_paths:
            _data, _quality, _metadata = importer(file, **kwargs)

            metadata.append(_metadata)

            data.append(_data)
            quality.append(_quality)

        data = numpy.asarray(data)
        quality = numpy.asarray(quality)
        if self.data is None:
            self.data = data
        else:
            self.data = numpy.concatenate((self.data, data), axis=0)

        if self.quality is None:
            self.quality = quality
        else:
            self.quality = numpy.concatenate((self.quality, quality), axis=0)

        # For now check that the metadata is the same for each time
        equal_metadata = numpy.all([metadata[0] == value
                                    for value in metadata])

        if not equal_metadata:
            print("Warning: Metadata differs for between some files!")

        self.metadata = metadata[0]


if __name__ == '__main__':
    #Usage example
    from os.path import join
    import glob

    data_dir = "../data/fmi/20170509"

    files = glob.glob(join(data_dir, '*.pgm.gz'))

    my_composites = RadarComposites()
    #Add composites to the collection
    my_composites.add(files, file_format='fmi', gzipped=True)

Another idea in the same line, is to create a sub-class from the xarray.Dataset class with the functionalities that I mention before.

[dnerini]

So far we have tried to avoid using classes. The reason was simply that we wanted to keep the code as simple and linear as possible. But as you nicely suggest, classes can offer a way to manage data and metadata, eventually improving the workflow. It is worth considering the option @aperezhortal suggests.

[wcwoo]

Hi everybody! I'm from the Hong Kong Observatory, currently responsible for the Community SWIRLS nowcast system, which is in operational use in meteorological services in China, India, Malaysia and South Africa. We're also working with those in the Philippines, Vietnam, Myanmar and soon other Southeast Asia countries.

From our (somewhat painful) experience, I fully agree that it is necessary to adopt a single and consistent data model, and that xarray is a good choice. You might be interested to know that we've recently developed some algorithms to read in radar data files and returns xarray.DataArray objects. Currently supported formats include ESRI ASCII, HDF5/ODIM, IRIS (raw and reflectivity), NetCDF (three variants used by Philippines and Vietnam), and UF. The project website is at https://com-swirls.org/. Kindly let me know if you would like to know more.

[dnerini]

Hi @wcwoo, many thanks for your comment. We can certainly learn a lot from your experience with the SWIRLS system!
I think xarray still represents to us the best way to go and you are confirming us this impression.

I wasn't aware of the community version com-SWIRLS, I will certainly have a look at the project website that you've provided.

[wcwoo]

Sorry that I missed your reply, as I forgot that I setup an auto email sorting rule some time ago. 😅

I wrote this page to describe Com-SWIRLS' data model. The gist is to make development divisible and modules swappable. I hope one day we can run various QPF algorithms simultaneously to benchmark performance or to generate super-ensemble.

[kmuehlbauer]

Is there some solution for this ahead? I would fully support the xarray approach. There is also no need to subclass Dataset but xarray Accessors can be used instead.

This would greatly simplify interaction with other packages.

Just today I had to go the long way importing/exporting data/results to/from pysteps to xarray to conveniently examine it using eg hvplot and feed it to my further processing chain.

If you consider xarray useful I would try to setup a Pull Request or example code.

[kmuehlbauer]

I've created a simple gist, which shows some capabilities of xarray and how it could be used for pysteps.

If there is real interest here, I can try to setup a simple workflow for some LK method etc. which uses xarray.pipe technique.