Skip to content

Commit

Permalink
Merge pull request #1368 from mraspaud/fix_eps_number_of_scanlines
Browse files Browse the repository at this point in the history 
Fix wrong number of scanlines in eps reader
  • Loading branch information
@mraspaud
mraspaud committed Sep 17, 2020
2 parents 84f0071 + 1571c9d commit aa65035
Show file tree
Hide file tree
Showing 2 changed files with 155 additions and 44 deletions.
61 changes: 28 additions & 33 deletions satpy/readers/eps_l1b.py
View file
Expand Up @@ -161,6 +161,7 @@ def _read_all(self):
self.scanlines = self['TOTAL_MDR']
if self.scanlines != len(self.sections[('mdr', 2)]):
logger.warning("Number of declared records doesn't match number of scanlines in the file.")
self.scanlines = len(self.sections[('mdr', 2)])
self.pixels = self["EARTH_VIEWS_PER_SCANLINE"]

def __getitem__(self, key):
Expand All @@ -186,20 +187,8 @@ def keys(self):
keys += val.dtype.fields.keys()
return keys

@delayed(nout=2, pure=True)
def _get_full_lonlats(self, lons, lats):
nav_sample_rate = self["NAV_SAMPLE_RATE"]
if nav_sample_rate == 20 and self.pixels == 2048:
from geotiepoints import metop20kmto1km
return metop20kmto1km(lons, lats)
else:
raise NotImplementedError("Lon/lat expansion not implemented for " +
"sample rate = " + str(nav_sample_rate) +
" and earth views = " +
str(self.pixels))

def get_full_lonlats(self):
"""Get the interpolated lons/lats."""
"""Get the interpolated longitudes and latitudes."""
if self.lons is not None and self.lats is not None:
return self.lons, self.lats

Expand All @@ -210,25 +199,39 @@ def get_full_lonlats(self):
raw_lons = np.hstack((self["EARTH_LOCATION_FIRST"][:, [1]],
self["EARTH_LOCATIONS"][:, :, 1],
self["EARTH_LOCATION_LAST"][:, [1]]))
self.lons, self.lats = self._get_full_lonlats(raw_lons, raw_lats)
self.lons = da.from_delayed(self.lons, dtype=self["EARTH_LOCATIONS"].dtype,
shape=(self.scanlines, self.pixels))
self.lats = da.from_delayed(self.lats, dtype=self["EARTH_LOCATIONS"].dtype,
shape=(self.scanlines, self.pixels))
self.lons, self.lats = self._interpolate(raw_lons, raw_lats)
return self.lons, self.lats

@delayed(nout=4, pure=True)
def _interpolate(self, lons_like, lats_like):
nav_sample_rate = self["NAV_SAMPLE_RATE"]
if nav_sample_rate == 20 and self.pixels == 2048:
lons_like_1km, lats_like_1km = self._interpolate_20km_to_1km(lons_like, lats_like)
lons_like_1km = da.from_delayed(lons_like_1km, dtype=lons_like.dtype,
shape=(self.scanlines, self.pixels))
lats_like_1km = da.from_delayed(lats_like_1km, dtype=lats_like.dtype,
shape=(self.scanlines, self.pixels))
return lons_like_1km, lats_like_1km
else:
raise NotImplementedError("Lon/lat and angle expansion not implemented for " +
"sample rate = " + str(nav_sample_rate) +
" and earth views = " +
str(self.pixels))

@delayed(nout=2, pure=True)
def _interpolate_20km_to_1km(self, lons, lats):
# Note: delayed will cast input dask-arrays to numpy arrays (needed by metop20kmto1km).
from geotiepoints import metop20kmto1km
return metop20kmto1km(lons, lats)

def _get_full_angles(self, solar_zenith, sat_zenith, solar_azimuth, sat_azimuth):

nav_sample_rate = self["NAV_SAMPLE_RATE"]
if nav_sample_rate == 20 and self.pixels == 2048:
from geotiepoints import metop20kmto1km
# Note: interpolation asumes lat values values between -90 and 90
# Note: interpolation assumes second array values between -90 and 90
# Solar and satellite zenith is between 0 and 180.
# Note: delayed will cast input dask-arrays to numpy arrays (needed by metop20kmto1km).
sun_azi, sun_zen = metop20kmto1km(solar_azimuth, solar_zenith - 90)
sun_azi, sun_zen = self._interpolate(solar_azimuth, solar_zenith - 90)
sun_zen += 90
sat_azi, sat_zen = metop20kmto1km(sat_azimuth, sat_zenith - 90)
sat_azi, sat_zen = self._interpolate(sat_azimuth, sat_zenith - 90)
sat_zen += 90
return sun_azi, sun_zen, sat_azi, sat_zen
else:
Expand All @@ -238,7 +241,7 @@ def _get_full_angles(self, solar_zenith, sat_zenith, solar_azimuth, sat_azimuth)
str(self.pixels))

