Coverage for runmacs/spec/retrieval/cloudside/asap.py : 20%

# Filename: asap.py 

######################################################################### 

# 

# cloudside/asap.py - This file is part of the Munich Aerosol Cloud Scanner package. 

# 

# Copyright (C) 2012-2018 Florian Ewald 

# 

# runMACS 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 2 

# of the License, or (at your option) any later version. 

# 

# Spectral Python 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 software; if not, write to 

# 

# Free Software Foundation, Inc. 

# 59 Temple Place, Suite 330 

# Boston, MA 02111-1307 

# USA 

# 

######################################################################### 

# 

# Send comments to: 

# Florian Ewald, florian.ewald@campus.lmu.de 

# 

import __builtin__ 

try: 

__builtin__.profile 

except AttributeError: 

# No line profiler, provide a pass-through version 

def profile(func): return func 

__builtin__.profile = profile 

import numpy as np 

import numpy.ma as ma 

from scipy.ndimage import gaussian_filter 

from itertools import izip 

import scipy 

try: 

import cPickle as pickle 

except: 

import pickle 

@profile 

def parse_config(path): 

config={} 

for folder in path.split('/'): 

if 'phi_' in folder: config['phi1'] = int(folder.split('_')[-1]) 

if 'phi0' in folder: config['phi0'] = int(folder.split('_')[-1]) 

if 'vza_' in folder: config['vza'] = int(folder.split('_')[-1]) 

if 'sza_' in folder: config['sza'] = int(folder.split('_')[-1]) 

if 'ch_' in folder: config['cha'] = int(folder.split('_')[-1]) 

if folder.isdigit(): config['case'] = int(folder) 

return config 

@profile 

def calc_bsca(phi1, phi0, vza, sza, mystic=True): 

import numpy as np 

# convert from mystic definition to spherical coordinates 

if mystic: 

phi0 = (-1*(phi0-450.))%360 

phi1 = (-1*(phi1-270.))%360 

rphi0 = np.deg2rad(phi0) 

rphi1 = np.deg2rad(phi1) 

rsza = np.deg2rad(sza) 

rvza = np.deg2rad(vza) 

sun_vec = (np.sin(rsza)*np.cos(rphi0),np.sin(rsza)*np.sin(rphi0),np.cos(rsza)) 

fov_vec = (np.sin(rvza)*np.cos(rphi1),np.sin(rvza)*np.sin(rphi1),np.cos(rvza)) 

bsca = np.rad2deg(np.arccos(np.einsum('i...,i...->...',fov_vec, sun_vec))) 

return bsca 

def dot_product(x, y): 

return np.einsum('...i,...i->...',x, y) 

def norm(x): 

return np.sqrt((x*x).sum(axis=2)) 

def normalize(x): 

return (x.T/norm(x).T).T 

def project_onto_plane(x, n): 

d = dot_product(x, n)/norm(n) 

p = (d.T * normalize(n).T).T 

return (x.T - p.T).T 

@profile 

def calc_cloudsurf(dx, dy, dz, phi1, phi0, vza, sza, panorama=None): 

import numpy as np 

from scipy.ndimage.filters import gaussian_filter 

nsize = len(dx.flat) 

nx, ny = dx.shape 

if panorama: 

vza = np.rot90(np.linspace(vza+panorama[2],vza+panorama[3],ny).repeat(nx).reshape(ny,nx),1) 

phi1 = np.linspace(phi1+panorama[0],phi1+panorama[1],nx).repeat(ny).reshape(nx,ny) 

else: 

vza = np.rot90(np.linspace(vza,vza,ny).repeat(nx).reshape(ny,nx),1) 

phi1 = np.linspace(phi1,phi1,nx).repeat(ny).reshape(nx,ny) 

phi0 = phi1*0. + phi0 

sza = vza*0. + sza 

# convert from mystic definition to spherical coordinates 

phi0sp = (-1*(phi0-450.))%360 

