Coverage for runmacs/processor/utils.py : 50%

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

# Filename: utils.py 

""" 

utils for macsProcessor 

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

This is part of the macsProcessor suite. 

Copyright (C) 2014 Tobias Kölling 

""" 

from sh import ssh 

import os 

from itertools import izip 

from datetime import datetime, date 

import dateutil.tz as datetz 

import numpy as np 

import itertools 

import ephem 

def getFileViaSSH(server, filename, user=None, blockRange=None, blockSize=1): 

if user is not None: 

server = '%s@%s'%(user,server) 

if blockRange is None: 

print "getting file", filename, "totally" 

return ssh(server, 'cat', filename).stdout 

else: 

print "getting file", filename, "in blockrange", blockRange 

return ssh(server, 

'dd', 

'if="%s"'%filename, 

'bs=%d'%blockSize, 

'skip=%d'%blockRange[0], 

'count=%d'%(blockRange[1]-blockRange[0])).stdout 

def dateFromTimestamp(ts): 

return datetime.utcfromtimestamp(ts).replace(tzinfo=datetz.tzutc()) 

epoch = dateFromTimestamp(0) 

def date2seconds(d): 

return (d-epoch).total_seconds() 

def centerDate(a, b): 

return ((a-epoch) + (b-epoch)) / 2 + epoch 

def try_parse_date(potential_date): 

try: 

date = datetime.strptime(potential_date, '%Y-%m-%dT%H:%M:%SZ').replace(tzinfo=datetz.tzutc()) 

except (ValueError, TypeError): 

return potential_date 

return date 

def parse_dates_in_tree(tree): 

if isinstance(tree, dict): 

return {k:parse_dates_in_tree(v) for k,v in tree.iteritems()} 

elif isinstance(tree, list): 

return [parse_dates_in_tree(v) for v in tree] 

else: 

return try_parse_date(tree) 

def getAngleBreak(angles): 

""" 

Returns an angle which is the most distant from any angle in the angles set. 

If the biggest distance between angles in the angles list is across 0, it will return 0. 

""" 

s = list(sorted(angles)) 

zeroCrossingDistance = s[0] - s[-1] + 2*np.pi 

delta = [b-a for a,b in zip(s[:-1], s[1:])] 

biggestDistancePosition = np.argmax(delta) 

a = s[biggestDistancePosition] 

b = s[biggestDistancePosition+1] 

biggestDistance = b-a 

if zeroCrossingDistance >= biggestDistance: 

return 0. 

return (a + b) / 2. 

def centerAngleList(angles): 

""" 

Returns a new angles array, with the central value set to 0. 

If the list hast an even number of elements, set the average of the two central values to 0. 

""" 

l = len(angles) 

iCenter = int(l/2) 

if len(angles) % 2 == 1: 

center = angles[iCenter] 

else: 

center = (angles[iCenter - 1] + angles[iCenter]) / 2. 

return angles - center 

def constrainAzEl(az_, el_): 

el = ((el_ + np.pi/2.)%(2.*np.pi)) - (np.pi/2.) 

elMask = el>(np.pi/2.) 

el = np.where(elMask,np.pi-el,el) 

az = np.where(elMask,az_+np.pi,az_)%(2.*np.pi) 

return az, el 

def Rx(alpha): 

alpha = np.array(alpha) 

res = np.zeros(alpha.shape + (3,3)) 

res[...,0,0] = 1 

res[...,1,1] = np.cos(alpha) 

res[...,1,2] = -np.sin(alpha) 

res[...,2,1] = np.sin(alpha) 

res[...,2,2] = np.cos(alpha) 

return res 

def Ry(alpha): 

alpha = np.array(alpha) 

res = np.zeros(alpha.shape + (3,3)) 

res[...,1,1] = 1 

res[...,0,0] = np.cos(alpha) 

res[...,0,2] = np.sin(alpha) 

res[...,2,0] = -np.sin(alpha) 

res[...,2,2] = np.cos(alpha) 

return res 

def Rz(alpha): 

alpha = np.array(alpha) 

res = np.zeros(alpha.shape + (3,3)) 

res[...,2,2] = 1 

res[...,0,0] = np.cos(alpha) 

res[...,0,1] = -np.sin(alpha) 

res[...,1,0] = np.sin(alpha) 

res[...,1,1] = np.cos(alpha) 

return res 

def shapeBroadcastCompare(a,b): 

if len(a) != len(b): 

return False 

for ai,bi in zip(a,b): 

