From 6541d2d9a09b2d941c054b47ec5dde50471d7bf5 Mon Sep 17 00:00:00 2001 From: Panu Lahtinen Date: Wed, 12 Aug 2015 13:56:54 +0300 Subject: [PATCH] Support a list of files which will be concatenated --- mpop/satin/aapp1b.py | 235 ++++++++++++++++++++++++++----------------- 1 file changed, 142 insertions(+), 93 deletions(-) diff --git a/mpop/satin/aapp1b.py b/mpop/satin/aapp1b.py index c8166f4a..56d050a0 100644 --- a/mpop/satin/aapp1b.py +++ b/mpop/satin/aapp1b.py @@ -425,7 +425,9 @@ class AAPP1b(object): """ def __init__(self, fname): - self.filename = fname + if not isinstance(fname, (list, tuple, set)): + fname = [fname] + self.filenames = fname self.channels = dict([(i, None) for i in AVHRR_CHANNEL_NAMES]) self.units = dict([(i, 'counts') for i in AVHRR_CHANNEL_NAMES]) @@ -439,24 +441,34 @@ def read(self): """Read the data. """ tic = datetime.datetime.now() - with open(self.filename, "rb") as fp_: - header = np.memmap(fp_, dtype=_HEADERTYPE, mode="r", - shape=(_HEADERTYPE.itemsize, )) - data = np.memmap(fp_, - dtype=_SCANTYPE, - offset=22016, mode="r") + all_data = [] + all_headers = [] + for fname in self.filenames: + with open(fname, "rb") as fp_: + header = np.memmap(fp_, dtype=_HEADERTYPE, mode="r", + shape=(_HEADERTYPE.itemsize, )) + all_headers.append(header) + data = np.memmap(fp_, + dtype=_SCANTYPE, + offset=22016, mode="r") + all_data.append(data) LOGGER.debug("Reading time %s", str(datetime.datetime.now() - tic)) - self._header = header - self._data = data + self._header = all_headers + self._data = all_data def navigate(self): """Return the longitudes and latitudes of the scene. """ tic = datetime.datetime.now() - lons40km = self._data["pos"][:, :, 1] * 1e-4 - lats40km = self._data["pos"][:, :, 0] * 1e-4 + lons40km, lats40km = [], [] + for data in self._data: + lons40km.append(data["pos"][:, :, 1] * 1e-4) + lats40km.append(data["pos"][:, :, 0] * 1e-4) + + lons40km = np.vstack(lons40km) + lats40km = np.vstack(lats40km) try: from geotiepoints import SatelliteInterpolator @@ -489,83 +501,115 @@ def calibrate(self, chns=("1", "2", "3A", "3B", "4", "5"), """ tic = datetime.datetime.now() - if "1" in chns: - if calib_coeffs is not None: - coeffs = calib_coeffs['ch1'] - else: - coeffs = None - self.channels['1'] = _vis_calibrate(self._data, 0, - calibrate, pre_launch_coeffs, - coeffs) - self.units['1'] = '%' - if calibrate == 0: - self.units['1'] = '' - - if "2" in chns: - if calib_coeffs is not None: - coeffs = calib_coeffs['ch2'] - else: - coeffs = None - self.channels['2'] = _vis_calibrate(self._data, 1, - calibrate, pre_launch_coeffs, - coeffs) - self.units['2'] = '%' - if calibrate == 0: - self.units['2'] = '' + channels = {'1': [], '2': [], '3A': [], '3B': [], '4': [], '5': []} + is3b_all = [] - if "3A" in chns or "3B" in chns: - # Is it 3A or 3B: - is3b = np.expand_dims(np.bitwise_and( - np.right_shift(self._data['scnlinbit'], 0), 1) == 1, 1) - self._is3b = is3b + i = 0 + for data in self._data: - if "3A" in chns: - if calib_coeffs is not None: - coeffs = calib_coeffs['ch3a'] - else: - coeffs = None - ch3a = _vis_calibrate(self._data, 2, - calibrate, pre_launch_coeffs, - coeffs) - self.channels['3A'] = np.ma.masked_array(ch3a, is3b * ch3a) - - self.units['3A'] = '%' - if calibrate == 0: - self.units['3A'] = '' + if "1" in chns: + LOGGER.debug("Calibrating channel 1.") + if calib_coeffs is not None: + coeffs = calib_coeffs['ch1'] + else: + coeffs = None + channels['1'].append( _vis_calibrate(data, 0, + calibrate, + pre_launch_coeffs, + coeffs)) + self.units['1'] = '%' + if calibrate == 0: + self.units['1'] = '' + + if "2" in chns: + LOGGER.debug("Calibrating channel 2.") + if calib_coeffs is not None: + coeffs = calib_coeffs['ch2'] + else: + coeffs = None + channels['2'].append(_vis_calibrate(data, 1, + calibrate, + pre_launch_coeffs, + coeffs)) + self.units['2'] = '%' + if calibrate == 0: + self.units['2'] = '' + + if "3A" in chns or "3B" in chns: + # Is it 3A or 3B: + is3b = np.expand_dims(np.bitwise_and( + np.right_shift(data['scnlinbit'], 0), 1) == 1, 1) + # self._is3b = is3b + is3b_all.append(is3b) + + if "3A" in chns: + LOGGER.debug("Calibrating channel 3a.") + if calib_coeffs is not None: + coeffs = calib_coeffs['ch3a'] + else: + coeffs = None + ch3a = _vis_calibrate(data, 2, + calibrate, pre_launch_coeffs, + coeffs) + channels['3A'].append(np.ma.masked_array(ch3a, is3b * ch3a)) + + self.units['3A'] = '%' + if calibrate == 0: + self.units['3A'] = '' + + if "3B" in chns: + LOGGER.debug("Calibrating channel 3b.") + ch3b = _ir_calibrate(self._header[i], data, 0, calibrate) + channels['3B'].append( + np.ma.masked_array(ch3b, + np.logical_or((is3b is False) * \ + ch3b, + ch3b < 0.1))) + if calibrate == 1: + self.units['3B'] = 'K' + elif calibrate == 2: + self.units['3B'] = 'W*m-2*sr-1*cm ?' + else: + self.units['3B'] = '' + + if "4" in chns: + LOGGER.debug("Calibrating channel 4.") + channels['4'].append(_ir_calibrate(self._header[i], + data, 1, calibrate)) + if calibrate == 1: + self.units['4'] = 'K' + elif calibrate == 2: + self.units['4'] = 'W*m-2*sr-1*cm ?' + else: + self.units['4'] = '' + + if "5" in chns: + LOGGER.debug("Calibrating channel 5.") + channels['5'].append(_ir_calibrate(self._header[i], + data, 2, calibrate)) + if calibrate == 1: + self.units['5'] = 'K' + elif calibrate == 2: + self.units['5'] = 'W*m-2*sr-1*cm ?' + else: + self.units['5'] = '' + i += 1 + + # transfer channel data to class attributes + for ch_ in channels: + try: + self.channels[ch_] = np.vstack(channels[ch_]) + except ValueError: + self.channels[ch_] = None + if "3A" in chns or "3B" in chns: + self._is3b = np.vstack(is3b_all) + if "3A" in chns: + self.channels['3A'].mask = self._is3b * self.channels['3A'] if "3B" in chns: - ch3b = _ir_calibrate(self._header, self._data, 0, calibrate) - self.channels['3B'] = \ - np.ma.masked_array(ch3b, - np.logical_or((is3b is False) * \ - ch3b, - ch3b < 0.1)) - if calibrate == 1: - self.units['3B'] = 'K' - elif calibrate == 2: - self.units['3B'] = 'W*m-2*sr-1*cm ?' - else: - self.units['3B'] = '' - - if "4" in chns: - self.channels['4'] = _ir_calibrate(self._header, - self._data, 1, calibrate) - if calibrate == 1: - self.units['4'] = 'K' - elif calibrate == 2: - self.units['4'] = 'W*m-2*sr-1*cm ?' - else: - self.units['4'] = '' - - if "5" in chns: - self.channels['5'] = _ir_calibrate(self._header, - self._data, 2, calibrate) - if calibrate == 1: - self.units['5'] = 'K' - elif calibrate == 2: - self.units['5'] = 'W*m-2*sr-1*cm ?' - else: - self.units['5'] = '' + self.channels['3B'].mask = np.logical_or((self._is3b is False) * \ + self.channels['3B'], + self.channels['3B'] < 0.1) LOGGER.debug("Calibration time %s", str(datetime.datetime.now() - tic)) @@ -732,23 +776,28 @@ def show(data, negate=False): debug_on() SCENE = AAPP1b(sys.argv[1]) SCENE.read() - for name, val in zip(SCENE._header.dtype.names, SCENE._header[0]): + for name, val in zip(SCENE._header[0].dtype.names, SCENE._header[0][0]): print name, val - starttime = datetime.datetime(SCENE._header[0]["startdatayr"], 1, 1, 0, 0) - starttime += datetime.timedelta(days=int(SCENE._header[0]["startdatady"]) - 1, - seconds=SCENE._header[0]["startdatatime"] / 1000.0) + starttime = datetime.datetime(SCENE._header[0][0]["startdatayr"], + 1, 1, 0, 0) + starttime += \ + datetime.timedelta(days=int(SCENE._header[0][0]["startdatady"]) - 1, + seconds=SCENE._header[0][0]["startdatatime"] / \ + 1000.0) print "starttime:", starttime - endtime = datetime.datetime(SCENE._header[0]["enddatayr"], 1, 1, 0, 0) - endtime += datetime.timedelta(days=int(SCENE._header[0]["enddatady"]) - 1, - seconds=SCENE._header[0]["enddatatime"] / 1000.0) + endtime = datetime.datetime(SCENE._header[-1][0]["enddatayr"], 1, 1, 0, 0) + endtime += \ + datetime.timedelta(days=int(SCENE._header[-1][0]["enddatady"]) - 1, + seconds=SCENE._header[-1][0]["enddatatime"] / \ + 1000.0) print "endtime:", endtime # print SCENE._data['hrpt'].shape #show(SCENE._data['hrpt'][:, :, 4].astype(np.float)) # raw_input() SCENE.calibrate() SCENE.navigate() - for i_ in AVHRR_CHANNEL_NAMES: - data_ = SCENE.channels[i_] + for i__ in AVHRR_CHANNEL_NAMES: + data_ = SCENE.channels[i__] print >> sys.stderr, "%-3s" % i_, \ "%6.2f%%" % (100. * (float(np.ma.count(data_)) / data_.size)), \ "%6.2f, %6.2f, %6.2f" % (data_.min(), data_.mean(), data_.max())