phi1sp = (-1*(phi1-270.))%360 

rphi0 = np.deg2rad(phi0sp) 

rphi1 = np.deg2rad(phi1sp) 

rsza = np.deg2rad(180.-sza) 

rvza = np.deg2rad(180.-vza) 

sun_vec = np.dstack((np.sin(rsza)*np.cos(rphi0),np.sin(rsza)*np.sin(rphi0),np.cos(rsza))) 

fov_vec = np.dstack((np.sin(rvza)*np.cos(rphi1),np.sin(rvza)*np.sin(rphi1),np.cos(rvza))) 

cld_vec = np.dstack((gaussian_filter(dx,2), gaussian_filter(dy,2), gaussian_filter(dz,2))) 

pmag = np.sqrt((cld_vec*cld_vec).sum(axis=2)) 

cld_vec = (cld_vec.T/pmag.T).T 

#print(fov_vec[250,150],sun_vec[250,150],cld_vec[250,150]) 

rbsc = np.einsum('...i,...i->...',fov_vec, sun_vec) 

rsza = np.einsum('...i,...i->...',cld_vec, sun_vec) 

rvza = np.einsum('...i,...i->...',cld_vec, fov_vec) 

cfov = normalize(project_onto_plane(fov_vec, cld_vec)) 

csun = normalize(project_onto_plane(sun_vec, cld_vec)) 

rphi = np.einsum('...i,...i->...',cfov, csun) 

#rvza = rvza * (rvza >= 0.0) - 1.0 * rvza * (np.bitwise_and(rvza < 0.0, rvza >= -1.0)) - 999.9 * (rvza < -1.0) 

#rvza = gaussian_filter(rvza,1) 

#rsza = gaussian_filter(rsza,1) 

#rphi = gaussian_filter(rphi,1) 

#rphi = np.pi - rphi #MYSTIC Definition 

return (rvza, rsza, rphi, rbsc, pmag, cld_vec, phi1, phi0, vza, sza) 

class LUTInterpolator(object): 

""" 

Linear interpolation in multidimensional Look-Up Table 

:param bins: Sequence of number of bins for each dimension of the LUT. 

""" 

@profile 

def __init__(self, bins): 

self.nbins = len(bins) 

self.bins = bins 

def setup(self, **kwargs): 

""" 

Setup of interpolation options. 

The arguments have to be a sequene of values. 

Each sequence must contain one value for each LUT dimension. 

:param outliers: If True: look-ups outside of the table will be 

extrapolated. Values will be set to NAN otherwise. 

Default is False. 

:param interpol: If True: interpolate, if False: use nearest neighbor. 

Default is True. 

""" 

self.outliers = kwargs.pop('outliers', [False] * self.nbins) 

self.interpol = kwargs.pop('interpol', [True] * self.nbins) 

@profile 

def locate_in_lut(self, fields): 

""" 

Find location in LUT according to the interpolation options. 

In the following, N is the number of dimensions in the LUT. 

:param fields: N-element sequence of positions along LUT axes 

to search for. 

Each element of the sequence can either be a number 

or an arbitrary dimensional numpy array. If they are 

numpy arrays, the results will have the axes of these 

arrays as additional axes 

(according to numpy broadcasting rules). 

:returns: 2-tuple containing `indexes` and `weights`. 

* The 0th axis of `indexes` will specify the LUT-dimension. 

* The axes 1...N axes of `indexes` specify multiple samples 

to take from the LUT in order to interpolate later. 

* The axes 0...N-1 of `weights` contain the interpolation 

weights for the samples taken with `indexes`. 

* The following axes are the sames as from the input. 

This method can be used as follows (assuming 2D LUT): 

.. code-block:: python 

indexes, weights = interpolator.locate_in_lut([1.7, 4.2]) 

result = lut.get_from_lut(*indexes) 

result_interpolated = np.einsum('ij...,ij...->...', result, weights) 

``result_interpolated`` will now contain the interpolated value at the 

location :math:`(1.7, 4.2)` of the Look-Up Table. 

The interpolator can also be used with less dimensions than the dimension 

of the LUT. Naturally, the result of the lookup will be an array containing 

the remaining dimensions. The user will be responsible to add missing axes 

for the ``einsum`` above as needed. 

""" 