def get_full_angles(self):
"""Get the interpolated lons/lats."""
"""Get the interpolated angles."""
if (self.sun_azi is not None and self.sun_zen is not None and
self.sat_azi is not None and self.sat_zen is not None):
return self.sun_azi, self.sun_zen, self.sat_azi, self.sat_zen
Expand All @@ -262,14 +265,6 @@ def get_full_angles(self):
sat_zenith,
solar_azimuth,
sat_azimuth)
self.sun_azi = da.from_delayed(self.sun_azi, dtype=self["ANGULAR_RELATIONS"].dtype,
shape=(self.scanlines, self.pixels))
self.sun_zen = da.from_delayed(self.sun_zen, dtype=self["ANGULAR_RELATIONS"].dtype,
shape=(self.scanlines, self.pixels))
self.sat_azi = da.from_delayed(self.sat_azi, dtype=self["ANGULAR_RELATIONS"].dtype,
shape=(self.scanlines, self.pixels))
self.sat_zen = da.from_delayed(self.sat_zen, dtype=self["ANGULAR_RELATIONS"].dtype,
shape=(self.scanlines, self.pixels))
return self.sun_azi, self.sun_zen, self.sat_azi, self.sat_zen

def get_bounding_box(self):
Expand Down
138 changes: 127 additions & 11 deletions satpy/tests/reader_tests/test_eps_l1b.py
View file
Expand Up @@ -24,10 +24,12 @@
from unittest import mock

import numpy as np
import xarray as xr
import pytest
import satpy
from satpy.tests.utils import make_dataid
import xarray as xr
from satpy.readers import eps_l1b as eps
from satpy.tests.utils import make_dataid

