Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix masking of AHI HSD space pixels #592

Merged
merged 3 commits into from
Jan 25, 2019
Merged

Fix masking of AHI HSD space pixels #592

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d6824ae
Use geostationary mask from readers.utils
sfinkens Jan 25, 2019
1ad5fc3
Refactor read_band() to facilitate testing
sfinkens Jan 25, 2019
9aec74e
Fix Stickler issues
sfinkens Jan 25, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
218 changes: 105 additions & 113 deletions satpy/readers/ahi_hsd.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@
from pyresample import geometry
from satpy import CHUNK_SIZE
from satpy.readers.file_handlers import BaseFileHandler
from satpy.readers.utils import get_geostationary_angle_extent, np2str
from satpy.readers.utils import get_geostationary_mask, np2str

AHI_CHANNEL_NAMES = ("1", "2", "3", "4", "5",
"6", "7", "8", "9", "10",
Expand Down Expand Up @@ -317,126 +317,115 @@ def get_area_def(self, dsid):
self.area = area
return area

def geo_mask(self):
"""Masking the space pixels from geometry info."""
cfac = np.uint32(self.proj_info['CFAC'])
lfac = np.uint32(self.proj_info['LFAC'])
coff = np.float32(self.proj_info['COFF'])
loff = np.float32(self.proj_info['LOFF'])
nlines = int(self.data_info['number_of_lines'])
ncols = int(self.data_info['number_of_columns'])

# count starts at 1
local_coff = 1
local_loff = (self.total_segments - self.segment_number) * nlines + 1

xmax, ymax = get_geostationary_angle_extent(self.area)

pixel_cmax = np.rad2deg(xmax) * cfac * 1.0 / 2**16
pixel_lmax = np.rad2deg(ymax) * lfac * 1.0 / 2**16

def ellipse(line, col):
return ((line / pixel_lmax) ** 2) + ((col / pixel_cmax) ** 2) <= 1
def _read_header(self, fp_):
"""Read header"""
header = {}

cols_idx = da.arange(-(coff - local_coff),
ncols - (coff - local_coff),
dtype=np.float, chunks=CHUNK_SIZE)
lines_idx = da.arange(nlines - (loff - local_loff),
-(loff - local_loff),
-1,
dtype=np.float, chunks=CHUNK_SIZE)
return ellipse(lines_idx[:, None], cols_idx[None, :])
header['block1'] = np.fromfile(
fp_, dtype=_BASIC_INFO_TYPE, count=1)
header["block2"] = np.fromfile(fp_, dtype=_DATA_INFO_TYPE, count=1)
header["block3"] = np.fromfile(fp_, dtype=_PROJ_INFO_TYPE, count=1)
header["block4"] = np.fromfile(fp_, dtype=_NAV_INFO_TYPE, count=1)
header["block5"] = np.fromfile(fp_, dtype=_CAL_INFO_TYPE, count=1)
logger.debug("Band number = " +
str(header["block5"]['band_number'][0]))
logger.debug('Time_interval: %s - %s',
str(self.start_time), str(self.end_time))
band_number = header["block5"]['band_number'][0]
if band_number < 7:
cal = np.fromfile(fp_, dtype=_VISCAL_INFO_TYPE, count=1)
else:
cal = np.fromfile(fp_, dtype=_IRCAL_INFO_TYPE, count=1)

header['calibration'] = cal

header["block6"] = np.fromfile(
fp_, dtype=_INTER_CALIBRATION_INFO_TYPE, count=1)
header["block7"] = np.fromfile(
fp_, dtype=_SEGMENT_INFO_TYPE, count=1)
header["block8"] = np.fromfile(
fp_, dtype=_NAVIGATION_CORRECTION_INFO_TYPE, count=1)
# 8 The navigation corrections:
ncorrs = header["block8"]['numof_correction_info_data'][0]
dtype = np.dtype([
("line_number_after_rotation", "<u2"),
("shift_amount_for_column_direction", "f4"),
("shift_amount_for_line_direction", "f4"),
])
corrections = []
for i in range(ncorrs):
corrections.append(np.fromfile(fp_, dtype=dtype, count=1))
fp_.seek(40, 1)
header['navigation_corrections'] = corrections
header["block9"] = np.fromfile(fp_,
dtype=_OBS_TIME_INFO_TYPE,
count=1)
numobstimes = header["block9"]['number_of_observation_times'][0]

dtype = np.dtype([
("line_number", "<u2"),
("observation_time", "f8"),
])
lines_and_times = []
for i in range(numobstimes):
lines_and_times.append(np.fromfile(fp_,
dtype=dtype,
count=1))
header['observation_time_information'] = lines_and_times
fp_.seek(40, 1)

header["block10"] = np.fromfile(fp_,
dtype=_ERROR_INFO_TYPE,
count=1)
dtype = np.dtype([
("line_number", "<u2"),
("numof_error_pixels_per_line", "<u2"),
])
num_err_info_data = header["block10"][
'number_of_error_info_data'][0]
err_info_data = []
for i in range(num_err_info_data):
err_info_data.append(np.fromfile(fp_, dtype=dtype, count=1))
header['error_information_data'] = err_info_data
fp_.seek(40, 1)

np.fromfile(fp_, dtype=_SPARE_TYPE, count=1)

return header

def _read_data(self, fp_, header):
"""Read data block"""
nlines = int(header["block2"]['number_of_lines'][0])
ncols = int(header["block2"]['number_of_columns'][0])
return da.from_array(np.memmap(self.filename, offset=fp_.tell(),
dtype='<u2', shape=(nlines, ncols), mode='r'),
chunks=CHUNK_SIZE)

def _mask_invalid(self, data, header):
"""Mask invalid data"""
invalid = da.logical_or(data == header['block5']["count_value_outside_scan_pixels"][0],
data == header['block5']["count_value_error_pixels"][0])
return da.where(invalid, np.float32(np.nan), data)

def _mask_space(self, data):
"""Mask space pixels"""
return data.where(get_geostationary_mask(self.area))

def read_band(self, key, info):
"""Read the data."""
# Read data
tic = datetime.now()
header = {}
with open(self.filename, "rb") as fp_:

header['block1'] = np.fromfile(
fp_, dtype=_BASIC_INFO_TYPE, count=1)
header["block2"] = np.fromfile(fp_, dtype=_DATA_INFO_TYPE, count=1)
header["block3"] = np.fromfile(fp_, dtype=_PROJ_INFO_TYPE, count=1)
header["block4"] = np.fromfile(fp_, dtype=_NAV_INFO_TYPE, count=1)
header["block5"] = np.fromfile(fp_, dtype=_CAL_INFO_TYPE, count=1)
logger.debug("Band number = " +
str(header["block5"]['band_number'][0]))
logger.debug('Time_interval: %s - %s',
str(self.start_time), str(self.end_time))
band_number = header["block5"]['band_number'][0]
if band_number < 7:
cal = np.fromfile(fp_, dtype=_VISCAL_INFO_TYPE, count=1)
else:
cal = np.fromfile(fp_, dtype=_IRCAL_INFO_TYPE, count=1)

header['calibration'] = cal

header["block6"] = np.fromfile(
fp_, dtype=_INTER_CALIBRATION_INFO_TYPE, count=1)
header["block7"] = np.fromfile(
fp_, dtype=_SEGMENT_INFO_TYPE, count=1)
header["block8"] = np.fromfile(
fp_, dtype=_NAVIGATION_CORRECTION_INFO_TYPE, count=1)
# 8 The navigation corrections:
ncorrs = header["block8"]['numof_correction_info_data'][0]
dtype = np.dtype([
("line_number_after_rotation", "<u2"),
("shift_amount_for_column_direction", "f4"),
("shift_amount_for_line_direction", "f4"),
])
corrections = []
for i in range(ncorrs):
corrections.append(np.fromfile(fp_, dtype=dtype, count=1))
fp_.seek(40, 1)
header['navigation_corrections'] = corrections
header["block9"] = np.fromfile(fp_,
dtype=_OBS_TIME_INFO_TYPE,
count=1)
numobstimes = header["block9"]['number_of_observation_times'][0]

dtype = np.dtype([
("line_number", "<u2"),
("observation_time", "f8"),
])
lines_and_times = []
for i in range(numobstimes):
lines_and_times.append(np.fromfile(fp_,
dtype=dtype,
count=1))
header['observation_time_information'] = lines_and_times
fp_.seek(40, 1)

header["block10"] = np.fromfile(fp_,
dtype=_ERROR_INFO_TYPE,
count=1)
dtype = np.dtype([
("line_number", "<u2"),
("numof_error_pixels_per_line", "<u2"),
])
num_err_info_data = header["block10"][
'number_of_error_info_data'][0]
err_info_data = []
for i in range(num_err_info_data):
err_info_data.append(np.fromfile(fp_, dtype=dtype, count=1))
header['error_information_data'] = err_info_data
fp_.seek(40, 1)

np.fromfile(fp_, dtype=_SPARE_TYPE, count=1)

nlines = int(header["block2"]['number_of_lines'][0])
ncols = int(header["block2"]['number_of_columns'][0])

res = da.from_array(np.memmap(self.filename, offset=fp_.tell(),
dtype='<u2', shape=(nlines, ncols), mode='r'),
chunks=CHUNK_SIZE)

invalid = da.logical_or(res == header['block5']["count_value_outside_scan_pixels"][0],
res == header['block5']["count_value_error_pixels"][0])
res = da.where(invalid, np.float32(np.nan), res)
header = self._read_header(fp_)
res = self._read_data(fp_, header)
res = self._mask_invalid(data=res, header=header)
self._header = header

logger.debug("Reading time " + str(datetime.now() - tic))

# Calibrate
res = self.calibrate(res, key.calibration)

# Update metadata
new_info = dict(units=info['units'],
standard_name=info['standard_name'],
wavelength=info['wavelength'],
Expand All @@ -453,7 +442,10 @@ def read_band(self, key, info):
satellite_altitude=float(self.nav_info['distance_earth_center_to_satellite'] -
self.proj_info['earth_equatorial_radius']) * 1000)
res = xr.DataArray(res, attrs=new_info, dims=['y', 'x'])
res = res.where(self.geo_mask())

# Mask space pixels
res = self._mask_space(res)

return res

def calibrate(self, data, calibration):
Expand Down
26 changes: 26 additions & 0 deletions satpy/tests/reader_tests/test_ahi_hsd.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,9 @@
import numpy as np
import dask.array as da
from datetime import datetime
from pyresample.geometry import AreaDefinition
from satpy.readers.ahi_hsd import AHIHSDFileHandler
from satpy.readers.utils import get_geostationary_mask


class TestAHIHSDNavigation(unittest.TestCase):
Expand Down Expand Up @@ -211,6 +213,30 @@ def test_calibrate(self, *mocks):
refl = fh.calibrate(data=counts, calibration='reflectance')
self.assertTrue(np.allclose(refl, refl_exp))

@mock.patch('satpy.readers.ahi_hsd.AHIHSDFileHandler._read_header')
djhoese marked this conversation as resolved.
Show resolved Hide resolved
@mock.patch('satpy.readers.ahi_hsd.AHIHSDFileHandler._read_data')
@mock.patch('satpy.readers.ahi_hsd.AHIHSDFileHandler._mask_invalid')
@mock.patch('satpy.readers.ahi_hsd.AHIHSDFileHandler.calibrate')
def test_read_band(self, calibrate, *mocks):
# Test masking of space pixels
nrows = 25
ncols = 100
self.fh.area = AreaDefinition(name='test', area_id='test', proj_id='test',
proj_dict={'a': '6378137.0', 'b': '6356752.3',
'h': '35785863.0', 'lon_0': '140.7',
'proj': 'geos', 'units': 'm'},
x_size=ncols, y_size=nrows,
area_extent=[-5499999.901174725, -4399999.92093978,
5499999.901174725, -3299999.9407048346])
calibrate.return_value = np.ones((nrows, ncols))
m = mock.mock_open()
with mock.patch('satpy.readers.ahi_hsd.open', m, create=True):
im = self.fh.read_band(info=mock.MagicMock(), key=mock.MagicMock())
# Note: Within the earth's shape get_geostationary_mask() is True but the numpy.ma mask
# is False
mask = im.to_masked_array().mask
ref_mask = np.logical_not(get_geostationary_mask(self.fh.area).compute())
self.assertTrue(np.all(mask == ref_mask))

def suite():
"""The test suite for test_scene."""
Expand Down