Coverage for runmacs/spec/calibration/calibrationdata.py : 64%

# -*- coding: utf-8 -*- 

""" 

calibrationdata 

=============== 

This module is responsible for storage and application of calibration data. 

""" 

import pint 

import numpy as np 

import sympy 

import inspect 

import datetime 

import dateutil.tz as datetz 

import dateutil.parser 

from runmacs.spec.io.envi import envi_to_dtype, dtype_to_envi, read_envi_header 

from runmacs.spec.calibration.badpixels import BadPixelFixer 

def getSymbols(expression): 

syms = set() 

for s in sympy.utilities.postorder_traversal(expression): 

if isinstance(s,sympy.Symbol): 

syms.add(s) 

return syms 

class CalibrationApplicator(object): 

def __init__(self, fields, rule): 

""" 

Can be called to apply calibration to raw data and integration time 

Can be sliced to create another instance which applies calibration to partial data 

""" 

self.fields = fields 

self.rule = rule 

def __getitem__(self, key): 

return CalibrationApplicator([f[key] for f in self.fields], self.rule) 

def __call__(self, rawData, tint): 

return self.rule(rawData, tint, *self.fields) 

class CalibrationData(object): 

""" 

Representation of specMACS calibration data. 

This class supports calibration data with arbitrary analytic nonlineary correction rules. 

""" 

def __init__(self): 

self.creationTime = datetime.datetime.utcnow().replace(tzinfo=datetz.tzutc()) 

self.dSymbol, self.tintSymbol, self.gainSymbol = map(sympy.Symbol, ['s0', 'tint', 'g1']) 

self.angleUnit = 'degree' 

self.wavelengthUnit = 'nm' 

self.calibratedUnit = 'mW m**-2 nm**-1 sr**-1' 

self.binningFactor = 1. 

self.unitFactor = 1. 

self.badPixels = [] 

@classmethod 

def fromFile(cls, filename): 

""" 

Create :py:class:`CalibrationData` from a given ``filename``. 

:param filename: Name of the calibration file (ENVI format, without file suffix) 

""" 

self = cls() 

header = read_envi_header(filename + '.hdr') 

self._fieldnames = header['fieldnames'] 

assert len(self._fieldnames) == int(header['lines']) 

self._sensorShape = (int(header['samples']), int(header['bands'])) 

self.creationTime = dateutil.parser.parse(header['creation time']) 

self.validFrom = dateutil.parser.parse(header['valid from']) 

self.validUntil = dateutil.parser.parse(header['valid until']) 

self.sensorId = header['sensorid'] 

self.symbolicNonlinearityRule = header['pixelwise nonlinearity rule'] 

self.calibratedUnit = header['calibrated unit'] 

self.angleUnit = header['angle unit'] 

self.wavelengthUnit = header['wavelength unit'] 

self.optics = header['optics'] 

try: 

hbp = header['bad pixels'] 

except KeyError: 

pass 

else: 

hbp = filter(lambda x: len(x) == 2, hbp) 

self.badPixels = [(int(a),int(b)) for a,b in hbp] 

self.data = np.fromfile(filename + '.cal', 

dtype=envi_to_dtype[header['data type']], 

count=(len(self._fieldnames) * self._sensorShape[0] * self._sensorShape[1]) 

).reshape(len(self._fieldnames), *reversed(self._sensorShape)).swapaxes(-2,-1) 

return self 

def toFile(self, filename): 

""" 

Saves the calibration data to a file with the given ``filename``. 

:param filename: Name of the calibration file (ENVI format, without file suffix) 

""" 

assert self.binningFactor == 1. #saving binned data is currently not supported 

assert self.unitFactor == 1. #unit factor other than 1 should be merged into g1 

with open(filename + '.hdr', 'w') as f: 

f.write(self.rawEnviHeader) 

with open(filename + '.cal', 'w') as f: 

f.write(self.data.swapaxes(-2,-1).tostring()) 

def genBinning(self, spatial=1, spectral=1): 

""" 

Create new calibration data for given binnings 

:param spatial: spatial binning factor 

:param spectral: spectral binning factor 

:returns: new :py:class:`CalibrationData` with given binning 

""" 

new = type(self)() 