130 ↛ 131line 130 didn't jump to line 131, because the condition on line 130 was never true if ai != bi and ai != 1 and bi != 1: 

return False 

return True 

def diagDot(r1, r2): 

""" 

Dot product for last two axes, all prepended axes must have the same shape and the diagonal is returned. 

If the shape indicates that one of the arguments contains a vector and not a matrix, 

rules of vector-matrix multiplication are applied. 

This is what is needed for the dot product of "lists" of rotation matrices. 

""" 

if len(r1.shape) >= 2 and len(r2.shape) >= 2 and shapeBroadcastCompare(r1.shape[:-2], r2.shape[:-2]) and r1.shape[-1] == r2.shape[-2]: 

return np.einsum('...ij,...jk->...ik',r1,r2) 

144 ↛ 146line 144 didn't jump to line 146, because the condition on line 144 was never false elif len(r1.shape) >= 2 and len(r2.shape) >= 1 and shapeBroadcastCompare(r1.shape[:-2], r2.shape[:-1]) and r1.shape[-1] == r2.shape[-1]: 

return np.einsum('...ij,...j->...i',r1,r2) 

elif len(r1.shape) >= 1 and len(r2.shape) >= 2 and shapeBroadcastCompare(r1.shape[:-1], r2.shape[:-2]) and r1.shape[-1] == r2.shape[-2]: 

return np.einsum('...j,...jk->...k',r1,r2) 

raise ValueError('shape missmatch %s and %s'%(r1.shape, r2.shape)) 

def azel2ned(azel): 

az,el = azel 

return np.array((np.cos(az)*np.cos(el),np.sin(az)*np.cos(el),-np.sin(el))) 

def ned2azel(ned): 

n,e,d = ned 

s = d/((ned**2).sum(axis=0)) 

if not isinstance(s, np.ndarray): 

s = np.array(s) 

el = -np.arcsin(np.where(s<-1., -1., np.where(s>1., 1., s))) 

return np.array((np.arctan2(e,n)%(np.pi*2), el)) 

def rotAzEl(inp, delta, mirror=False): 

""" 

Makes a forward rotation in azel coordinates (e.g. sun-system to observer system) 

:param inp:nd-array with az,el as 0-th axis 

""" 

R = diagDot(Rz(delta[0]),Ry(delta[1])) 

if mirror: 

R = -R 

return ned2azel(diagDot(R,azel2ned(inp))) 

def rotAzElRev(inp, delta, mirror=False): 

""" 

Makes a reverse rotation in azel coordinates (e.g. observer system to sun system) 

:param inp:nd-array with az,el as 0-th axis 

""" 

R = diagDot(Ry(-delta[1]),Rz(-delta[0])) 

if mirror: 

R = -R 

return ned2azel(diagDot(R,azel2ned(inp))) 

def rpy2R(roll, pitch, yaw): 

""" 

Converts roll, pitch, yaw to a rotation matrix which transforms from RPY coordinates to NED coordinates. 

""" 

return reduce(diagDot, (Rx(roll),Ry(pitch),Rz(yaw))) 

186 ↛ exitline 186 didn't return from function 'pos', because the return on line 186 wasn't executeddef pos(body): return body.alt, body.az 

def sunPosition(time, lat, lon, alt): 

observer = ephem.Observer() 

observer.date = time 

observer.lat = str(lat) 

observer.lon = str(lon) 

observer.elevation = alt 

sun = ephem.Sun(observer) 

return sun 

def calcScatGeo(time, lat, lon, alt, rolls, pits, yaws, vzas): 

import numpy as np 

sun=[] 

for t,la,lo,al in izip(time, lat, lon, alt): 

sun.append(pos(sunPosition(t,la,lo,al))) 

sun = np.array(sun) 

sun_ned = azel2ned((sun[:,1],sun[:,0])) 

fov_uav = (np.array([[0.]*len(vzas),np.cos(vzas),-1.*np.sin(vzas)])[np.newaxis,...]) 

Rot = rpy2R(rolls,pits,yaws) 

fov_ned = diagDot(Rot, fov_uav) 

scat = np.einsum('...ij,i...->...j',fov_ned, sun_ned) 

fov_sca = np.arccos(scat) 

fov_phi = np.arctan2(fov_ned[:,1,:],fov_ned[:,0,:]) 

fov_vza = np.arccos(fov_ned[:,2,:]) 

sun_phi = np.ones(fov_vza.shape)*sun[:,1][:,np.newaxis] 

sun_sza = np.ones(fov_vza.shape)*sun[:,0][:,np.newaxis] 

