Rewrite grib implementation

[brancomat]

Rewrite grib encoding:

• adding grib2 support (template 4.32 - suggestion 15/10/2013)
• adding the possibility to compute reflectivity values from grib

[spanezz]

grib_api keys seem to be documented here: http://www.cosmo-model.org/content/model/documentation/grib/pdtemplate_4.32.htm

[spanezz]

It looks like template 4.32 is for synthetic satellite data, and for satellite products it's template 4.31 (see http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table4-0.shtml ).

Unfortunately, I am having a hard time implementing this using the documentation in http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp4-31.shtml or the one in http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp4-32.shtml and I'm now looking for something more comprehensive.

Help is welcome.

[spanezz]

I'll now try to proceed by reading backwards from what is expected in grib_api definition files, and trying to cross-reference and guess from there.

[dcesari]

Official WMO documentation is here (linked from http://www.wmo.int/pages/prog/www/WMOCodes/WMO306_vI2/LatestVERSION/LatestVERSION.html grib2 templates); yes it is correct: templates 4.32-34 are for model-simulated satellite data, analog of 4.0, 4.1, 4.11 respectively. Real satellite data are written using template 4.31, AFAIK 4.30 is deprecated in favour of 4.31.

[spanezz]

Done in commit:09d7eea85da16c39a584890a91fa7520fe8e40ad.

The tests of conversions from SAFH5 currently fail: since SAF test data is for derived products that do not correspond to Meteosat channels, and use non-standard channel numbers, and since 4.31 product definition templates do not store a channel number but the channel central wave numbers, I have no way to map arbitrary SAF channels to central wave numbers and back. I do not know how to solve this: currently meteosatlib exits with an error if an image with a nonstandard satellite chamber is exported to GRIB2.

[spanezz]

Still missing: verify if and how reflectance data can also be encoded in grib, and used when processing.

[spanezz]

The reflectance of channels vis006, vis008, ir1.6 and hrv seems to be computed with the source channel only. The reflectance of channel ir0.39 however also needs channels ir1.08 and ir1.34, so we need a way to tell meteosatlib all the source files involved in the computation of reflectance. I'll now try to find out how it can be done.

[spanezz]

One option is that I introduce a new .msat file format, which is a text file that references various meteosatlib-openable data sources, and is opened as a multi-channel image. The resulting multichannel image will always have the 13 standard meteosat channels, and will automatically assign the input files to the right channel based on their channel number.

At that point, I can implement AddBand to request adding reflectance raster bands for given channels.

[spanezz]

Alternatively, I can remove the hackish reflectance virtual raster bands, and implement a command that just takes source files on the command line and outputs GDAL datasets with the reflectance values.

[spanezz]

At this stage I have a preference for the external tool, that would simplify the GDAL side of meteosatlib by taking reflectance calculations out of it. In the future we can spend some more time thinking how to integrate it with GDAL again, if we have a need to do so.

[brancomat]

The virtual raster bands are used in production so I'd like to keep them... our issue is the possibility to archive (portions of) satellite data in grib using arkimet mantaining the possibility of accessing that data in a fashion similar to python "products" script. All this, possibly without having impacts on the operational procedures. Also, a pony.

My 2 cents: if it's possible to mantain the possibility to compute the reflectance from vis006, vis008, ir1.6 and hrv grib data, it would be ok to compute ir 0.39 reflectance from raw data and archive it as a grib.
So maybe understanding if and how reflectance data can also be encoded in grib (cit) would be enough. But I'm not attached to this solution in particular as long as we can archive data subsets and generate products from the archive. We can negotiate about the pony.

[spanezz]

Ok, I'll now think about how it can be done. The pony would probably be the easiest part.

[spanezz]

The products script currently accesses satellite channels using the HRIT hackish file names, and as such is only able to open HRIT files. The trick with reflectance is that since the HRIT hackish file names are used to compute actual file names but do not exist in the file systems, they can be hacked even more to compute reflectance if an 'r' is added to the channel name part of them:

    def gdal_dataset(self, channel_name):
        fname = "H:MSG3:%s:%s" % (
                channel_name,
                self.dt.strftime("%Y%m%d%H%M")
        )

and

    @property
    def reflectance(self):
        if self.cached_refl is None:·
            self.cached_refl = self.sat.get_array(self.name + "r")
        return self.cached_refl

If we want to have products work using GRIB data, the way it generates the names of files to load as GDAL datasets needs to change.

[spanezz]

I've implemented reflectance calculations for all meteosatlib-supported dataset types (those in gdalinfo --formats|grep Msat).

To enable it, pass MSAT_COMPUTE=reflectance as a dataset open option, as in gdalwarp -oo MSAT_COMPUTE=reflectance. Open options can also be passed to gdal.OpenEx in python, but I still haven't checked how.

This works for all channels that have reflectance, except IR 3.9, because IR 3.9 cannot be computed from just one channel of information, and it also needs IR 10.8 and IR 13.4.

Computing reflectance for IR 3.9 using channel information from all sort of formats is actually implemented in the code, you can even use IR 3.9 data from GRIB, IR 10.8 from NetCDF and IR 13.4 from HRIT. However, there is currently no user interface to enable that.

Shall I create a new executable to compute reflectance for IR 3.9, and we see how it goes?

[brancomat]

I'm still behind in testing all this, but an executable to compute IR 3.9 reflectance sounds good

[spanezz]

Ho implementato il calcolo della riflettanza da dataset arbitrari via dataset GDAL virtuali.

In tests/data/reflectance_ir039.vrt c'è un esempio. I sorgenti nell'esempio sono HRIT, ma possono essere anche GRIB o qualunque altra cosa letta da meteosatlib.

Scrivere quel file .vrt non è al momento una cosa semplice. L'ordine dei sorgenti è importante, per esempio, e c'è da ripeterne uno 3 volte per calcolare 3 dati che servono per la riflettanza, perché quei dati non sono resi disponibili da GDAL alla PixelFunction msat_reflectance_ir039.

Ora ragiono di fare un programma che genera quel .vrt dati 3 dataset.

[spanezz]

In 885dfc0 ho aggiunto tools/msat-reflectance-ir39.

Ora, dati 3 file ir039.grib, ir108.nc, e ir134.tif, per avere la riflettanza si può fare:

msat-reflectance-ir39 ir039.grib ir108.nc ir134.tif > reflectance.vrt

E a quel punto reflectance.vrt è un file che quando viene aperto con GDAL calcola la riflettanza usando gli input con cui è stato generato.

msat-reflectance-ir39 è abbastanza furbo da trovare quali canali ci sono nei file in input a prescindere dall'ordine in cui vengono dati, e può leggere qualunque file supportato da GDAL a patto che contenga i metadati di meteosatlib.

Ti passo la palla per fare il punto.

[brancomat]

Ok, quindi il giro completo diventerà:
i file hrit vengono convertiti in grib (magari riproiettati e ritagliati) e archiviati in arkimet, tramite qualche script (da creare) in estrazione saranno poi ricalcolate dai grib le riflettanze del caso per poter rigenerare i prodotti.
Come step 1 appena riesco a verificare che dopo riproiezione e ritaglio si riescono a rigenerare le riflettanze, aggiorno il bug