indexes = [] 

weights = [] 

for i, (field, ibin, outlier, inter) in enumerate(izip(fields, self.bins, self.outliers, self.interpol)): 

left, right, w = self.digitize(field, np.min(ibin), np.max(ibin), len(ibin), outlier=outlier) 

indexes.append(self.idx(left, right, i, self.nbins)) 

weights.append(self.w(w, i, self.nbins, interpol=inter)) 

weights = np.prod(weights) 

return indexes, weights 

@profile 

def digitize(self, vals, bin_min, bin_max, nbins, outlier=True): 

wout = 0 

if outlier is False: wout = np.nan 

bin_step = (bin_max-bin_min)/np.float(nbins-1) 

fb = (vals-(bin_min+bin_step/2.))/(bin_step) 

bb = fb.astype('int') 

ba = bb+1 

fr = 1.-(fb-bb) 

mask_low = fb<0 

mask_lbo = ba==0 

mask_hig = fb>=(nbins-2) 

mask_hbo = ba==(nbins-2) 

fw = np.array((fr, 1.-fr)) 

bb[mask_low] = 0 

ba[mask_low] = 0 

bb[mask_hig] = nbins-2 

ba[mask_hig] = nbins-2 

270 ↛ 271line 270 didn't jump to line 271, because the condition on line 270 was never true if outlier is True: 

fw[0,mask_low] = 0 

fw[0,mask_hig] = 0 

fw[1,mask_low] = 1 

fw[1,mask_hig] = 1 

else: 

fw[:,mask_low] = 0 

fw[:,mask_hig] = 0 

fw[1,mask_lbo] = 1 

fw[1,mask_hbo] = 1 

return bb, ba, fw 

@profile 

def s(self, i,n): 

a = [1]*n 

a[i] = 2 

return tuple(a) 

@profile 

def w(self, f, i, n, interpol=True): 

if interpol is False: f = np.round(f) 

#f = np.array((f, 1.-f)) 

return f.reshape(self.s(i,n)+f.shape[1:]) 

@profile 

def idx(self, b,a,i,n): 

ba = np.array((b,a)) 

return ba.reshape(self.s(i,n)+ba.shape[1:]) 

class MYSTICrun(object): 

""" 

Class describing MYSTIC output with auxiliary information 

""" 

def __init__(self): 

self.rad = None 

self.cld = None 

self.conf = None 

self.reff = None 

self.tau = None 

self.geo = None 

@classmethod 

def fromPath(cls, radpath, cldpath=None, pattern='reptran.cldprp', panorama=[-45,+45,-23,+23]): 

from runmacs.rtm.uvspec import Output 

import os 

new = cls() 

out = Output() 

if panorama is not False: 

cls.panorama = panorama 

radpath = os.path.abspath(radpath) 

new.rad = out.read_mystic_image(radpath) 

new.conf = parse_config(radpath) 

if cldpath: 

cldpath = os.path.abspath(cldpath) 

new.cld = out.read_mystic_image(cldpath, pattern=pattern, output="mc.mean.cldprp") 

new.calc_aux() 

return new 

def calc_aux(self): 

""" 

Read MYSTIC Output 

""" 

from scipy.ndimage.morphology import binary_erosion 

reffwc = self.cld[:,:,0] 

reffic = self.cld[:,:,1] 

hrefwc = self.cld[:,:,2] 

hrefic = self.cld[:,:,3] 

tauswc = self.cld[:,:,4] 

tausic = self.cld[:,:,5] 

dxwc = self.cld[:,:,6] 

dywc = self.cld[:,:,7] 