grh_dtype = np.dtype([("record_class", "|i1"),
("INSTRUMENT_GROUP", "|i1"),
("RECORD_SUBCLASS", "|i1"),
Expand Down Expand Up @@ -57,23 +59,38 @@ def create_sections(structure):
return sections


class TestEPSL1B(TestCase):
"""Test the filehandler."""
class BaseTestCaseEPSL1B(TestCase):
"""Base class for EPS l1b test case."""

def setUp(self):
"""Set up the tests."""
# ipr is not present in the xml format ?
def _create_structure(self):
structure = [(1, ('mphr', 0)), (1, ('sphr', 0)), (11, ('ipr', 0)),
(1, ('geadr', 1)), (1, ('geadr', 2)), (1, ('geadr', 3)),
(1, ('geadr', 4)), (1, ('geadr', 5)), (1, ('geadr', 6)),
(1, ('geadr', 7)), (1, ('giadr', 1)), (1, ('giadr', 2)),
(1, ('veadr', 1)), (1080, ('mdr', 2))]

(1, ('veadr', 1)), (self.scan_lines, ('mdr', 2))]
sections = create_sections(structure)
sections[('mphr', 0)]['TOTAL_MDR'] = b'TOTAL_MDR = 1080\n'
return sections


class TestEPSL1B(BaseTestCaseEPSL1B):
"""Test the filehandler."""

def setUp(self):
"""Set up the tests."""
# ipr is not present in the xml format ?
self.scan_lines = 1080
self.earth_views = 2048

sections = self._create_structure()
sections[('mphr', 0)]['TOTAL_MDR'] = (b'TOTAL_MDR = ' +
bytes(str(self.scan_lines), encoding='ascii') +
b'\n')
sections[('mphr', 0)]['SPACECRAFT_ID'] = b'SPACECRAFT_ID = M03\n'
sections[('mphr', 0)]['INSTRUMENT_ID'] = b'INSTRUMENT_ID = AVHR\n'
sections[('sphr', 0)]['EARTH_VIEWS_PER_SCANLINE'] = b'EARTH_VIEWS_PER_SCANLINE = 2048\n'
sections[('sphr', 0)]['EARTH_VIEWS_PER_SCANLINE'] = (b'EARTH_VIEWS_PER_SCANLINE = ' +
bytes(str(self.earth_views), encoding='ascii') +
b'\n')
sections[('sphr', 0)]['NAV_SAMPLE_RATE'] = b'NAV_SAMPLE_RATE = 20\n'

_fd, fname = mkstemp()
fd = open(_fd)
Expand Down Expand Up @@ -160,7 +177,106 @@ def mock_getitem(key):
sun_zen_np2 = np.array(avhrr_reader.sun_zen)
assert np.allclose(sun_zen_np1, sun_zen_np2)


class TestWrongScanlinesEPSL1B(BaseTestCaseEPSL1B):
"""Test the filehandler on a corrupt file."""

@pytest.fixture(autouse=True)
def inject_fixtures(self, caplog):
"""Inject caplog."""
self._caplog = caplog

def setUp(self):
"""Set up the tests."""
# ipr is not present in the xml format ?
self.scan_lines = 1080
self.earth_views = 2048

sections = self._create_structure()
sections[('mphr', 0)]['TOTAL_MDR'] = (b'TOTAL_MDR = ' +
bytes(str(self.scan_lines - 2), encoding='ascii') +
b'\n')
sections[('mphr', 0)]['SPACECRAFT_ID'] = b'SPACECRAFT_ID = M03\n'
sections[('mphr', 0)]['INSTRUMENT_ID'] = b'INSTRUMENT_ID = AVHR\n'
sections[('sphr', 0)]['EARTH_VIEWS_PER_SCANLINE'] = (b'EARTH_VIEWS_PER_SCANLINE = ' +
bytes(str(self.earth_views), encoding='ascii') +
b'\n')
sections[('sphr', 0)]['NAV_SAMPLE_RATE'] = b'NAV_SAMPLE_RATE = 20\n'
_fd, fname = mkstemp()
fd = open(_fd)

self.filename = fname
for _, arr in sections.items():
arr.tofile(fd)
fd.close()
self.fh = eps.EPSAVHRRFile(self.filename, {'start_time': 'now',
'end_time': 'later'}, {})

def test_read_all_return_right_number_of_scan_lines(self):
"""Test scanline assignment."""
self.fh._read_all()
assert self.fh.scanlines == self.scan_lines

def test_read_all_warns_about_scan_lines(self):
"""Test scanline assignment."""
self.fh._read_all()
assert "scanlines" in self._caplog.records[0].message
assert self._caplog.records[0].levelname == 'WARNING'

def test_read_all_assigns_int_scan_lines(self):
"""Test scanline assignment."""
self.fh._read_all()
assert isinstance(self.fh.scanlines, int)

def test_get_dataset_longitude_shape_is_right(self):
"""Test that the shape of longitude is 1080."""
key = make_dataid(name="longitude")
longitudes = self.fh.get_dataset(key, dict())
assert longitudes.shape == (self.scan_lines, self.earth_views)

def tearDown(self):
"""Tear down the tests."""
with suppress(OSError):
os.remove(self.filename)


class TestWrongSamplingEPSL1B(BaseTestCaseEPSL1B):
"""Test the filehandler on a corrupt file."""

@pytest.fixture(autouse=True)
def inject_fixtures(self, caplog):
"""Inject caplog."""
self._caplog = caplog

def setUp(self):
"""Set up the tests."""
self.scan_lines = 1080
self.earth_views = 2048
self.sample_rate = 23
sections = self._create_structure()
sections[('mphr', 0)]['TOTAL_MDR'] = (b'TOTAL_MDR = ' +
bytes(str(self.scan_lines), encoding='ascii') +
b'\n')
sections[('mphr', 0)]['SPACECRAFT_ID'] = b'SPACECRAFT_ID = M03\n'
sections[('mphr', 0)]['INSTRUMENT_ID'] = b'INSTRUMENT_ID = AVHR\n'
sections[('sphr', 0)]['EARTH_VIEWS_PER_SCANLINE'] = (b'EARTH_VIEWS_PER_SCANLINE = ' +
bytes(str(self.earth_views), encoding='ascii') +
b'\n')
sections[('sphr', 0)]['NAV_SAMPLE_RATE'] = (b'NAV_SAMPLE_RATE = ' +
bytes(str(self.sample_rate), encoding='ascii') +
b'\n')
_fd, fname = mkstemp()
fd = open(_fd)

self.filename = fname
for _, arr in sections.items():
arr.tofile(fd)
fd.close()
self.fh = eps.EPSAVHRRFile(self.filename, {'start_time': 'now',
'end_time': 'later'}, {})

def test_get_dataset_fails_because_of_wrong_sample_rate(self):
"""Test that lons fail to be interpolate."""
key = make_dataid(name="longitude")
with pytest.raises(NotImplementedError):
self.fh.get_dataset(key, dict())

0 comments on commit aa65035

Please sign in to comment.