pytroll/pyspectral

GK2A AMI RSR wavelengths are reversed

djhoese opened this issue · 4 comments

Code Sample, a minimal, complete, and verifiable piece of code

# requires AMI PR from Satpy

from satpy import Scene
from glob import glob
scn = Scene(reader='ami_l1b', filenames=glob('/data/satellite/ami/20190930_0300/*.nc'))
scn.load(['day_microphysics_eum'])

Problem description

AHI is able to make the day_microphysics_eum composite, but AMI is not. Tracking it down I found that the nir_reflectance modifier in Satpy was failing during some calls to pyspectral. Specifically the piece of the solar.py code that uses the wvl from the RSR files (I think). If I print out the values of wvl for AMI I see:

[3.9797828 3.979624  3.979466  ... 3.6831057 3.68297   3.6828344] 

This shows that the values are decreasing instead of increasing. The code is failing (see traceback below) when we try to tell np.linspace to give us -593 elements in our array. This number is calculated by doing end - start where end and start are the last and first wavelength in the above array. This shows that AMI's RSR data is stored backwards (90% sure).

If I add this right before the linspace call then it works:

        if start > end:
            start, end = end, start
            wvl = wvl[::-1]

But this is an ugly hack and I feel like pyspectral should probably require that RSR files have increasing wavelengths in its file.

Expected Output

No errors and some new data returned.

Actual Result, Traceback if applicable

Traceback

~/repos/git/satpy/satpy/composites/__init__.py in _init_refl3x(self, projectables)
    613             raise ImportError("No module named pyspectral.near_infrared_reflectance")
    614         _nir, _tb11 = projectables
--> 615         self._refl3x = Calculator(_nir.attrs['platform_name'], _nir.attrs['sensor'], _nir.attrs['name'])
    616 
    617     def _get_reflectance(self, projectables, optional_datasets):

~/repos/git/pyspectral/pyspectral/near_infrared_reflectance.py in __init__(self, platform_name, instrument, band, **kwargs)
     87         self.solar_flux = kwargs.get('solar_flux', None)
     88         if self.solar_flux is None:
---> 89             self._get_solarflux()
     90 
     91         self._rad3x = None

~/repos/git/pyspectral/pyspectral/near_infrared_reflectance.py in _get_solarflux(self)
    162                                     dlambda=0.0005,
    163                                     wavespace=self.wavespace)
--> 164         self.solar_flux = solar_spectrum.inband_solarflux(self.rsr[self.bandname])
    165 
    166     def emissive_part_3x(self, tb=True):

~/repos/git/pyspectral/pyspectral/solar.py in inband_solarflux(self, rsr, scale, **options)
    119         spectral response valid for an earth-sun distance of one AU.
    120         """
--> 121         return self._band_calculations(rsr, True, scale, **options)
    122 
    123     def inband_solarirradiance(self, rsr, scale=1.0, **options):

~/repos/git/pyspectral/pyspectral/solar.py in _band_calculations(self, rsr, flux, scale, **options)
    166         LOG.debug("Begin and end wavelength/wavenumber: %f %f ", start, end)
    167         dlambda = self._dlambda
--> 168         xspl = np.linspace(start, end, round((end - start) / self._dlambda) + 1)
    169 
    170         ius = InterpolatedUnivariateSpline(wvl, resp)

<__array_function__ internals> in linspace(*args, **kwargs)

~/miniconda3/envs/satpy_py37/lib/python3.7/site-packages/numpy/core/function_base.py in linspace(start, stop, num, endpoint, retstep, dtype, axis)
    128     num = _index_deprecate(num)
    129     if num < 0:
--> 130         raise ValueError("Number of samples, %s, must be non-negative." % num)
    131     div = (num - 1) if endpoint else num
    132 

ValueError: Number of samples, -593, must be non-negative.

Versions of Python, package at hand and relevant dependencies

Current master of everything.

Thank you for reporting an issue !

With my above fix I get:

image

Does that look right-ish? @simonrp84 @adybbroe @mraspaud

The image looks a bit strange, the green channel should be more present. It's a daytime scene right?

Yeah, I agree with @mraspaud missing the green color.
And I agree @djhoese that pyspectral should require that RSR files have increasing wavelengths. It should be fixed in the ami specific reader. Let me see if I can find the original RSR data...

@djhoese I have reproduced your error above. I didn't look at trying to generate the above image with the fix you tried in Satpy.
But with the above PR I get this image with the following piece of code:

from satpy import Scene
from satpy.utils import debug_on
from glob import glob
debug_on()
composite_name = 'day_microphysics_eum'
scn = Scene(reader='ami_l1b', filenames=GK2A_FILES)
scn.load([composite_name])
new_scn = scn.resample(scn.min_area(), resampler='native')
new_scn.show(composite_name)

GK2A_FILES is a list of filename urls for one slot/scene

gk2a_day_microphysics_eum_20201204_0520_thumb