dzwc = self.cld[:,:,8] 

dxic = self.cld[:,:,9] 

dyic = self.cld[:,:,10] 

dzic = self.cld[:,:,11] 

wc = tauswc>1.0 

ic = tausic>1.0 

self.phase = tauswc*0. 

self.phase[wc>0] = 1 

self.phase[ic>0] = 2 

self.dx = dxwc*(self.phase==1) + dxic*(self.phase==2) 

self.dy = dywc*(self.phase==1) + dyic*(self.phase==2) 

self.dz = dzwc*(self.phase==1) + dzic*(self.phase==2) 

self.reff = reffwc*(self.phase==1) + reffic*(self.phase==2) 

self.href = hrefwc*(self.phase==1) + hrefic*(self.phase==2) 

self.taus = tauswc+ tausic 

self.geo = calc_cloudsurf(self.dx, self.dy, self.dz,self.conf['phi1'], 

self.conf['phi0'],self.conf['vza'],self.conf['sza'],panorama=self.panorama) 

self.cumu = self.geo[0] 

self.csza = np.rad2deg(np.arccos(self.geo[1])) 

self.cphi = 180.-np.rad2deg(np.arccos(self.geo[2])) 

self.bsca = np.rad2deg(np.arccos(self.geo[3])) 

self.rad = np.rot90(self.rad) 

self.reff = np.rot90(self.reff) 

self.href = np.rot90(self.href) 

self.taus = np.rot90(self.taus) 

self.phase = np.rot90(self.phase) 

self.cumu = np.rot90(self.cumu) 

self.csza = np.rot90(self.csza) 

self.cphi = np.rot90(self.cphi) 

self.bsca = np.rot90(self.bsca) 

self.phi1 = np.rot90(self.geo[6]) 

self.phi0 = np.rot90(self.geo[7]) 

self.vza = np.rot90(self.geo[8]) 

self.sza = np.rot90(self.geo[9]) 

class ASAPModel(object): 

def __init__(self, **kwargs): 

import os 

from netCDF4 import Dataset 

self.rundir = os.path.abspath(os.path.dirname('__file__')) 

self.savstr = kwargs.pop('savstr', 'RCE') 

self.savdir = kwargs.pop('savdir', self.rundir+'/sav') 

channel = kwargs.pop('channel', None) 

self.I0 = kwargs.pop('I0', None) 

lut1d_wc = kwargs.pop('lut1d_wc', None) 

lut1d_ic = kwargs.pop('lut1d_ic', None) 

if 'refl008' in channel: 

self.I0 = 9.672943e+02 

lut1d_wc = kwargs.pop('lut1d_wc', self.rundir+'/luts/'+'ch008_wc_luts_mie_afglms_cloudside.dat.cdf') 

lut1d_ic = kwargs.pop('lut1d_ic', self.rundir+'/luts/'+'ch008_ic_luts_baum_afglms_cloudside.dat.cdf') 

if 'refl021' in channel: 

self.I0 = 9.793503e+01 

lut1d_wc = kwargs.pop('lut1d_wc', self.rundir+'/luts/'+'ch021_wc_luts_mie_afglms_cloudside.dat.cdf') 

lut1d_ic = kwargs.pop('lut1d_ic', self.rundir+'/luts/'+'ch021_ic_luts_baum_afglms_cloudside.dat.cdf') 

self.channel = channel 

self.lut_wc = Dataset(lut1d_wc) 

self.lut_ic = Dataset(lut1d_ic) 

self.lut_wc_refl = np.array(self.lut_wc.variables[channel]) 

self.lut_ic_refl = np.array(self.lut_ic.variables[channel]) 

self.bins_sza = self.lut_wc.variables['sza'][:] 

self.bins_umu = self.lut_wc.variables['umu'][:] 

self.bins_phi = self.lut_wc.variables['phi'][:] 

self.bins_tau = self.lut_wc.variables['tau'][:] 

