Error in pygac data zenith angle

[carloshorn]

Hello,

For our reprocessing activities at EUMETSAT, we are using your AVHR GAC data from ECMWF FS. Unfortunately, some files show an artifact on the zenith angle.
My college Kristina Petraityte has already been in contact with Stephan Finkensieper who advised us to raise a ticket here.

The Error:
Not all files affected, but for example in ECC_GAC_sunsatangles_metopb_99999_20130101T2257085Z_20130102T0039455Z.h5 on image2 at row index 10081, the values for all 409 pixels are about twice the expected value (~12000 instead of ~6000).
We tried to reproduce the file using pygac, but surprisingly in our case, the feature disappeared. This means that the software seems to be correct. However, I would recommend to add the TLE lines into the output, because we are not sure, if we are using exactly the same input as you did.

More details on our reproduction of ECC_GAC_sunsatangles_metopb_99999_20130101T2257085Z_20130102T0039455Z.h5
branch: develop
git hash: bdd31d0
TLE: '1 38771U 12049A 13001.30683292 -.00000109 00000-0 -29853-4 0 9993\n',
'2 38771 098.7267 063.7881 0001633 317.7328 042.3723 14.21473300 15006\n'
pyorbital version: 1.5.0
python: latest anaconda 2.7 64 bit installer

Many thanks in advance!

[mraspaud]

Hi,

Indeed, it could be so that the data was to blame. Have you looked at the quality flags in the data ?

[carloshorn]

Sorry for asking, but could you point me to the location of the quality flag. I could not find it in the hdf5 file. Sill, it would be interesting to understand why we could not reproduce the error.

[mraspaud]

you probably should have a corresponding l1c file for that scene I suppose, with the quality flags in it.

[carloshorn]

Hi again,

we had a closer look at the differences between the files. Please see the attached figures.
I am not very familiar with the content of the quality flag file, but the values in the data array are the same for the downloaded files and our reproduced ones. The panels in the quality flag figure show the arrays in qual_flags/data.

The more concerning figure is the one revealing the differences between the downloaded data and the reproduced one. It shows the comparison of data for an arbitrarily chosen pixel index (100). We already removed all data filled with the missing data wildcard (which is the same for the downloaded and the reproduced data). The panels in every second line shows the difference between our reproduction and the downloaded data of the above panels.
This figure clearly exhibits the erroneous value in the downloaded data for the satellite zenith angle which initially attracted our attention. The differences in the lon/lat values are at the rounding precision which therefore should have only little impact. However, now also concerning are relatively large differences for example in AVHRR ch3a.

The thing that puzzles us is that we cannot reproduce these outliers and that we find differences at all. The only thing that might be different in our input might be the TLE.

Maybe also worth mentioning, not everything cannot be reproduced. We find the exact agreement in AVHRR ch3b, AVHRR ch4, AVHRR ch5.

[mraspaud]

Ok, we talked about this with @abhaydd , and he pointed out that the version of the code used for CLARA-A2 and within the ESA-CCI project are what is now the 'feature-clock' branch. It seems that the files that you are using matches is from one of these projects, so you could try to change the version of pygac. The latest updates, including updating the calibration routines is available on the 'develop' branch. The master branch hasn't been updated yet as some features are still work in progress (like merging lac and gac reading). So please try the 'feature-clock' branch and give us an update :)

[carloshorn]

Hi @mraspaud and @adybbroe,

partially good news. Using the 'feature-clock' branch (55d2107) gives correct results for the different channels, but still we do not reproduce the features in the angles.

Looking at the code GACReader.get_angles (gac_reader.py line 313), we see that these angles are calculated using the external library 'pyorbital'. Likely, there was a bug in an older version. Do you still know which version you have been using?

[mraspaud]

I'll drag in @sfinkens in the conversation, he might have some ideas.

[sfinkens]

@mattescarlos Interesting analysis. I assume what you refer to as download is the level 1c dataset which CLARA-A2 and Cloud_cci v3 are based on (probably downloaded from ec:/sf7/data/AVHRR_GAC_L1c_archive_v3)? We used pyorbital-0.3.2 for CLARA-A2 and Cloud_cci v3. Hope this helps.

