Skip to content

Commit

Permalink
pySICEv1.3 (#33)
Browse files Browse the repository at this point in the history 
 Latest update of python scripts: 20-04-2020 ([email protected])
 From Baptiste:
- reorganized sice_lib.py
- prevented code to crash when no pixels are suitable for retrieval
  • Loading branch information
@BaptisteVandecrux
BaptisteVandecrux authored Apr 20, 2020
1 parent c0c3ef7 commit 44f2cf2
Show file tree
Hide file tree
Showing 2 changed files with 265 additions and 487 deletions.
167 changes: 83 additions & 84 deletions sice.py
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,12 @@
# pySICEv1.2
# pySICEv1.3
#
# from FORTRAN VERSION 5
# March 31, 2020
#
# Update 16-04-2020 ([email protected])
# Latest update of python script: 20-04-2020 ([email protected])
# From Baptiste:
#- prevented caluclation on nan values in snow_impurities function
#- now use clean snow BBA caluclation for clear plluted pixels (isnow = 7)
# From Alex's side:
#- approximation of the bba clean snow
#- reorganized sice_lib.py
#- prevented code to crash when no pixels are suitable for retrieval
#
# BAV 10-10-2019 ([email protected])
# Changes:
Expand Down Expand Up @@ -104,7 +102,7 @@

start_time = time.process_time()

InputFolder = sys.argv[1] + '/'
InputFolder = sys.argv[1] + '/'

#%% ========= input tif ================
Oa01 = rio.open(InputFolder+'r_TOA_01.tif')
Expand Down Expand Up @@ -157,18 +155,18 @@ def WriteOutput(var,var_name,in_folder):
isnow[sza>75] = 100
isnow[toa_cor_o3[20, :,:] < 0.1] = 102
for i_channel in range(21):
toa_cor_o3[i_channel, np.logical_not(np.isnan(isnow))] = np.nan
toa_cor_o3[i_channel, ~np.isnan(isnow)] = np.nan

vaa[np.logical_not(np.isnan(isnow))] = np.nan
saa[np.logical_not(np.isnan(isnow))] = np.nan
sza[np.logical_not(np.isnan(isnow))] = np.nan
vza[np.logical_not(np.isnan(isnow))] = np.nan
vaa[ ~np.isnan(isnow)] = np.nan
saa[ ~np.isnan(isnow)] = np.nan
sza[ ~np.isnan(isnow)] = np.nan
vza[ ~np.isnan(isnow)] = np.nan
height = height.astype(float)
height[np.logical_not(np.isnan(isnow))] = np.nan
height[ ~np.isnan(isnow)] = np.nan

# =========== view geometry and atmosphere propeties ==============
raa, am1, am2, ak1, ak2, amf, co = sl.view_geometry(vaa, saa, sza, vza, aot, height)
tau, p, g = sl.aerosol_properties(aot, height, co)
tau, p, g,gaer,taumol,tauaer = sl.aerosol_properties(aot, height, co)

# =========== snow properties ====================================
D, area, al, r0, bal = sl.snow_properties(toa_cor_o3, ak1, ak2)
Expand All @@ -189,7 +187,8 @@ def WriteOutput(var,var_name,in_folder):
# =========== clean snow ====================================
# for that we calculate the theoretical reflectance at band 1 of a surface with:
# r0 = 1, a (albedo) = 1, ak1 = 1, ak2 = 1
t1, t2, ratm, r, astra, rms = sl.prepare_coef(tau, g, p, am1, am2, amf)
# t1 and t2 are the backscattering fraction
t1, t2, ratm, r, astra, rms = sl.prepare_coef(tau, g, p, am1, am2, amf,gaer,taumol,tauaer)
rs_1 = sl.alb2rtoa(1, t1[0,:,:], t2[0,:,:], np.ones_like(r0), np.ones_like(ak1),
np.ones_like(ak2), ratm[0,:,:], r[0,:,:])

Expand All @@ -202,12 +201,12 @@ def WriteOutput(var,var_name,in_folder):
def mult_channel(c,A):
tmp = A.T*c
return tmp.T

alb_sph = np.exp(-np.sqrt(1000.*4.*np.pi* mult_channel(bai/w, np.tile(al,(21,1,1)))))
alb_sph[alb_sph>0.999]=1

# ========== very dirty snow ====================================
ind_pol = toa_cor_o3[0,:,:] < rs_1

isnow[ind_pol] = 1

ind_very_dark = np.logical_and(toa_cor_o3[20]<0.4, ind_pol)
Expand All @@ -228,75 +227,75 @@ def mult_channel(c,A):

# =========== polluted snow ====================================
ind_pol = np.logical_or(ind_very_dark, ind_pol)
subs_pol = np.argwhere(ind_pol)

# approximation of the transcendental equation allowing closed-from solution
#alb_sph[:,ind_pol] = (toa_cor_o3[:,ind_pol] - r[:,ind_pol])/(t1[:,ind_pol]*t2[:,ind_pol]*r0[ind_pol] + ratm[:,ind_pol]*(toa_cor_o3[:,ind_pol] - r[:,ind_pol]))

# solving iteratively the transcendental equation
alb_sph[:,ind_pol] = 1

def solver_wrapper(toa_cor_o3,tau, t1, t2, r0, ak1, ak2, ratm, r):
def func_solv(albedo):
return toa_cor_o3 - sl.alb2rtoa(albedo, t1, t2, r0, ak1, ak2, ratm, r)
# it is assumed that albedo is in the range 0.1-1.0
return sl.zbrent(func_solv, 0.1, 1, 100, 1.e-6)

solver_wrapper_v = np.vectorize(solver_wrapper)

# loop over all bands except band 19, 20
for i_channel in np.append(np.arange(18), [20]):
alb_sph[i_channel,ind_pol] = solver_wrapper_v(
toa_cor_o3[i_channel,ind_pol],
tau[i_channel,ind_pol],
t1[i_channel,ind_pol],
t2[i_channel,ind_pol],
r0[ind_pol], ak1[ind_pol],
ak2[ind_pol],
ratm[i_channel,ind_pol],
r[i_channel,ind_pol])
if np.any(ind_pol):
subs_pol = np.argwhere(ind_pol)

ind_bad = alb_sph[i_channel,:,:]==-999
alb_sph[i_channel,ind_bad] = np.nan
isnow[ind_bad]= -i_channel

# INTERNal CHECK FOR CLEAN PIXELS
# Are reprocessed as clean
ind_clear_pol1 = np.logical_and(ind_pol, alb_sph[0,:,:]>0.98)
ind_clear_pol2 = np.logical_and(ind_pol, alb_sph[1,:,:]>0.98)
ind_clear_pol = np.logical_or(ind_clear_pol1, ind_clear_pol2)
isnow[ind_clear_pol]= 7
for i_channel in range(21):
alb_sph[i_channel,ind_clear_pol] = np.exp(-np.sqrt(4.*1000.*al[ind_clear_pol] * np.pi * bai[i_channel] / w[i_channel] ))
# approximation of the transcendental equation allowing closed-from solution
#alb_sph[:,ind_pol] = (toa_cor_o3[:,ind_pol] - r[:,ind_pol])/(t1[:,ind_pol]*t2[:,ind_pol]*r0[ind_pol] + ratm[:,ind_pol]*(toa_cor_o3[:,ind_pol] - r[:,ind_pol]))

# solving iteratively the transcendental equation
alb_sph[:,ind_pol] = 1

def solver_wrapper(toa_cor_o3,tau, t1, t2, r0, ak1, ak2, ratm, r):
def func_solv(albedo):
return toa_cor_o3 - sl.alb2rtoa(albedo, t1, t2, r0, ak1, ak2, ratm, r)
# it is assumed that albedo is in the range 0.1-1.0
return sl.zbrent(func_solv, 0.1, 1, 100, 1.e-6)

solver_wrapper_v = np.vectorize(solver_wrapper)
# loop over all bands except band 19, 20
for i_channel in np.append(np.arange(18), [20]):
alb_sph[i_channel,ind_pol] = solver_wrapper_v(
toa_cor_o3[i_channel,ind_pol],
tau[i_channel,ind_pol],
t1[i_channel,ind_pol],
t2[i_channel,ind_pol],
r0[ind_pol], ak1[ind_pol],
ak2[ind_pol],
ratm[i_channel,ind_pol],
r[i_channel,ind_pol])

ind_bad = alb_sph[i_channel,:,:]==-999
alb_sph[i_channel,ind_bad] = np.nan
isnow[ind_bad]= -i_channel

# INTERNal CHECK FOR CLEAN PIXELS
# Are reprocessed as clean
ind_clear_pol1 = np.logical_and(ind_pol, alb_sph[0,:,:]>0.98)
ind_clear_pol2 = np.logical_and(ind_pol, alb_sph[1,:,:]>0.98)
ind_clear_pol = np.logical_or(ind_clear_pol1, ind_clear_pol2)
isnow[ind_clear_pol]= 7
for i_channel in range(21):
alb_sph[i_channel,ind_clear_pol] = np.exp(-np.sqrt(4.*1000.*al[ind_clear_pol] * np.pi * bai[i_channel] / w[i_channel] ))

# re-defining polluted pixels
ind_pol = np.logical_or(isnow==6, isnow==1)

# re-defining polluted pixels
ind_pol = np.logical_or(isnow==6, isnow==1)

#retrieving snow impurities
ntype, bf, conc = sl.snow_impurities(alb_sph, bal)

# alex 09.06.2019
# reprocessing of albedo to remove gaseous absorption using linear polynomial approximation in the range 753-778nm.
# Meaning: alb_sph[12],alb_sph[13] and alb_sph[14] are replaced by a linear interpolation between alb_sph[11] and alb_sph[15]
afirn=(alb_sph[15,ind_pol]-alb_sph[11,ind_pol])/(w[15]-w[11])
bfirn=alb_sph[15,ind_pol]-afirn*w[15]
alb_sph[12,ind_pol] = bfirn + afirn*w[12]
alb_sph[13,ind_pol] = bfirn + afirn*w[13]
alb_sph[14,ind_pol] = bfirn + afirn*w[14]

# BAV 09-02-2020: 0.5 to 0.35
# pixels that are clean enough in channels 18 19 20 and 21 are not affected by pollution, the analytical equation can then be used
ind_ok = np.logical_and(ind_pol, toa_cor_o3[20,:,:]>0.35)
for i_channel in range(17,21):
alb_sph[i_channel,ind_ok] = np.exp(-np.sqrt(4.*1000.*al[ind_ok] * np.pi * bai[i_channel] / w[i_channel] ))
# Alex, SEPTEMBER 26, 2019
# to avoid the influence of gaseous absorption (water vapor) we linearly interpolate in the range 885-1020nm for bare ice cases only (low toa[20])
# Meaning: alb_sph[18] and alb_sph[19] are replaced by a linear interpolation between alb_sph[17] and alb_sph[20]
delx=w[20]-w[17]
bcoef=(alb_sph[20,ind_pol]-alb_sph[17,ind_pol])/delx
acoef=alb_sph[20,ind_pol]-bcoef*w[20]
alb_sph[18,ind_pol] = acoef + bcoef*w[18]
alb_sph[19,ind_pol] = acoef + bcoef*w[19]
#retrieving snow impurities
ntype, bf, conc = sl.snow_impurities(alb_sph, bal)

# alex 09.06.2019
# reprocessing of albedo to remove gaseous absorption using linear polynomial approximation in the range 753-778nm.
# Meaning: alb_sph[12],alb_sph[13] and alb_sph[14] are replaced by a linear interpolation between alb_sph[11] and alb_sph[15]
afirn=(alb_sph[15,ind_pol]-alb_sph[11,ind_pol])/(w[15]-w[11])
bfirn=alb_sph[15,ind_pol]-afirn*w[15]
alb_sph[12,ind_pol] = bfirn + afirn*w[12]
alb_sph[13,ind_pol] = bfirn + afirn*w[13]
alb_sph[14,ind_pol] = bfirn + afirn*w[14]

# BAV 09-02-2020: 0.5 to 0.35
# pixels that are clean enough in channels 18 19 20 and 21 are not affected by pollution, the analytical equation can then be used
ind_ok = np.logical_and(ind_pol, toa_cor_o3[20,:,:]>0.35)
for i_channel in range(17,21):
alb_sph[i_channel,ind_ok] = np.exp(-np.sqrt(4.*1000.*al[ind_ok] * np.pi * bai[i_channel] / w[i_channel] ))
# Alex, SEPTEMBER 26, 2019
# to avoid the influence of gaseous absorption (water vapor) we linearly interpolate in the range 885-1020nm for bare ice cases only (low toa[20])
# Meaning: alb_sph[18] and alb_sph[19] are replaced by a linear interpolation between alb_sph[17] and alb_sph[20]
delx=w[20]-w[17]
bcoef=(alb_sph[20,ind_pol]-alb_sph[17,ind_pol])/delx
acoef=alb_sph[20,ind_pol]-bcoef*w[20]
alb_sph[18,ind_pol] = acoef + bcoef*w[18]
alb_sph[19,ind_pol] = acoef + bcoef*w[19]

# ========= derivation of plane albedo and reflectance ===========
rp = np.power (alb_sph, ak1)
Expand Down Expand Up @@ -356,4 +355,4 @@ def func_solv(albedo):
WriteOutput(rp[i,:,:], 'albedo_spectral_planar_'+str(i+1).zfill(2), InputFolder)
WriteOutput(refl[i,:,:], 'rBRR_'+str(i+1).zfill(2), InputFolder)

print("End SICE.py %s --- %s seconds ---" % (InputFolder, time.process_time() - start_time))
print("End SICE.py %s --- %s CPU seconds ---" % (InputFolder, time.process_time() - start_time))
Loading

0 comments on commit 44f2cf2

Please sign in to comment.