self.bins_ref = self.lut_wc.variables['reff'][:] 

self.radfile1 = None 

self.radfile2 = None 

def load_sza(self, rad1_dir, cld1_dir, rad2_dir, cld2_dir, panorama=True): 

from runmacs.rtm.uvspec import Output 

from scipy.ndimage.morphology import binary_erosion 

out = Output() 

self.rad1 = out.read_mystic_image(rad1_dir) 

self.rad2 = out.read_mystic_image(rad2_dir) 

self.cld1 = out.read_mystic_image(cld1_dir, pattern='reptran.cldprp',output="mc.mean.cldprp") 

self.cld2 = out.read_mystic_image(cld2_dir, pattern='reptran.cldprp',output="mc.mean.cldprp") 

self.conf1 = parse_config(cld1_dir) 

self.conf2 = parse_config(cld2_dir) 

if panorama: 

self.panorama1 = [-45,+45,-23,+23] 

self.panorama2 = [-45,+45,-23,+23] 

else: 

self.panorama1 = None 

self.panorama2 = None 

self.reff1 = self.cld1[:,:,0] 

self.reff2 = self.cld2[:,:,0] 

self.taus1 = self.cld1[:,:,4] 

self.taus2 = self.cld2[:,:,4] 

self.dx1 = self.cld1[:,:,6] 

self.dy1 = self.cld1[:,:,7] 

self.dz1 = self.cld1[:,:,8] 

self.dx2 = self.cld2[:,:,6] 

self.dy2 = self.cld2[:,:,7] 

self.dz2 = self.cld2[:,:,8] 

self.geo1 = calc_cloudsurf(self.dx1, self.dy1, self.dz1,self.conf1['phi1'], 

self.conf1['phi0'],self.conf1['vza'],self.conf1['sza'],panorama=self.panorama1) 

self.cumu1 = self.geo1[0] 

self.csza1 = np.rad2deg(np.arccos(self.geo1[1])) 

self.cphi1 = np.rad2deg(np.arccos(self.geo1[2])) 

self.bsc1 = np.rad2deg(np.arccos(self.geo1[3])) 

self.geo2 = calc_cloudsurf(self.dx2, self.dy2, self.dz2,self.conf2['phi1'], 

self.conf2['phi0'],self.conf2['vza'],self.conf2['sza'],panorama=self.panorama2) 

self.cumu2 = self.geo2[0] 

self.csza2 = np.rad2deg(np.arccos(self.geo2[1])) 

self.cphi2 = np.rad2deg(np.arccos(self.geo2[2])) 

self.bsc2 = np.rad2deg(np.arccos(self.geo2[3])) 

self.rads1 = self.get_from_lut(self.csza1, self.cumu1, self.cphi1, self.reff1, self.taus1) 

self.rads2 = self.get_from_lut(self.csza2, self.cumu2, self.cphi2, self.reff2, self.taus2) 

#self.rad1_match = np.argmin(np.abs(self.rads1.T-self.rad1.T).T,axis=2) 

#self.rad2_match = np.argmin(np.abs(self.rads2.T-self.rad2.T).T,axis=2) 

#nx, ny = self.rad1.shape 

#a=np.arange(nx).reshape(nx,1) 

#b=np.arange(ny).reshape(1,ny) 

#self.rad1_match_idx = (a,b,self.rad1_match) 

#self.rad2_match_idx = (a,b,self.rad2_match) 

kernel=np.array([[0.,0.5,0.5,0],[0.5,1.,1.,0.5],[0.5,1.,1.,0.5],[0.,0.5,0.5,0.]]) 

self.no_shadow = np.bitwise_and(self.csza2<90.,self.csza1<90).astype('float') 

self.no_shadow = gaussian_filter(binary_erosion(self.no_shadow,np.ones((6,6))).astype('float'), 4) 

def init_szas(self, path): 

import os 

import fnmatch 

