Source code for satpy.readers.aapp_l1b

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) 2012-2019 Satpy developers
#
# This file is part of satpy.
#
# satpy is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# satpy is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along with
# satpy.  If not, see <http://www.gnu.org/licenses/>.
"""Reader for aapp level 1b data.

Options for loading:

 - pre_launch_coeffs (False): use pre-launch coefficients if True, operational
   otherwise (if available).

http://research.metoffice.gov.uk/research/interproj/nwpsaf/aapp/
NWPSAF-MF-UD-003_Formats.pdf
"""

import logging
from datetime import datetime, timedelta

import numpy as np
import xarray as xr

import dask.array as da
from satpy.readers.file_handlers import BaseFileHandler
from satpy import CHUNK_SIZE

logger = logging.getLogger(__name__)

CHANNEL_NAMES = ['1', '2', '3a', '3b', '4', '5']

ANGLES = {'sensor_zenith_angle': 'satz',
          'solar_zenith_angle': 'sunz',
          'sun_sensor_azimuth_difference_angle': 'azidiff'}

PLATFORM_NAMES = {4: 'NOAA-15',
                  2: 'NOAA-16',
                  6: 'NOAA-17',
                  7: 'NOAA-18',
                  8: 'NOAA-19',
                  11: 'Metop-B',
                  12: 'Metop-A',
                  13: 'Metop-C',
                  14: 'Metop simulator'}


