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 !
The image looks a bit strange, the green channel should be more present. It's a daytime scene right?
@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