radfiles = np.unique([dirpath 

for dirpath, dirnames, files in os.walk(path) 

for f in fnmatch.filter(files, '*rad.spc')]) 

geofiles = np.unique([dirpath 

for dirpath, dirnames, files in os.walk(path) 

for f in fnmatch.filter(files, '*cldprp')]) 

if 'refl008' in self.channel: 

chan='ch_870' 

if 'refl021' in self.channel: 

chan='ch_2100' 

self.radfiles = [r for r in radfiles if chan in r] 

self.geofiles = [g for g in geofiles if 'ch_2100' in g] 

self.rad_configs = [parse_config(r) for r in self.radfiles] 

self.geo_configs = [parse_config(g) for g in self.geofiles] 

def get_sza(self, sza): 

import os 

import numpy as np 

import shelve 

import bisect 

sza = np.float(sza) 

szas = np.array([conf['sza'] for conf in self.rad_configs], dtype='float') 

if sza in szas: 

idx_sza1 = np.where(szas==sza)[0][0] 

idx_sza2 = np.where(szas==sza)[0][0] 

else: 

idx_sza2 = bisect.bisect(szas, sza) 

idx_sza1 = idx_sza2-1 

if idx_sza2>=len(szas): 

idx_sza2=len(szas)-1 

if idx_sza1<0: 

idx_sza1=0 

radfile1 = self.radfiles[idx_sza1] 

radfile2 = self.radfiles[idx_sza2] 

geofile1 = self.geofiles[idx_sza1] 

geofile2 = self.geofiles[idx_sza2] 

config = parse_config(radfile1) 

resfile = str(self.savdir)+'/'+str(self.savstr)+'_case_'+str(config['case'])+'_phi_'+str(config['phi1'])+'_phi0_'+str(config['phi0'])+'_vza_'+str(config['vza'])+'_sza_'+str(sza)+'_ch_'+str(config['cha'])+'.npz' 

picfile = resfile.rsplit('.',1)[0] + '.png' 

if os.path.isfile(resfile): 

print ' .. already calculated, opening: %s ' % resfile 

#img_set = dict(np.load(resfile)) 

tmp_set = np.load(resfile) 

img_set = tmp_set['arr_0'][()] 

else: 

if radfile1!=self.radfile1 or radfile2!=self.radfile2: 

self.radfile1=radfile1 

self.radfile2=radfile2 

self.load_sza(radfile1, geofile1, radfile2, geofile2) 

img_set = self.interp_sza(sza) 

np.savez_compressed(resfile,img_set) 

#scipy.misc.imsave(picfile, np.rot90(img_set['rads'])) 

return img_set 

def interp_sza(self, sza): 

geo_inp1 = calc_cloudsurf(self.dx1, self.dy1, self.dz1, self.conf1['phi1'], 

self.conf1['phi0'], self.conf1['vza'], sza, panorama=self.panorama1) 

self.icumu1 = geo_inp1[0] 

self.icsza1 = np.rad2deg(np.arccos(geo_inp1[1])) 

self.icphi1 = np.rad2deg(np.arccos(geo_inp1[2])) 

self.ibsca1 = np.rad2deg(np.arccos(geo_inp1[3])) 

self.iphi11 = geo_inp1[6] 

self.iphi01 = geo_inp1[7] 

self.ivza1 = geo_inp1[8] 

self.isza1 = geo_inp1[9] 

geo_inp2 = calc_cloudsurf(self.dx2, self.dy2, self.dz2, self.conf2['phi1'], 

self.conf2['phi0'], self.conf2['vza'], sza, panorama=self.panorama2) 

self.icumu2 = geo_inp2[0] 

self.icsza2 = np.rad2deg(np.arccos(geo_inp2[1])) 

self.icphi2 = np.rad2deg(np.arccos(geo_inp2[2])) 

self.ibsca2 = np.rad2deg(np.arccos(geo_inp2[3])) 

self.iphi12 = geo_inp2[6] 