[docs]def create_xarray(arr): """Create xarray DataArray from numpy array.""" res = da.from_array(arr, chunks=(CHUNK_SIZE, CHUNK_SIZE)) res = xr.DataArray(res, dims=['y', 'x']) return res
[docs]class AVHRRAAPPL1BFile(BaseFileHandler): """Reader for AVHRR L1B files created from the AAPP software.""" def __init__(self, filename, filename_info, filetype_info): """Initialize object information by reading the input file.""" super(AVHRRAAPPL1BFile, self).__init__(filename, filename_info, filetype_info) self.channels = {i: None for i in AVHRR_CHANNEL_NAMES} self.units = {i: 'counts' for i in AVHRR_CHANNEL_NAMES} self._data = None self._header = None self._is3b = None self._is3a = None self._shape = None self.lons = None self.lats = None self.area = None self.sensor = 'avhrr-3' self.read() self.platform_name = PLATFORM_NAMES.get(self._header['satid'][0], None) if self.platform_name is None: raise ValueError("Unsupported platform ID: %d" % self.header['satid']) self.sunz, self.satz, self.azidiff = None, None, None @property def start_time(self): """Get the time of the first observation.""" return datetime(self._data['scnlinyr'][0], 1, 1) + timedelta( days=int(self._data['scnlindy'][0]) - 1, milliseconds=int(self._data['scnlintime'][0])) @property def end_time(self): """Get the time of the final observation.""" return datetime(self._data['scnlinyr'][-1], 1, 1) + timedelta( days=int(self._data['scnlindy'][-1]) - 1, milliseconds=int(self._data['scnlintime'][-1]))
[docs] def get_dataset(self, key, info): """Get a dataset from the file.""" if key.name in CHANNEL_NAMES: dataset = self.calibrate(key) elif key.name in ['longitude', 'latitude']: if self.lons is None or self.lats is None: self.navigate() if key.name == 'longitude': dataset = create_xarray(self.lons) else: dataset = create_xarray(self.lats) dataset.attrs = info else: # Get sun-sat angles if key.name in ANGLES: if isinstance(getattr(self, ANGLES[key.name]), np.ndarray): dataset = create_xarray(getattr(self, ANGLES[key.name])) else: dataset = self.get_angles(key.name) else: raise ValueError("Not a supported sun-sensor viewing angle: %s", key.name) dataset.attrs.update({'platform_name': self.platform_name, 'sensor': self.sensor}) dataset.attrs.update(key.to_dict()) for meta_key in ('standard_name', 'units'): if meta_key in info: dataset.attrs.setdefault(meta_key, info[meta_key]) if not self._shape: self._shape = dataset.shape return dataset
[docs] def read(self): """Read the data.""" tic = datetime.now() with open(self.filename, "rb") as fp_: header = np.memmap(fp_, dtype=_HEADERTYPE, mode="r", shape=(1, )) data = np.memmap(fp_, dtype=_SCANTYPE, offset=22016, mode="r") logger.debug("Reading time %s", str(datetime.now() - tic)) self._header = header self._data = data
[docs] def get_angles(self, angle_id): """Get sun-satellite viewing angles.""" sunz40km = self._data["ang"][:, :, 0] * 1e-2 satz40km = self._data["ang"][:, :, 1] * 1e-2 azidiff40km = self._data["ang"][:, :, 2] * 1e-2 try: from geotiepoints.interpolator import Interpolator except ImportError: logger.warning("Could not interpolate sun-sat angles, " "python-geotiepoints missing.") self.sunz, self.satz, self.azidiff = sunz40km, satz40km, azidiff40km else: cols40km = np.arange(24, 2048, 40) cols1km = np.arange(2048) lines = sunz40km.shape[0] rows40km = np.arange(lines) rows1km = np.arange(lines) along_track_order = 1 cross_track_order = 3 satint = Interpolator( [sunz40km, satz40km, azidiff40km], (rows40km, cols40km), (rows1km, cols1km), along_track_order, cross_track_order) self.sunz, self.satz, self.azidiff = satint.interpolate() return create_xarray(getattr(self, ANGLES[angle_id]))
[docs] def navigate(self): """Get the longitudes and latitudes of the scene.""" lons40km = self._data["pos"][:, :, 1] * 1e-4 lats40km = self._data["pos"][:, :, 0] * 1e-4 try: from geotiepoints import SatelliteInterpolator except ImportError: logger.warning("Could not interpolate lon/lats, " "python-geotiepoints missing.") self.lons, self.lats = lons40km, lats40km else: cols40km = np.arange(24, 2048, 40) cols1km = np.arange(2048) lines = lons40km.shape[0] rows40km = np.arange(lines) rows1km = np.arange(lines) along_track_order = 1 cross_track_order = 3 satint = SatelliteInterpolator( (lons40km, lats40km), (rows40km, cols40km), (rows1km, cols1km), along_track_order, cross_track_order) self.lons, self.lats = satint.interpolate()
[docs] def calibrate(self, dataset_id, pre_launch_coeffs=False, calib_coeffs=None): """Calibrate the data.""" if calib_coeffs is None: calib_coeffs = {} units = {'reflectance': '%', 'brightness_temperature': 'K', 'counts': '', 'radiance': 'W*m-2*sr-1*cm ?'} if dataset_id.name in ("3a", "3b") and self._is3b is None: # Is it 3a or 3b: is3b = np.expand_dims(np.bitwise_and(self._data['scnlinbit'], 3) == 1, 1) self._is3b = np.repeat(is3b, self._data['hrpt'][0].shape[0], axis=1) is3a = np.expand_dims(np.bitwise_and(self._data['scnlinbit'], 3) == 0, 1) self._is3a = np.repeat(is3a, self._data['hrpt'][0].shape[0], axis=1) try: vis_idx = ['1', '2', '3a'].index(dataset_id.name) ir_idx = None except ValueError: vis_idx = None ir_idx = ['3b', '4', '5'].index(dataset_id.name) if vis_idx is not None: coeffs = calib_coeffs.get('ch' + dataset_id.name) ds = create_xarray( _vis_calibrate(self._data, vis_idx, dataset_id.calibration, pre_launch_coeffs, coeffs, mask=(dataset_id.name == '3a' and np.logical_not(self._is3a)))) else: ds = create_xarray( _ir_calibrate(self._header, self._data, ir_idx, dataset_id.calibration, mask=(dataset_id.name == '3b' and np.logical_not(self._is3b)))) if dataset_id.name == '3a' and np.all(np.isnan(ds)): raise ValueError("Empty dataset for channel 3A") if dataset_id.name == '3b' and np.all(np.isnan(ds)): raise ValueError("Empty dataset for channel 3B") ds.attrs['units'] = units[dataset_id.calibration] ds.attrs.update(dataset_id._asdict()) return ds
AVHRR_CHANNEL_NAMES = ("1", "2", "3a", "3b", "4", "5") # AAPP 1b header _HEADERTYPE = np.dtype([("siteid", "S3"), ("blank", "S1"), ("l1bversnb", "<i2"), ("l1bversyr", "<i2"), ("l1bversdy", "<i2"), ("reclg", "<i2"), ("blksz", "<i2"), ("hdrcnt", "<i2"), ("filler0", "S6"), ("dataname", "S42"), ("prblkid", "S8"), ("satid", "<i2"), ("instid", "<i2"), ("datatype", "<i2"), ("tipsrc", "<i2"), ("startdatajd", "<i4"), ("startdatayr", "<i2"), ("startdatady", "<i2"), ("startdatatime", "<i4"), ("enddatajd", "<i4"), ("enddatayr", "<i2"), ("enddatady", "<i2"), ("enddatatime", "<i4"), ("cpidsyr", "<i2"), ("cpidsdy", "<i2"), ("filler1", "S8"), # data set quality indicators ("inststat1", "<i4"), ("filler2", "S2"), ("statchrecnb", "<i2"), ("inststat2", "<i4"), ("scnlin", "<i2"), ("callocscnlin", "<i2"), ("misscnlin", "<i2"), ("datagaps", "<i2"), ("okdatafr", "<i2"), ("pacsparityerr", "<i2"), ("auxsyncerrsum", "<i2"), ("timeseqerr", "<i2"), ("timeseqerrcode", "<i2"), ("socclockupind", "<i2"), ("locerrind", "<i2"), ("locerrcode", "<i2"), ("pacsstatfield", "<i2"), ("pacsdatasrc", "<i2"), ("filler3", "S4"), ("spare1", "S8"), ("spare2", "S8"), ("filler4", "S10"), # Calibration ("racalind", "<i2"), ("solarcalyr", "<i2"), ("solarcaldy", "<i2"), ("pcalalgind", "<i2"), ("pcalalgopt", "<i2"), ("scalalgind", "<i2"), ("scalalgopt", "<i2"), ("irttcoef", "<i2", (4, 6)), ("filler5", "<i4", (2, )), # radiance to temperature conversion ("albcnv", "<i4", (2, 3)), ("radtempcnv", "<i4", (3, 3)), ("filler6", "<i4", (3, )), # Navigation ("modelid", "S8"), ("nadloctol", "<i2"), ("locbit", "<i2"), ("filler7", "S2"), ("rollerr", "<i2"), ("pitcherr", "<i2"), ("yawerr", "<i2"), ("epoyr", "<i2"), ("epody", "<i2"), ("epotime", "<i4"), ("smaxis", "<i4"), ("eccen", "<i4"), ("incli", "<i4"), ("argper", "<i4"), ("rascnod", "<i4"), ("manom", "<i4"), ("xpos", "<i4"), ("ypos", "<i4"), ("zpos", "<i4"), ("xvel", "<i4"), ("yvel", "<i4"), ("zvel", "<i4"), ("earthsun", "<i4"), ("filler8", "S16"), # analog telemetry conversion ("pchtemp", "<i2", (5, )), ("reserved1", "<i2"), ("pchtempext", "<i2", (5, )), ("reserved2", "<i2"), ("pchpow", "<i2", (5, )), ("reserved3", "<i2"), ("rdtemp", "<i2", (5, )), ("reserved4", "<i2"), ("bbtemp1", "<i2", (5, )), ("reserved5", "<i2"), ("bbtemp2", "<i2", (5, )), ("reserved6", "<i2"), ("bbtemp3", "<i2", (5, )), ("reserved7", "<i2"), ("bbtemp4", "<i2", (5, )), ("reserved8", "<i2"), ("eleccur", "<i2", (5, )), ("reserved9", "<i2"), ("motorcur", "<i2", (5, )), ("reserved10", "<i2"), ("earthpos", "<i2", (5, )), ("reserved11", "<i2"), ("electemp", "<i2", (5, )), ("reserved12", "<i2"), ("chtemp", "<i2", (5, )), ("reserved13", "<i2"), ("bptemp", "<i2", (5, )), ("reserved14", "<i2"), ("mhtemp", "<i2", (5, )), ("reserved15", "<i2"), ("adcontemp", "<i2", (5, )), ("reserved16", "<i2"), ("d4bvolt", "<i2", (5, )), ("reserved17", "<i2"), ("d5bvolt", "<i2", (5, )), ("reserved18", "<i2"), ("bbtempchn3b", "<i2", (5, )), ("reserved19", "<i2"), ("bbtempchn4", "<i2", (5, )), ("reserved20", "<i2"), ("bbtempchn5", "<i2", (5, )), ("reserved21", "<i2"), ("refvolt", "<i2", (5, )), ("reserved22", "<i2"), ]) # AAPP 1b scanline _SCANTYPE = np.dtype([("scnlin", "<i2"), ("scnlinyr", "<i2"), ("scnlindy", "<i2"), ("clockdrift", "<i2"), ("scnlintime", "<i4"), ("scnlinbit", "<i2"), ("filler0", "S10"), ("qualind", "<i4"), ("scnlinqual", "<i4"), ("calqual", "<i2", (3, )), ("cbiterr", "<i2"), ("filler1", "S8"), # Calibration ("calvis", "<i4", (3, 3, 5)), ("calir", "<i4", (3, 2, 3)), ("filler2", "<i4", (3, )), # Navigation ("navstat", "<i4"), ("attangtime", "<i4"), ("rollang", "<i2"), ("pitchang", "<i2"), ("yawang", "<i2"), ("scalti", "<i2"), ("ang", "<i2", (51, 3)), ("filler3", "<i2", (3, )), ("pos", "<i4", (51, 2)), ("filler4", "<i4", (2, )), ("telem", "<i2", (103, )), ("filler5", "<i2"), ("hrpt", "<i2", (2048, 5)), ("filler6", "<i4", (2, )), # tip minor frame header ("tipmfhd", "<i2", (7, 5)), # cpu telemetry ("cputel", "S6", (2, 5)), ("filler7", "<i2", (67, )), ]) def _vis_calibrate(data, chn, calib_type, pre_launch_coeffs=False, calib_coeffs=None, mask=False): """Calibrate visible channel data. ``calib_type`` in count, reflectance, radiance. """ # Calibration count to albedo, the calibration is performed separately for # two value ranges. if calib_type not in ['counts', 'radiance', 'reflectance']: raise ValueError('Calibration ' + calib_type + ' unknown!') arr = data["hrpt"][:, :, chn] mask |= arr == 0 channel = arr.astype(np.float) if calib_type == 'counts': return channel if calib_type == 'radiance': logger.info("Radiances are not yet supported for " + "the VIS/NIR channels!") if pre_launch_coeffs: coeff_idx = 2 else: # check that coeffs are valid if np.all(data["calvis"][:, chn, 0, 4] == 0): logger.info( "No valid operational coefficients, fall back to pre-launch") coeff_idx = 2 else: coeff_idx = 0 intersection = data["calvis"][:, chn, coeff_idx, 4] if calib_coeffs is not None: logger.info("Updating from external calibration coefficients.") # intersection = np.expand_dims slope1 = np.expand_dims(calib_coeffs[0], 1) intercept1 = np.expand_dims(calib_coeffs[1], 1) slope2 = np.expand_dims(calib_coeffs[2], 1) intercept2 = np.expand_dims(calib_coeffs[3], 1) else: slope1 = np.expand_dims(data["calvis"][:, chn, coeff_idx, 0] * 1e-10, 1) intercept1 = np.expand_dims(data["calvis"][:, chn, coeff_idx, 1] * 1e-7, 1) slope2 = np.expand_dims(data["calvis"][:, chn, coeff_idx, 2] * 1e-10, 1) intercept2 = np.expand_dims(data["calvis"][:, chn, coeff_idx, 3] * 1e-7, 1) if chn == 2: slope2[slope2 < 0] += 0.4294967296 mask1 = channel <= np.expand_dims(intersection, 1) mask2 = channel > np.expand_dims(intersection, 1) channel[mask1] = (channel * slope1 + intercept1)[mask1] channel[mask2] = (channel * slope2 + intercept2)[mask2] channel = channel.clip(min=0) return np.where(mask, np.nan, channel) def _ir_calibrate(header, data, irchn, calib_type, mask=False): """Calibrate for IR bands. ``calib_type`` in brightness_temperature, radiance, count """ count = data["hrpt"][:, :, irchn + 2].astype(np.float) if calib_type == 0: return count # Mask unnaturally low values mask |= count == 0.0 k1_ = np.expand_dims(data['calir'][:, irchn, 0, 0] / 1.0e9, 1) k2_ = np.expand_dims(data['calir'][:, irchn, 0, 1] / 1.0e6, 1) k3_ = np.expand_dims(data['calir'][:, irchn, 0, 2] / 1.0e6, 1) # Count to radiance conversion: rad = k1_ * count * count + k2_ * count + k3_ all_zero = np.logical_and( np.logical_and( np.equal(k1_, 0), np.equal(k2_, 0)), np.equal(k3_, 0)) idx = np.indices((all_zero.shape[0], )) suspect_line_nums = np.repeat(idx[0], all_zero[:, 0]) if suspect_line_nums.any(): logger.info("Suspicious scan lines: %s", str(suspect_line_nums)) if calib_type == 2: mask |= rad <= 0.0 return np.where(mask, np.nan, rad) # Central wavenumber: cwnum = header['radtempcnv'][0, irchn, 0] if irchn == 0: cwnum = cwnum / 1.0e2 else: cwnum = cwnum / 1.0e3 bandcor_2 = header['radtempcnv'][0, irchn, 1] / 1e5 bandcor_3 = header['radtempcnv'][0, irchn, 2] / 1e6 ir_const_1 = 1.1910659e-5 ir_const_2 = 1.438833 t_planck = (ir_const_2 * cwnum) / \ np.log(1 + ir_const_1 * cwnum * cwnum * cwnum / rad) # Band corrections applied to t_planck to get correct # brightness temperature for channel: if bandcor_2 < 0: # Post AAPP-v4 tb_ = bandcor_2 + bandcor_3 * t_planck else: # AAPP 1 to 4 tb_ = (t_planck - bandcor_2) / bandcor_3 # Mask unnaturally low values # mask |= tb_ < 0.1 return np.where(mask, np.nan, tb_)