En raison du déménagement des baies serveurs, les services gitlab.irstea.fr et mattermost.irstea.fr seront interrompus le samedi 2 octobre 2021 au matin. Ils devraient revenir à la normale dans la journée.

Commit fa2140c0 authored by Harmel Tristan's avatar Harmel Tristan
Browse files

add use of xarray for spectral rho

parent 2ea143dc
......@@ -119,8 +119,6 @@ for name, raw in df.resample('1H'):
# Rrs, rho = awr.process_wrapper(wl_common, raw, raw.sza, ws=ws, vza = [vza] * N, azi=raw.azi)
#
#
#
# awr.get_rho_values(raw.sza-90,[vza]*N, raw.azi.values, ws=[2])
#
# rho = awr.get_rho_mobley(awr.rhoM1999,raw.sza-90,[vza]*N, azi=raw.azi, ws=[2])[0,:,0]
......@@ -20,11 +20,16 @@ import plotly.graph_objects as go
import plotly.express as px
import plotly.offline as po
opb = os.path.basename
opj = os.path.join
idir = '/DATA/projet/hydrosim/data/trios'
odir = opj(idir, 'L2/awr')
figdir = opj(idir, 'fig/L2')
coordf = os.path.join(idir, '../metadata', 'Coordenadas_Rrs__2016_2019.xlsx')
metaf = os.path.join(idir, 'Metadata_hydroSim.xlsx')
metaf = os.path.join(idir, 'Metadata_hydroSim_awr.xlsx')
infofile = os.path.join(idir, 'triosdatainfo.csv')
L1dir = os.path.join(idir, 'L1')
......@@ -41,74 +46,63 @@ def get_meta(file):
def generate_database_info(idirs, infofile):
opb = os.path.basename
meta_attrs = ['%IDDevice', '%IntegrationTime', '%IDDataBack', '%IDDataCal', '%CalFactor']
df = pd.DataFrame(columns=np.concatenate([['ID', 'date', 'lat', 'lon', 'Ed_file', 'Lsky_file', 'Lt_file'],
df = pd.DataFrame(columns=np.concatenate([['ID', 'date', 'lat', 'lon', 'wind', 'vza', 'azi',
'Ed_file', 'Lsky_file', 'Lt_file'],
['Ed' + s for s in meta_attrs],
['Lsky' + s for s in meta_attrs],
['Lt' + s for s in meta_attrs]]))
# loop on directories (for each date)
for i, meta in metas.iterrows():
#print(i)
# print(i)
date_dir = meta[0].strftime('%Y%m%d')
full_date = dt.datetime.combine(meta[0], meta[1])
ID, site, comment, lat, lon, wind = meta[2:8]
ID, site, comment, lat, lon, wind, vza, azi = meta[2:10]
idir_ = os.path.join(idir, 'L1', date_dir)
date = opb(idir_)
print(idir_ + '/' + ID.lower() + ID.replace('id',''))
print(idir_ + '/' + ID.lower() + ID.replace('id', ''))
files = pd.Series(
np.unique(
np.append(glob.glob(idir_ + '/awr/' + ID.lower() + '_*.mlb'),
glob.glob(idir_ + '/awr/ID' + ID.replace('id','').replace('ID','') + '_*.mlb'))))
print('nb files ',len(files))
np.unique(
np.append(glob.glob(idir_ + '/awr/' + ID.lower() + '_*.mlb'),
glob.glob(idir_ + '/awr/ID' + ID.replace('id', '').replace('ID', '') + '_*.mlb'))))
if len(files) == 0:
continue
print('nb files ', len(files))
for file in files:
name = opb(file).replace('.mlb','')
name = opb(file).replace('.mlb', '')
print(name.split('_'))
ID = name.split('_')[0]
Edf = opj(idir_, 'awr', ID + '_Ed.mlb')
Lskyf = opj(idir_, 'awr', ID + '_Ld.mlb')
Ltf = opj(idir_, 'awr', ID + '_Lu.mlb')
# loop on data files (for each acquisition sequence)
for Edf in files:
ID_ = ID.lower() #Edf.split('_')[-1].replace('.mlb', '')
pattern = Edf.split('_')[0:-2]
try:
Lskyf = glob.glob(Edf.replace('Ed', 'Ld').replace('.mlb','*.csv'))[0]
except:
print('problem with: ', date, ID_, pattern, Edf)
continue
try:
Ltf = glob.glob(Edf.replace('Ed', 'Lu').replace(Edf.split('_')[-2], '*'))[0]
except:
print('problem with: ', date, ID_, pattern, Edf)
continue
# read metadata from files
Edmeta, Lskymeta, Ltmeta = get_meta(Edf), get_meta(Lskyf), get_meta(Ltf)
# save into dataframe df
info = np.concatenate([[ID_, date, lat, lon, opb(Edf), opb(Lskyf), opb(Ltf)], \
Edmeta[meta_attrs].values[0], Lskymeta[meta_attrs].values[0],
Ltmeta[meta_attrs].values[0]])
df.loc[i] = info
i += 1
# read metadata from files
Edmeta, Lskymeta, Ltmeta = get_meta(Edf), get_meta(Lskyf), get_meta(Ltf)
# save into dataframe df
info = np.concatenate([[ID, date, lat, lon, wind, vza, azi, opb(Edf), opb(Lskyf), opb(Ltf)], \
Edmeta[meta_attrs].values[0], Lskymeta[meta_attrs].values[0],
Ltmeta[meta_attrs].values[0]])
df.loc[i] = info
df.to_csv(infofile, index=False)
def plot_map(coords, datafile=''):
def plot_map( idir, datafile=''):
''' Plot for band 'wl' and 'method' '''
wl = 659
wl = 860
#wl = 860
method = 'M99'
if datafile == '':
datafile = os.path.join(idir, 'L2/all/Rrs_stat.csv')
datafile = opj(idir, 'Rrs_stat.csv')
df = pd.read_csv(datafile, parse_dates=True)
df['datetime'] = df['date']
......@@ -137,9 +131,9 @@ def plot_map(coords, datafile=''):
size='size',
range_color=[0, df_['0.5'].max() * 0.7],
animation_frame='month',
color_continuous_scale=px.colors.cyclical.IceFire,
zoom=5.5, opacity=0.45, height=900, width=1100,
title='Sample points from Rogerio; Rrs at ' + str(wl) + ' nm')
color_continuous_scale=px.colors.diverging.RdYlBu_r,
zoom=5.5, opacity=0.65, height=900, width=1100,
title='Sample points; Rrs at ' + str(wl) + ' nm')
# fig.update_traces(marker=dict(line_width=2),selector=dict(mode='markers'))
fig.update_traces(
......@@ -162,7 +156,7 @@ def plot_map(coords, datafile=''):
])
# fig.update_layout(margin={"r":10,"t":0,"l":10,"b":0})
po.plot(fig, filename=os.path.join(idir, '../fig/map_rogerio_Rrs' + str(wl) + '.html'))
po.plot(fig, filename=opj(figdir, 'map_Rrs' + str(wl) + '.html'))
def call(command):
......@@ -177,12 +171,10 @@ def call(command):
def process(infofile, method='M99', ncore=10):
odir = '/DATA/OBS2CO/data/rogerio/data/L2'
figdir = '/DATA/OBS2CO/data/rogerio/fig/L2'
info = pd.read_csv(infofile)
utc_conv = '4'
command = []
commands = []
for i, info_ in info.iterrows():
......@@ -195,24 +187,25 @@ def process(infofile, method='M99', ncore=10):
if not os.path.exists(odir_):
os.makedirs(odir_)
lat, lon = str(info_.lat), str(info_.lon)
lat, lon,vza,azi,wind = str(info_.lat), str(info_.lon),str(info_.vza), str(info_.azi),str(info_.wind)
Edf, Lskyf, Ltf = info_.Ed_file, info_.Lsky_file, info_.Lt_file
command.append('trios_processing ' + os.path.join(L1dir, date) + ' ' + ID + \
' awr --lat ' + lat + ' --lon ' + lon + ' --name _' + method + \
' --data_files "' + Edf + ' ' + Lskyf + ' ' + Ltf + \
'" --odir ' + odir_ + ' --utc_conv=' + utc_conv + \
' --method ' + method + ' --plot --figdir ' + figdir + ' --format mlb --no_clobber')
command = 'trios_processing ' + os.path.join(L1dir, date,'awr') + ' ' + ID + \
' awr --lat ' + lat + ' --lon ' + lon + ' --name _' + method + \
' --vza ' + vza + ' --azi ' + azi + ' --ws ' + wind + \
' --data_files "' + Edf + ' ' + Lskyf + ' ' + Ltf + \
'" --odir ' + odir_ + ' --utc_conv=' + utc_conv + \
' --method ' + method + ' --plot --figdir ' + figdir + ' --format mlb --no_clobber'
commands.append(command)
print(command)
with Pool(processes=ncore) as pool:
print(pool.map(call, command))
print(pool.map(call, commands))
pool.close
def post_process(create_stat_file=False):
idir = '/DATA/OBS2CO/data/rogerio/data/L2'
def post_process(idir,figdir,create_stat_file=False):
stat_file = os.path.join(idir, 'all/Rrs_stat.csv')
if create_stat_file:
......@@ -240,14 +233,14 @@ def post_process(create_stat_file=False):
df.sort_values(['date', 'ID', 'method', 'wl'], inplace=True)
import matplotlib.pyplot as plt
plt.ioff()
ncols = 4
obj = df.groupby(['date', 'ID'])
N = obj.groups.__len__()
rows = N // ncols + (1 if N % ncols else 0)
aspect_ratio = 1 * ncols
fig, axs = plt.subplots(nrows=rows, ncols=ncols, figsize=(20, rows * aspect_ratio))
aspect_ratio = 1.3 * ncols
fig, axs = plt.subplots(nrows=rows, ncols=ncols, figsize=(25, rows * aspect_ratio))
fig.subplots_adjust(left=0.1, right=0.9, hspace=.5, wspace=0.45)
c = dict(M99='red', osoaa='black', temp_opt='blue')
......@@ -260,7 +253,8 @@ def post_process(create_stat_file=False):
ax.legend(loc='best', frameon=False)
ax.set_xlabel(r'Wavelength (nm)')
ax.set_title('ID ' + idx[1] + ', ' + idx[0].date().__str__())
plt.savefig('/DATA/OBS2CO/data/rogerio/fig/compare_method.pdf', bbox_inches='tight')
plt.tight_layout()
plt.savefig(figdir, bbox_inches='tight',dpi=300)
plt.close()
......@@ -276,7 +270,7 @@ process(infofile, method=method, ncore=10)
for method in ['M99', 'osoaa', 'temp_opt']:
process(infofile, method=method, ncore=10)
post_process(create_stat_file=True)
post_process(odir,figdir,create_stat_file=True)
# ----
# END
......@@ -60,7 +60,7 @@ irr = ua.irradiance()
irr.load_F0()
# TODO check noise values (e.g., NEI from Trios), should it be spectral?
noise = 0.1
idpr = '169'
idpr = '150'
idx_list_for_plot=(28, 37, 51, 71, 91, 105, 130, 140, 170)
#idx_list_for_plot=(165,170,175,180,185,190,195,200,205)
......@@ -109,13 +109,15 @@ for idpr in idprs[0:]:
mean = df.groupby('rounded_depth').mean()
median = df.groupby('rounded_depth').median()
std = df.groupby('rounded_depth').std()
q25 = df.groupby('rounded_depth').quantile(0.25)
q75 = df.groupby('rounded_depth').quantile(0.75)
df_=df.drop(df.columns[df.dtypes=='object'],axis=1)
q25 = df_.groupby('rounded_depth').quantile(0.25)
q75 = df_.groupby('rounded_depth').quantile(0.75)
# ---------------------------
# Data processing
# ---------------------------
# load process object
process = iwr_process(wl=wl_)
# process data
......
......@@ -32,7 +32,7 @@ class awr_process:
self.M1999_file = M1999_file
self.M2015_file = M2015_file
self.load_rho_lut()
self.rho = self.rhosoaa_fine
self.rho = self.rhosoaa_coarse
def load_rho_lut(self):
self.rhosoaa_fine = pd.read_csv(self.rhosoaa_fine_file, index_col=[0, 1, 2, 3, 4, 5])
......@@ -83,7 +83,7 @@ class awr_process:
:return:
'''
# Warning SZA set to 90 if Sun below the horizon
#sza[sza>90]=90
# sza[sza>90]=90
rhodf = rhodf.query('sza<75 & vza >0')
rhodf.index = rhodf.index.remove_unused_levels()
......@@ -104,7 +104,7 @@ class awr_process:
# filtering
# ------------------
ind = self.filtering(self.df.Lt, self.df.Lsky, self.df.Ed)
clean = self.df #[ind]
clean = self.df # [ind]
Lt, Lsky, Ed, sza = clean.Lt.values, clean.Lsky.values, clean.Ed.values, clean.sza.values
# -----------------------------
......@@ -140,11 +140,12 @@ class awr_process:
Ltm = Lt.mean(axis=0)
Edm = Ed.mean(axis=0)
def add_envelope(ax,wl,values,label='',**kwargs):
up.add_curve(ax, wl, values.mean(axis=0), label=label, c='black',**kwargs)
ax.fill_between(wl,np.quantile(values,0.05, axis=0),np.quantile(values,0.95, axis=0), alpha=0.25, color='grey')
ax.fill_between(wl,np.quantile(values,0.25, axis=0),np.quantile(values,0.75, axis=0), alpha=0.35, color='red')
def add_envelope(ax, wl, values, label='', **kwargs):
up.add_curve(ax, wl, values.mean(axis=0), label=label, c='black', **kwargs)
ax.fill_between(wl, np.quantile(values, 0.05, axis=0), np.quantile(values, 0.95, axis=0), alpha=0.25,
color='grey')
ax.fill_between(wl, np.quantile(values, 0.25, axis=0), np.quantile(values, 0.75, axis=0), alpha=0.35,
color='red')
mpl.rcParams.update({'font.size': 18})
fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(20, 12))
......@@ -152,26 +153,26 @@ class awr_process:
# ---- Ed
ax = axs[0, 0]
add_envelope(ax,wl,Ed,label=r'$E_s$')
#up.add_curve(ax, wl, Ed.mean(axis=0), Ed.std(axis=0), label=r'$E_s$', c='black')
add_envelope(ax, wl, Ed, label=r'$E_s$')
# up.add_curve(ax, wl, Ed.mean(axis=0), Ed.std(axis=0), label=r'$E_s$', c='black')
ax.set_ylabel(r'$E_{d}(0^{+})\ (mW\ m^{-2}\ nm^{-1})$')
ax.set_xlabel(r'Wavelength (nm)')
ax.legend(loc='best', frameon=False)
# ---- Lsky
ax = axs[0, 1] #ax2 = ax.twinx()
add_envelope(ax,wl,Lsky,label=r'$L_{sky}$')
ax = axs[0, 1] # ax2 = ax.twinx()
add_envelope(ax, wl, Lsky, label=r'$L_{sky}$')
ax.set_ylabel(r'$L_{sky}\ (mW\ m^{-2}\ nm^{-1}\ sr^{-1})$')#, color='r')
#ax.tick_params('y', colors='r')
ax.set_ylabel(r'$L_{sky}\ (mW\ m^{-2}\ nm^{-1}\ sr^{-1})$') # , color='r')
# ax.tick_params('y', colors='r')
ax.set_xlabel(r'Wavelength (nm)')
ax.legend(loc='best', frameon=False)
# ---- Lt vs Lsurf
ax = axs[0, 2]
add_envelope(ax,wl,Lt,label=r'$L_t$')
#up.add_curve(ax, wl, Lt.mean(axis=0), Lt.std(axis=0), label=r'$L_t$', c='black')
add_envelope(ax,wl,Lsky* rho,label=r'$L_{surf}$',linestyle=':')
add_envelope(ax, wl, Lt, label=r'$L_t$')
# up.add_curve(ax, wl, Lt.mean(axis=0), Lt.std(axis=0), label=r'$L_t$', c='black')
add_envelope(ax, wl, Lsky * rho, label=r'$L_{surf}$', linestyle=':')
# up.add_curve(ax, wl, Lsky.mean(axis=0) * rho, Lsky.std(axis=0) * rho,
# label=method + ' (' + str(round(rho, 4)) + ')', c='violet')
ax.set_ylabel(r'$L_t\ or\ L_{surf}\ (mW\ m^{-2}\ nm^{-1}\ sr^{-1})$')
......@@ -180,7 +181,7 @@ class awr_process:
# ---- Proportion o(Lt - Lsurf ) /Lt
ax = axs[1, 0]
add_envelope(ax,wl,Lsky * rho / Ltm,label=r'$L_{surf}/L_t$')
add_envelope(ax, wl, Lsky * rho / Ltm, label=r'$L_{surf}/L_t$')
# up.add_curve(ax, wl, Lsky.mean(axis=0) * rho / Ltm, Lsky.std(axis=0) * rho,
# label=method + ' (' + str(round(rho, 4)) + ')', c='violet')
ax.set_ylabel(r'$L_{surf}/L_t$')
......@@ -188,7 +189,7 @@ class awr_process:
# ---- Lw
ax = axs[1, 1]
add_envelope(ax,wl,Rrs* Edm)
add_envelope(ax, wl, Rrs * Edm)
# up.add_curve(ax, wl, Rrs.mean(axis=0) * Edm, Rrs.std(axis=0) * Edm,
# label=method + ' (' + str(round(rho, 4)) + ')', c='blue')
......@@ -197,7 +198,7 @@ class awr_process:
# ---- Rrs
ax = axs[1, 2]
add_envelope(ax,wl,Rrs)
add_envelope(ax, wl, Rrs)
# up.add_curve(ax, wl, Rrs.transpose().mean(axis=1), Rrs.transpose().std(axis=1),
# label=method + ' (' + str(round(rho, 4)) + ')', c='blue')
ax.set_ylabel(r'$R_{rs}\ (sr^{-1})$')
......@@ -442,6 +443,16 @@ class iwr_process:
self.name = name
self.idpr = idpr
################
# load aux data
iopw = ua.iopw()
iopw.load_iopw()
irr = ua.irradiance()
irr.load_F0()
self.aw, self.bbw = iopw.get_iopw(wl)
self.F0 = irr.get_F0(wl)
def call_process(self, ofile="", plot_file=""):
# ---------------------------
......@@ -460,17 +471,17 @@ class iwr_process:
mean = df.groupby('rounded_depth').mean()
median = df.groupby('rounded_depth').median()
std = df.groupby('rounded_depth').std()
df_=df.drop(df.columns[df.dtypes=='object'],axis=1)
df_ = df.drop(df.columns[df.dtypes == 'object'], axis=1)
q25 = df_.groupby('rounded_depth').quantile(0.25)
q75 = df_.groupby('rounded_depth').quantile(0.75)
# ---------------------------
# Data processing
# ---------------------------
res = self.process(mean, std)
res_w = self.process(mean, std, mode='log')
res_med = self.process(median, std)
res_lsq = self.process(mean, std, mode='lsq')
res, res_Lu = self.process(mean, std)
res_w, void = self.process(mean, std, mode='log')
res_med, void = self.process(median, std)
res_lsq, void = self.process(mean, std, mode='lsq')
Kd = res_lsq.popt[:, 0]
Edm = res_lsq.popt[:, 1]
......@@ -512,29 +523,36 @@ class iwr_process:
fig, axs = plt.subplots(nrows=3, ncols=3, figsize=(16, 10))
fig.subplots_adjust(left=0.1, right=0.9, hspace=.5, wspace=0.25)
print(std.prof_Edz)
for idx in idx_list_for_plot:
ax = axs.flat[i]
i += 1
y = mean.Edz.iloc[:, idx]
sigma = std.Edz.iloc[:, idx]
yerr = [median.Edz.iloc[:, idx] - q25.Edz.iloc[:, idx], q75.Edz.iloc[:, idx] - median.Edz.iloc[:, idx]]
ax.errorbar(x, y, yerr=yerr, ecolor='lightgray', alpha=0.43, fmt='none', elinewidth=3, capsize=0,
xerr = std.prof_Edz.values # [mean.prof_Edz - std.prof_Edz, mean.prof_Edz + std.prof_Edz]
yerr = [mean.Edz.iloc[:, idx] - std.Edz.iloc[:, idx], mean.Edz.iloc[:, idx] + std.Edz.iloc[:, idx]]
ax.errorbar(x, y, xerr=xerr, yerr=yerr, ecolor='gray', alpha=0.73, fmt='none', elinewidth=2, capsize=0,
label='std')
ax.scatter(x, y,
# s = std.Edz.iloc[:,50]/std.Edz.iloc[:,50].mean()+100,
c=mean.Ed.iloc[:, idx],
alpha=0.6, cmap=cmocean.cm.thermal, label=None
)
Ed_sim = iwr_process().f_Edz(depth_, *res.popt[idx, :])
Ed_sim = self.f_Edz(depth_, *res.popt[idx, :])
ax.plot(depth_, Ed_sim, linestyle='-', c='black', label='mean')
Ed_sim = iwr_process().f_Edz(depth_, *res_w.popt[idx, :])
Ed_sim = self.f_Edz(depth_, *res_w.popt[idx, :])
ax.plot(depth_, Ed_sim, linestyle=':', c='black', label='log-space')
# Ed_sim = iwr_process.f_Edz(depth_, *res_rw[idx][0])
# Ed_sim = self.f_Edz(depth_, *res_rw[idx][0])
# ax.plot(depth_, Ed_sim, linestyle=':', c='red', label='mean, relativ weighted')
Ed_sim = iwr_process().f_Edz(depth_, *res_med.popt[idx, :])
Ed_sim = self.f_Edz(depth_, *res_med.popt[idx, :])
ax.plot(depth_, Ed_sim, linestyle='--', c='black', label='median')
ax.semilogy()
# ax.semilogx()
# ax.colorbar()
ax.set_ylabel(r'$E_d\ ({mW\cdot m^{-2}\cdot nm^{-1}})$')
ax.set_xlabel('Depth (m)')
......@@ -571,16 +589,18 @@ class iwr_process:
c=mean.Ed.iloc[:, idx],
alpha=0.6, cmap=cmocean.cm.thermal, label=None
)
Lu_sim = iwr_process().f_Lu(depth_, *res_lsq.popt[idx, -2:])
Lu_sim = self.f_Lu(depth_, *res_lsq.popt[idx, -3:-1])
ax.plot(depth_, Lu_sim, linestyle='-', c='black', label='mean')
# Lu_sim = iwr_process().f_Luz(depth_, *res_w.popt[idx, :])
# Lu_sim = self.f_Luz(depth_, *res_w.popt[idx, :])
# ax.plot(depth_, Lu_sim, linestyle=':', c='black', label='log-space')
# # Lu_sim = iwr_process.f_Luz(depth_, *res_rw[idx][0])
# # Lu_sim = self.f_Luz(depth_, *res_rw[idx][0])
# # ax.plot(depth_, Lu_sim, linestyle=':', c='red', label='mean, relativ weighted')
# Lu_sim = iwr_process().f_Luz(depth_, *res_med.popt[idx, :])
# Lu_sim = self.f_Luz(depth_, *res_med.popt[idx, :])
# ax.plot(depth_, Lu_sim, linestyle='--', c='black', label='median')
ax.semilogy()
# ax.semilogx()
# ax.colorbar()
ax.set_ylabel(r'$L_u\ ({mW\cdot m^{-2}\cdot sr^{-1}\cdot nm^{-1}})$')
ax.set_xlabel('Depth (m)')
......@@ -633,7 +653,9 @@ class iwr_process:
iparam = 2
ax = axs.flat[iparam]
y, sigma = res_lsq.popt[:, iparam], np.sqrt(res_lsq.pcov[:, iparam, iparam]) * 0
y, sigma = res_Lu.popt[:, 0], np.sqrt(res_Lu.pcov[:, 0,0])
up.add_curve(ax, wl, y, sigma, c='red', label='mean')
y, sigma = res_lsq.popt[:, iparam], np.sqrt(res_lsq.pcov[:, iparam, iparam]) * 1e-3
sigma[sigma > 10 * np.nanmedian(sigma)] = np.nan
up.add_curve(ax, wl, y, sigma, c='blue', label='lsq')
ax.set_ylabel(r'$K_{Lu}\ ({m^{-1}})$')
......@@ -642,10 +664,13 @@ class iwr_process:
iparam = 3
ax = axs.flat[iparam]
y, sigma = res_lsq.popt[:, iparam], np.sqrt(res_lsq.pcov[:, iparam, iparam]) * 0
y, sigma = res_lsq.popt[:, iparam], np.sqrt(res_lsq.pcov[:, iparam, iparam]) * 1e-3
sigma[sigma > 10 * np.nanmedian(sigma)] = np.nan
up.add_curve(ax, wl, y, sigma, c='blue', label='lsq')
ax.set_ylabel(r'$L_{w}(0^{-})\ ({mW\cdot m^{-2}\cdot sr^{-1}\cdot nm^{-1})$')
y, sigma = res_Lu.popt[:, 1], np.sqrt(res_Lu.pcov[:, 1, 1])
up.add_curve(ax, wl, y, sigma, c='red', label='mean')
ax.legend(loc='best', frameon=False)
ax.set_ylabel(r'$L_{w}(0^{-})\ (mW\cdot m^{-2}\cdot sr^{-1}\cdot nm^{-1})$')
ax.set_xlabel(r'Wavelength (nm)')
fig.suptitle('trios_iwr ' + self.name + ' idpr' + self.idpr, fontsize=16)
......@@ -681,7 +706,20 @@ class iwr_process:
ax.set_xlabel(r'Wavelength (nm)')
ax.legend(loc='best', frameon=False)
ax = axs.flat[3]
iparam = 0
y, sigma = res.popt[:, iparam], res.pcov[:, iparam, iparam]
ax.plot(wl, np.gradient(np.gradient(y, wl), wl), c='red', label='mean')
y, sigma = res_w.popt[:, iparam], res_w.pcov[:, iparam, iparam]
ax.plot(wl, np.gradient(np.gradient(y, wl), wl), c='orange', label='log-space')
y, sigma = res_med.popt[:, iparam], res_med.pcov[:, iparam, iparam]
ax.plot(wl, np.gradient(np.gradient(y, wl), wl), c='green', label='median')
y, sigma = res_lsq.popt[:, iparam], np.sqrt(res_lsq.pcov[:, iparam, iparam])
ax.plot(wl, np.gradient(np.gradient(y, wl), wl), c='blue', label='lsq')
ax.set_xlabel(r'Wavelength (nm)')
ax.set_ylabel(r'$\delta K_{d}/\delta\lambda\ ({m^{-1}})$')
fig.suptitle('trios_iwr ' + self.name + ' idpr' + self.idpr, fontsize=16)
pdf.savefig()
......@@ -691,60 +729,72 @@ class iwr_process:
def process(self, meas, std, mode='linear'):
wl_ = self.wl
################
# load aux data
iopw = ua.iopw()
iopw.load_iopw()
irr = ua.irradiance()
irr.load_F0()
# TODO check noise values (e.g., NEI from Trios), should it be spectral?
noise = 0.1
N = len(wl_)
x = meas.prof_Edz # - 0.56
res = uu.fit(N)
res_Lu = uu.fit(N)
if mode == 'lsq':
res = uu.fit(N, 4)
res = uu.fit(N, 5)
first =True
for idx, wl in enumerate(wl_[:-10]):
aw, bbw = iopw.get_iopw(wl)
F0 = irr.get_F0(wl)
aw, bbw = self.aw[idx], self.bbw[idx]
F0 = self.F0[idx]
y = meas.Edz.iloc[:, idx]
sigma = std.Edz.iloc[:, idx]
sigma[sigma < noise] = noise
sigma.fillna(np.inf, inplace=True)
if mode == 'linear':
res.popt[idx, :], res.pcov[idx, ...] = so.curve_fit(self.f_Edz, x, y, [1.1 * aw, 100],
res.popt[idx, :], res.pcov[idx, ...] = so.curve_fit(self.f_Edz, x, y, [1.1 * aw, 0.9*F0],
bounds=([aw, 0], [np.inf, F0]))
xx = meas.prof_Luz
yy = meas.Luz.iloc[:, idx]
res_Lu.popt[idx, :], res_Lu.pcov[idx, ...] = so.curve_fit(self.f_Lu, xx, yy, [1.1 * aw, 0.005],
bounds=([aw/2, 0], [np.inf, np.inf]))
elif mode == 'log':
res.popt[idx, :], res.pcov[idx, ...] = so.curve_fit(self.f_logEdz, x, np.log(1 + y),
[1.1 * aw, 100], bounds=(
[aw, 0], [np.inf, F0])) # , sigma=sigma, absolute_sigma=True
elif mode == 'lsq':
z = (meas.prof_Edz, meas.prof_Luz)
y = (meas.Edz.iloc[:, idx], meas.Luz.iloc[:, idx])
y = (meas.Edz.iloc[:, idx], meas.Luz.iloc[:, idx], meas.Ed.iloc[:, idx].mean())
sig_Edz = self.format_sigma(std.Edz.iloc[:, idx], meas.Edz.iloc[:, idx], 0.1)
sig_Luz = self.format_sigma(std.Luz.iloc[:, idx], meas.Luz.iloc[:, idx], 1e-3)
sigma = (sig_Edz, sig_Luz)
sigma = (1, 1)
x0 = [1.1 * aw, min(F0, meas.Ed.iloc[:, idx].mean()), 1.1 * aw, meas.Luz.iloc[0, idx]]
sigma = (meas.Ed.iloc[:, idx].values,meas.Ed.iloc[:, idx].values)
#sigma = (1, 1)
x0 = [1.1 * aw, min(F0, meas.Ed.iloc[:, idx].mean()), 1.1 * aw, meas.Luz.iloc[0, idx], 0]
print('x0', wl, F0, aw, x0)
print(([aw, 0, aw / 2, 0,-0.1], [np.inf, F0, np.inf, np.inf,0.1]))
lsq = so.least_squares(self.cost_func, x0, args=(z, y, sigma),