new.setFieldNames(*self._fieldnames) 

new.setSensorShape(self._sensorShape[0] / spatial, self._sensorShape[1] / spectral) 

new.creationTime = self.creationTime 

new.validFrom = self.validFrom 

new.validUntil = self.validUntil 

new.optics = self.optics 

new.sensorId = self.sensorId 

new.symbolicNonlinearityRule = self.symbolicNonlinearityRule 

new.calibratedUnit = self.calibratedUnit 

new.angleUnit = self.angleUnit 

new.wavelengthUnit = self.wavelengthUnit 

new.data = self.data.reshape(self.data.shape[0], 

self._sensorShape[0]/spatial, spatial, 

self._sensorShape[1]/spectral, spectral).mean(axis=-1).mean(axis=-2) 

new.binningFactor = self.binningFactor * spatial * spectral 

new.unitFactor = self.unitFactor 

new.badPixels = [(int(a/spatial), int(b/spectral)) for a,b in self.badPixels] 

return new 

def convertUnit(self, targetUnit, ureg=None): 

""" 

Create new calibration data for another radiance unit. 

:param targetUnit: desired radiance unit 

:param ureg: :py:mod:`pint` unit registry, if ``None``, a temporary one is created 

:returns: new :py:class:`CalibrationData` with given radiance unit 

""" 

137 ↛ 138line 137 didn't jump to line 138, because the condition on line 137 was never true if ureg is None: 

ureg = pint.UnitRegistry() 

conversionFactor = ureg(self.calibratedUnit).to(ureg(targetUnit)).magnitude 

new = type(self)() 

new.setFieldNames(*self._fieldnames) 

new.setSensorShape(*self._sensorShape) 

new.creationTime = self.creationTime 

new.validFrom = self.validFrom 

new.validUntil = self.validUntil 

new.optics = self.optics 

new.sensorId = self.sensorId 

new.symbolicNonlinearityRule = self.symbolicNonlinearityRule 

new.calibratedUnit = targetUnit 

new.angleUnit = self.angleUnit 

new.wavelengthUnit = self.wavelengthUnit 

new.data = self.data 

new.binningFactor = self.binningFactor 

new.unitFactor = conversionFactor*self.unitFactor 

new.badPixels = self.badPixels 

return new 

@property 

def rawEnviHeader(self): 

""" 

The raw ENVI header which describes the current calibration data. 

""" 

return """ENVI 

description = { specMACS calibration data } 

file type = specMACS calibration data v2 

sensor type = specMACS 

acquitision date = DATE(dd-mm-yyy): %(date)s 

acquitision time = UTC TIME: %(time)s 

creation time = %(creationTime)s 

valid from = %(validFrom)s 

valid until = %(validUntil)s 

optics = %(optics)s 

sensorid = %(sensorId)s 

interleave = bil 

samples = %(spatialPixels)d 

lines = %(frames)d 

bands = %(spectralPixels)d 

data type = %(enviDatatype)s 

header offset = 0 

fieldnames = { %(enviFieldnames)s } 

pixelwise nonlinearity rule = %(nonlinearityrule)s 

pixelwise calibration rule = %(calibrationrule)s 

calibrated unit = %(calibratedUnit)s 

wavelength unit = %(wavelengthUnit)s 

angle unit = %(angleUnit)s 

bad pixels = { %(enviBadPixels)s } 

"""%self.metaDict 

@property 

def metaDict(self): 

assert self.data.shape[0] == len(self._fieldnames) 

return {'date': self.creationTime.strftime('%d-%m-%Y'), 

'time': self.creationTime.strftime('%H:%M:%S'), 

'creationTime': self.creationTime.isoformat(), 

'validFrom': self.validFrom.isoformat(), 

'validUntil': self.validUntil.isoformat(), 

'optics': self.optics, 

'sensorId': self.sensorId, 

'spatialPixels': self.shape[1], 

'spectralPixels': self.shape[2], 

'frames': self.shape[0], 

'enviDatatype': dtype_to_envi[self.data.dtype], 

'enviFieldnames': ', '.join(self._fieldnames), 

'nonlinearityrule': self.symbolicNonlinearityRule, 

'calibrationrule': self.symbolicCalibrationRule, 

'calibratedUnit': self.calibratedUnit, 

'wavelengthUnit': self.wavelengthUnit, 

'angleUnit': self.angleUnit, 

'enviBadPixels': ', '.join('{%d, %d}'%tuple(bp) for bp in self.badPixels) 

} 