[carloshorn]

Hi @sfinkens,

yes, it is the data from ec:/sf7/data/AVHRR_GAC_L1c_archive_v3. We now also did the calculation with pyorbital-0.3.2, but still the feature in the angles remain. Since it is not a software problem, we cannot explain how this outlier in the data could occur. It is unlikely, but maybe there was an error introduced when downloading/copying the data. Could you maybe check if you also find this strange spike in the angles in your data?

From the git logs of pyorbital, I was very confident that this might explain the issue, but maybe it is an artifact of the specific TLE you have been using for processing the affected files, e.g. NSS.GHRR.M1.D13001.S2257.E0039.B0150910.SV.

Many thanks in advance!

[sfinkens]

@mattescarlos Sure, but it might take some time. One question though: Do you really need to reproduce the faulty behaviour? I think you'd be better off with the most recent version (or at least the feature-clock branch).

[carloshorn]

Hi @sfinkens,

Thank you for the swift reply. Using the most recent version might be the way to proceed.

Well, from the software engineering point of view, reproducing errors is crucial to check the deterministic behavior of software. Only understood issues can reliably be fixed. The quality of our products depends on that.

Initially, we wanted to use the already existing hdf5 files, but our reprocessing software started to raise exceptions. The subsequent back tracing finally brought us here.

[mraspaud]

@sfinkens Does someone in the CMSAF have a complete list of the software versions used for CLARA-A2 ?

[sfinkens]

Yes, there is an internal document listing all the versions. I cannot post it here, but I can paste the version numbers:

Main Modules
============
pyGAC v0.1.0-py2.7 (https://github.com/pytroll/pygac/releases/tag/ClaraA2)
PPS 2014-patch20150327
PPS_probcma 1.0
SAL 4.0
SIS 0.8
GAC_averaging 2.0

Dependencies
============
Proj 4.8.0
szip 2.1
HDF5 library 1.8.11
Hlhdf 0.8.1
Netcdf4 4.3.0
Python 2.7.5-01
Numpy 1.7
Jasper 1.900.1
RTTOV 11.2
Ahamap v2014-patch20150206
Acpg v2014-patch20150327
GRIB-API 1.11.0
pyproj 1.9.3-py2.7
pyresample 1.1.0-py2.7
pygrib 1.9.9-py2.7
configobj 5.0.5-py2.7
pyorbital V0.3.2-py2.7

[sfinkens]

Sorry, @mattescarlos, I cannot reproduce your results. I processed NSS.GHRR.M1.D13001.S2257.E0039.B0150910.SV with the above software versions. The generated ECC_GAC_sunsatangles_metopb_99999_20130101T2257085Z_20130102T0039455Z.h5 is identical to the one in our ECFS archive.

But I can reproduce similar plots if I don't exclude fill-values. Are you sure you mask fill-values before applying gain & offset? Here is an example:

import h5py
import matplotlib.pyplot as plt
import numpy as np


h5 = h5py.File('ECC_GAC_sunsatangles_metopb_99999_20130101T2257085Z_20130102T0039455Z.h5')

variables = {'solzen': '/image1',
             'satzen': '/image2',
             'relazi': '/image3',
             'solazi': '/image4',
             'satazi': '/image5'}

ncols = 3
nrows = 2
fig, axes = plt.subplots(nrows=nrows, ncols=ncols)

for i, (variable, path) in enumerate(variables.items()):
    # Read metadata
    what = h5.get(path + '/what')
    fv = what.attrs['missingdata']
    gain = what.attrs['gain']
    offset = what.attrs['offset']

    # Mask fill values and decompress data
    data = gain * np.ma.masked_equal(h5.get(path + '/data')[:], fv) + offset

    # Plot 100th pixel in all scanlines
    row, col = i // ncols, i % ncols
    ax = axes[row, col]
    ax.plot(data[:, 100])
    ax.set_ylabel(variable)

plt.show()

[carloshorn]

Hi @sfinkens,

Sorry for the confusion, the angle error is in the v2 data. The v3 data is fine. From our side, we can close the ticket.

Thank you!

[sfinkens]

@mattescarlos Glad you figured it out! You can close the issue now.