return (fov_sca, fov_phi, fov_vza, sun_phi, sun_sza) 

223 ↛ exitline 223 didn't run the lambda on line 223dthandler = lambda obj: ( 

obj.isoformat() 

if isinstance(obj, datetime) 

or isinstance(obj, date) 

else None) 

def transformTime(t): 

return np.array(map(date2seconds, t)) 

class presentAttrsAsItems(object): 

def __init__(self, obj): 

self.obj = obj 

def __getitem__(self, name): 

return getattr(self.obj, name) 

def dictZip(**dct): 

keys, values = zip(*dct.items()) 

for v in itertools.izip(*values): 

yield dict(zip(keys, v)) 

defaultPartDescriptions = { 

'radianceData': {'long_name': 'Spectral radiance', 

'standard_name': 'spectral_radiance', 

'units': 'mW m**-2 nm**-1 sr**-1'}, 

'rawRadianceData': {'long_name': 'Uncalibrated Spectral radiance', 

'standard_name': 'uncalibrated_spectral_radiance'}, 

'rawRadianceError': {'long_name': 'Absolute Error of Uncalibrated Spectral radiance', 

'standard_name': 'uncalibrated_spectral_radiance_error'}, 

'previewdata': {'long_name': 'Preview data', 

'standard_name': 'preview_data'}, 

'az': {'long_name': 'Azimuth angle', 

'standard_name': 'azimuth', 

'units': 'radian'}, 

'el': {'long_name': 'Elevation angle', 

'standard_name': 'elevation', 

'units': 'radian'}, 

'positions': {'long_name': 'Angular pixel positions (az/el)', 

'standard_name': 'angular_pixel_positions', 

'units': 'radian'}, 

'sh': {'long_name': 'Sun relative azimuth angle (around the sun, 0=up)', 

'standard_name': 'sun_relative_azimuth', 

'units': 'radian'}, 

'sv': {'long_name': 'Sun relative elevation angle (90deg = center of sun)', 

'standard_name': 'sun_relative_elevation', 

'units': 'radian'}, 

'shv': {'long_name': 'Angular pixel positions sun relative (circular / radial)', 

'standard_name': 'angular_pixel_positions_sun_relative', 

'units': 'radian'}, 

'sca': {'long_name': 'Solar scattering angle', 

'standard_name': 'solar_scattering_angle', 

'units': 'radian'}, 

'wvlns': {'dimensions': ('spatial',), 

'long_name': 'Wavelength', 

'standard_name': 'wavelength', 

'units': 'nm'}, 

'wavelength': {'long_name': 'Wavelength', 

'standard_name': 'wavelength', 

'units': 'nm'}, 

'time': {'long_name': 'time', 

'standard_name': 'time', 

'units': 'Seconds since 1970-1-1 00:00:00', 

'_transform': transformTime}, 

'rawtime': {'long_name': 'time', 

'standard_name': 'time', 

'units': 'Seconds since 1970-1-1 00:00:00'}, 

'alpha': {'long_name': 'Deviation of vertical pixel angle from mount axis', 

'standard_name': 'camera_alpha', 

'units': 'radian'}, 

'act': {'long_name': 'across track angle', 

'standard_name': 'act', 

'units': 'radian'}, 

'alt': {'long_name': 'along track angle', 

'standard_name': 'alt', 

'units': 'radian'}, 

'forwardmap': {'long_name': 'Mapping from source pixels to target grid coordinates', 

'standard_name': 'projection_forward_map'}, 

'reversemap': {'long_name': 'Mapping from target pixels to source pixels', 

'standard_name': 'projection_reverse_map'}, 

'grid': {'long_name': 'Target image coordinate grid', 

'standard_name': 'projection_target_grid'}, 

'transmission': {'long_name': 'Transmission', 

'standard_name': 'transmission'}, 

'lat': {'long_name': 'latitude position above WGS84', 

'standard_name': 'latitude', 

'units': 'degree'}, 

'lon': {'long_name': 'longitude position above WGS84', 

'standard_name': 'longitude', 

'units': 'degree'}, 

'heightWGS84': {'long_name': 'height above WGS84', 

'standard_name': 'height', 

'units': 'm'}, 

'roll': {'long_name': 'roll angle', 

'standard_name': 'roll', 

'units': 'radian'}, 

'pitch': {'long_name': 'pitch angle', 

'standard_name': 'pitch', 

'units': 'radian'}, 

'yaw': {'long_name': 'yaw angle', 

'standard_name': 'yaw', 

'units': 'radian'}, 

}