@property 

def shape(self): 

return self.data.shape 

@property 

def symbolicNonlinearityRule(self): 

return self._symbolicNonlinearityRule 

@symbolicNonlinearityRule.setter 

def symbolicNonlinearityRule(self, rule): 

if isinstance(rule, (str, unicode)): 

rule = sympy.sympify(rule) 

syms = getSymbols(rule) 

self._nonlinearitySymbols = [self.dSymbol, self.tintSymbol] 

otherSymbols = [self.gainSymbol] + map(sympy.Symbol, self._fieldnames) 

self._nonlinearitySymbols += [s for s in otherSymbols if s in syms] 

for s in syms: 

assert s in self._nonlinearitySymbols 

self._symbolicNonlinearityRule = rule 

self._nonlinearityRule = sympy.lambdify(self._nonlinearitySymbols, rule) 

@property 

def symbolicCalibrationRule(self): 

return self.gainSymbol * self._symbolicNonlinearityRule 

@property 

def nonlinearityRule(self): 

return self._nonlinearityRule 

@nonlinearityRule.setter 

def nonlinearityRule(self, rule): 

argNames = inspect.getargspec(rule).args 

self._nonlinearitySymbols = map(sympy.Symbol, argNames) 

self._symbolicNonlinearityRule = rule(*self._nonlinearitySymbols) 

self._nonlinearityRule = rule 

@property 

def _calibrationSymbols(self): 

syms = [self.dSymbol, self.tintSymbol, self.gainSymbol] 

syms += [s for s in self._nonlinearitySymbols if s not in syms] 

return syms 

@property 

def calibrationRule(self): 

return sympy.lambdify(self._calibrationSymbols, self.symbolicCalibrationRule.subs(self.dSymbol, self.dSymbol*self.unitFactor/self.binningFactor)) 

@property 

def apply(self): 

""" 

Applies the calibration data to given raw data. 

This property can be used as a function, to calibrate from ``S0`` and ``tint`` to radiance units:: 

>>> caldata = CalibrationData.fromFile('cal...') 

>>> calibrated = caldata.apply(signal - dark_signal, tint) 

Alternatively, ``apply`` can be sliced, if only a part of the sensor needs to be calibrated:: 

>>> part_calibrated = caldata.apply[7,10:20](signal0[7,10:20], tint) 

""" 

return CalibrationApplicator([getattr(self, sym.name) for sym in self._calibrationSymbols[2:]], self.calibrationRule) 

def setValiditySpan(self, td): 

self.validFrom = self.creationTime - td 

self.validUntil = self.creationTime + td 

def setFieldNames(self, *fieldnames): 

self._fieldnames = list(fieldnames) 

self.resetData() 

def setSensorShape(self, spatial, spectral): 

self._sensorShape = (spatial, spectral) 

self.resetData() 

def resetData(self): 

try: 

self.data = np.zeros((len(self._fieldnames),) + self._sensorShape) 

except AttributeError: 

self.data = None 

def getBadPixelFixer(self, strategy=None): 

""" 

Get a :py:class:`BadPixelFixer` from the bad pixels given in the calibration data. 

:param strategy: Fixing strategy, gets passed to the :py:class:`BadPixelFixer`. 

""" 

if strategy is not None: 

return BadPixelFixer(self._sensorShape, self.badPixels, strategy) 

else: 

return BadPixelFixer(self._sensorShape, self.badPixels) 

def __getattr__(self, key): 

if key.startswith('_'): 

raise AttributeError 

try: 

fieldId = self._fieldnames.index(key) 

except ValueError: 

raise AttributeError 

return self.data[fieldId] 

def __setattr__(self, key, value): 

if key.startswith('_'): 

super(CalibrationData, self).__setattr__(key, value) 

try: 

fieldId = self._fieldnames.index(key) 

except (AttributeError, ValueError): 

super(CalibrationData, self).__setattr__(key, value) 

else: 

self.data[fieldId] = value