Skip to content

Commit

Permalink
Merge pull request #592 from sfinkens/fix_531
Browse files Browse the repository at this point in the history 
Fix masking of AHI HSD space pixels
  • Loading branch information
@djhoese
djhoese authored Jan 25, 2019
2 parents e473030 + 9aec74e commit 5c5f0c5
Show file tree
Hide file tree
Showing 2 changed files with 131 additions and 113 deletions.
218 changes: 105 additions & 113 deletions satpy/readers/ahi_hsd.py
View file
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
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')
@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

0 comments on commit 5c5f0c5

Please sign in to comment.