#!/usr/bin/env python # -*- coding: utf-8 -*- # Copyright (c) 2012, 2013 SMHI # Author(s): # Martin Raspaud # Adam Dybbroe # Nina Håkansson # Oana Nicola # Lars Ørum Rasmussen # This program 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. # This program 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 this program. If not, see . """Reader for aapp level 1b data. http://research.metoffice.gov.uk/research/interproj/nwpsaf/aapp/ NWPSAF-MF-UD-003_Formats.pdf """ import numpy as np import os import logging import datetime import glob from ConfigParser import ConfigParser from mpop import CONFIG_PATH LOG = logging.getLogger('aapp1b') def load(satscene, *args, **kwargs): """Read data from file and load it into *satscene*. A possible *calibrate* keyword argument is passed to the AAPP reader. Should be 0 for off (counts), 1 for default (brightness temperatures and reflectances), and 2 for radiances only. """ del args conf = ConfigParser() conf.read(os.path.join(CONFIG_PATH, satscene.fullname + ".cfg")) options = {} for option, value in conf.items(satscene.instrument_name + "-level2", raw = True): options[option] = value options["calibrate"] = kwargs.get("calibrate", True) LOG.info("Loading instrument '%s'" % satscene.instrument_name) try: CASES[satscene.instrument_name](satscene, options) except KeyError: raise KeyError("Unknown instrument '%s'" % satscene.instrument_name) def load_avhrr(satscene, options): """Read avhrr data from file and load it into *satscene*. """ if "filename" not in options: raise IOError("No filename given, cannot load.") chns = satscene.channels_to_load & set(AVHRR_CHANNEL_NAMES) LOG.info("Loading channels " + str(sorted(list(chns)))) if len(chns) == 0: return values = {"orbit": satscene.orbit, "satname": satscene.satname, "number": satscene.number, "instrument": satscene.instrument_name, "satellite": satscene.fullname } filename = os.path.join(satscene.time_slot.strftime(options["dir"])%values, satscene.time_slot.strftime(options["filename"]) %values) file_list = glob.glob(filename) if len(file_list) > 1: raise IOError("More than one l1b file matching!") elif len(file_list) == 0: raise IOError("No l1b file matching!: " + filename) filename = file_list[0] LOG.debug("Loading from " + filename) scene = AAPP1b(filename) scene.read() scene.calibrate(chns, calibrate=options.get('calibrate', 1)) scene.navigate() for chn in chns: if scene.channels.has_key(chn): satscene[chn].data = scene.channels[chn] satscene[chn].info['units'] = scene.units[chn] try: from pyresample import geometry except ImportError, ex_: LOG.debug("Could not load pyresample: " + str(ex_)) satscene.lat = scene.lats satscene.lon = scene.lons else: satscene.area = geometry.SwathDefinition(lons=scene.lons, lats=scene.lats) AVHRR_CHANNEL_NAMES = ("1", "2", "3A", "3B", "4", "5") # AAPP 1b header _HEADERTYPE = np.dtype([("siteid", "S3"), ("blank", "S1"), ("l1bversnb", " np.expand_dims(data["calvis"][:, chn, 2, 4], 1) channel[mask1] = (channel * np.expand_dims(data["calvis"][:, chn, 2, 0] * 1e-10, 1) + np.expand_dims(data["calvis"][:, chn, 2, 1] * 1e-7, 1))[mask1] channel[mask2] = (channel * np.expand_dims(data["calvis"][:, chn, 2, 2] * 1e-10, 1) + np.expand_dims(data["calvis"][:, chn, 2, 3] * 1e-7, 1))[mask2] channel[channel < 0] = np.nan return np.ma.masked_array(channel, np.isnan(channel)) def _ir_calibrate(header, data, irchn, calib_type): """IR calibration *calib_type* = 0: Counts *calib_type* = 1: Reflectances *calib_type* = 2: Radiances """ count = data['hrpt'][:, :, irchn + 2].astype(np.float) if calib_type == 0: return count ir_const_1 = 1.1910659e-5 ir_const_2 = 1.438833 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) # 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 # Count to radiance conversion: rad = k1_ * count*count + k2_*count + k3_ if calib_type == 2: return rad 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(): LOG.info("Suspect scan lines: " + str(suspect_line_nums)) 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 tb_[tb_ <= 0] = np.nan # Data with count=0 are often related to erroneous (bad) lines, but in case # of saturation (channel 3b) count=0 can be observed and associated to a # real measurement. So we leave out this filtering to the user! # tb_[count == 0] = np.nan tb_[rad == 0] = np.nan return np.ma.masked_array(tb_, np.isnan(tb_)) def show(data, negate=False): """Show the stetched data. """ import Image as pil data = np.array((data - data.min()) * 255.0 / (data.max() - data.min()), np.uint8) if negate: data = 255 - data img = pil.fromarray(data) img.show() CASES = { "avhrr": load_avhrr, } if __name__ == "__main__": import sys from mpop.utils import debug_on debug_on() SCENE = AAPP1b(sys.argv[1]) SCENE.read() SCENE.calibrate() SCENE.navigate() 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()) show(SCENE.channels['4'], negate=True)