self.iphi02 = geo_inp2[7] 

self.ivza2 = geo_inp2[8] 

self.isza2 = geo_inp2[9] 

rads_inp1 = self.get_from_lut(self.icsza1, self.icumu1, self.icphi1, self.reff1, self.taus1) 

rads_inp2 = self.get_from_lut(self.icsza2, self.icumu2, self.icphi2, self.reff2, self.taus2) 

self.rads1_dif = rads_inp1-self.rads1 

self.rads2_dif = rads_inp2-self.rads2 

#self.rads1_dif = rads1_difs[self.rad1_match_idx] 

#self.rads2_dif = rads2_difs[self.rad2_match_idx] 

self.rads1_dif[np.isnan(self.rads1_dif)] = 0 

self.rads2_dif[np.isnan(self.rads2_dif)] = 0 

self.rads1_dif = gaussian_filter(self.rads1_dif,2)*self.no_shadow 

self.rads2_dif = gaussian_filter(self.rads2_dif,2)*self.no_shadow 

self.weight = np.abs(self.conf1['sza']-sza)/np.abs(self.conf1['sza']-self.conf2['sza']) 

if self.weight>1: 

self.weight=1 

if np.isnan(self.weight): 

self.weight=1 

self.ireff = self.reff1*(1.-self.weight) + (self.reff2)*(self.weight) 

self.itaus = self.taus1*(1.-self.weight) + (self.taus2)*(self.weight) 

self.icumu = self.icumu1*(1.-self.weight) + (self.icumu2)*(self.weight) 

self.icsza = self.icsza1*(1.-self.weight) + (self.icsza2)*(self.weight) 

self.icphi = self.icphi1*(1.-self.weight) + (self.icphi2)*(self.weight) 

self.ibsca = self.ibsca1*(1.-self.weight) + (self.ibsca2)*(self.weight) 

self.iphi1 = self.iphi11*(1.-self.weight) + (self.iphi12)*(self.weight) 

self.iphi0 = self.iphi01*(1.-self.weight) + (self.iphi02)*(self.weight) 

self.ivza = self.ivza1*(1.-self.weight) + (self.ivza2)*(self.weight) 

self.isza = self.isza1*(1.-self.weight) + (self.isza2)*(self.weight) 

rads1_new = self.rad1+self.rads1_dif 

rads2_new = self.rad2+self.rads2_dif 

rads1_new[rads1_new<0.] = 0. 

rads2_new[rads2_new<0.] = 0. 

self.irads = rads1_new*(1.-self.weight) + rads2_new*(self.weight) 

Dataset = {'rads': self.irads, 

'reff': self.ireff, 

'taus': self.itaus, 

'vzas': self.ivza, 

'szas': self.isza, 

'phi1': self.iphi1, 

'phi0': self.iphi0, 

'bsca': self.ibsca, 

'cumu': self.icumu, 

'csza': self.icsza, 

'cphi': self.icphi} 

return Dataset 

def get_from_lut(self, sza, umu, phi, reff, tau): 

idx_sza, idx_umu, idx_phi, idx_reff, idx_tau, weight = self.locate_in_1dlut(sza, umu, phi, reff, tau) 

refl = self.lut_wc_refl[idx_sza, idx_umu, idx_phi, 0, idx_tau, idx_reff] 

self.refl_interp = np.einsum('ijklm...,ijklm...->...',refl,weight) 

self.radiances = self.L(self.refl_interp, sza, self.I0) 

return self.radiances 

def locate_in_1dlut(self, sza, umu, phi, reff, tau): 

itau, ntau, rtau = self.calc_index_weight(tau, self.bins_tau, collect_out=True, no_interpol=False) 

ireff, nreff, rreff = self.calc_index_weight(reff, self.bins_ref, collect_out=True, no_interpol=False) 

iphi, nphi, rphi = self.calc_index_weight(phi, self.bins_phi, collect_out=True, no_interpol=False)