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 bfcbf70f authored by Harmel Tristan's avatar Harmel Tristan
Browse files

implementation of in-water radiometry (optimization scheme Edz)

parent 0c732c28
......@@ -4,11 +4,13 @@ import scipy.optimize as so
import plotly.graph_objs as go
import trios.utils.utils as uu
from trios.utils.utils import reshape as r
import trios.utils.auxdata as ua
from trios.config import *
class awr_process:
def __init__(self, df=None, wl=None):
self.df = df
......@@ -215,7 +217,7 @@ class awr_process:
geom = [sza[i], vza, azi]
meas = [Lt[i], Lsky[i], Ed[i]]
x0 = ws
res_lsq = so.least_squares(self.cost_func, x0, args=(geom, meas, Rrs_bar))
res_lsq = so.least_squares(self.cost_func, x0, bounds=(0,15), args=(geom, meas, Rrs_bar))
res.append(res_lsq)
x_est.append(res_lsq.x[0])
Rrs, rho = self.process(wl, Lt[i], Lsky[i], Ed[i], sza[i], ws=res_lsq.x[0], azi=azi)
......@@ -312,30 +314,73 @@ class iwr_process:
self.df = df
self.wl = wl
def process(self):
wl = self.wl
df = self.df
def process(self, meas, std, mode ='linear'):
wl_ = self.wl
reflectance = df.loc[:, ("Luz")] / df.loc[:, ("Edz")]
reflectance.columns = pd.MultiIndex.from_product([['reflectance'], reflectance.columns], names=['param', 'wl'])
self.reflectance = reflectance
################
# 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)
for idx, wl in enumerate(wl_[:-10]):
aw, bbw = iopw.get_iopw(wl)
F0 = irr.get_F0(wl)
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],
bounds=([aw, 0], [np.inf, F0]))
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 == 'lmsq':
z = (meas.prof_Edz, meas.prof_Luz)
y = (meas.Edz.iloc[:, idx], meas.Luz.iloc[:, idx])
x0 = [1.1 * aw, 100, 1.1 * aw, meas.Luz.iloc[0, idx]]
so.least_squares(self.cost_func, x0, args=(z, y))
return res
#@staticmethod
def f_Edz(self, depth, Kd, Ed0):
'''simple Edz model for homogeneous water column'''
return Ed0 * np.exp(-Kd * depth)
df['rounded_depth', ''] = df.prof_Edz.round(1)
df.groupby('rounded_depth').mean()
#@staticmethod
def f_logEdz(self, depth, Kd, Ed0):
'''simple Edz model for homogeneous water column'''
return np.log(1 + self.f_Edz(depth, Kd, Ed0)) # Ed0) -Kd*depth
# TODO add central part of run_iwr.py here
def f_Lu(self, depth, KLu, Lw0minus):
'''simple Edz model for homogeneous water column'''
return Lw0minus * np.exp(-KLu * depth)
def cost_func(self, x, z, mes):
z_Edz = z[0]
z_Lu = z[1]
Edz = mes[0]
Lu = mes[1]
cost_f1 = Edz - self.f_Edz( z_Edz, x[0], x[1])
cost_f2 = Lu - self.f_Lu( z_Lu, x[2], x[3])
return np.append(cost_f1,cost_f2)
return self.reflectance
@staticmethod
def f_Edz(depth, Kd, Ed0):
'''simple Edz model for homogeneous water column'''
return Ed0 * np.exp(-Kd * depth)
@staticmethod
def f_logEdz(depth, Kd, Ed0):
'''simple Edz model for homogeneous water column'''
return np.log(1 + iwr_process.f_Edz(depth, Kd, Ed0)) # Ed0) -Kd*depth
def Kd(self, depth, Edz):
Kd = np.diff(Edz) / np.diff(depth)
......
import glob
import re
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 18})
from scipy.interpolate import interp1d
from trios.utils.sunposition import sunpos
......@@ -18,6 +19,7 @@ swrfiles = glob.glob("/DATA/OBS2CO/data/trios/raw/Lu0*idpr*.csv")
iopw = ua.iopw()
iopw.load_iopw()
def add_curve(ax, x, mean, std, c='red', label=''):
ax.plot(x, mean, linestyle='solid', c=c, lw=2.5,
alpha=0.8, label=label)
......@@ -30,7 +32,7 @@ idpr = '167'
# get idpr numbers
idprs = np.unique([re.findall(r'idpr(\d+)', x)[0] for x in swrfiles])
#idprs = np.array(['170'])
# idprs = np.array(['170'])
# loop over idpr
for idpr in idprs:
c = coords[coords.ID_prel == int(idpr)] # .values[0]
......@@ -58,7 +60,7 @@ for idpr in idprs:
Rrs_swr = swr.call_process()
add_curve(ax, wl_swr, Rrs_swr.transpose().mean(axis=1), Rrs_swr.transpose().std(axis=1), label='swr', c='black')
Rrs_swr = swr.call_process(shade_corr=True)
add_curve(ax, wl_swr, Rrs_swr.transpose().mean(axis=1), Rrs_swr.transpose().std(axis=1), label='swr', c='red')
# add_curve(ax, wl_swr, Rrs_swr.transpose().mean(axis=1), Rrs_swr.transpose().std(axis=1), label='swr', c='red')
# -----------------------------------------------
# AWR processing
......@@ -102,37 +104,57 @@ for idpr in idprs:
# Lt.reset_index(level=[1,2],inplace=True)
# merge sensor data on time
df = pd.merge_asof(Lt, Ed, left_index=True, right_index=True, tolerance=pd.Timedelta("2 seconds"),
direction="nearest")
df = pd.merge_asof(df, Lsky, left_index=True, right_index=True, tolerance=pd.Timedelta("2 seconds"),
direction="nearest")
raw = pd.merge_asof(Lt, Ed, left_index=True, right_index=True, tolerance=pd.Timedelta("2 seconds"),
direction="nearest")
raw = pd.merge_asof(raw, Lsky, left_index=True, right_index=True, tolerance=pd.Timedelta("2 seconds"),
direction="nearest")
# add solar angle data and idpr
# compute solar angle (mean between fisrt and last aqcuisition time
df['sza', ''] = np.nan
for index, row in df.iterrows():
raw['sza', ''] = np.nan
for index, row in raw.iterrows():
# print index
sza = sunpos(index, lat, lon, alt)[1]
df.at[index, 'sza'] = sza
rho_h = awr.get_rho_values([df.sza.mean()], [vza], [azi], wl=wl)
rho15 = awr.get_rho_mobley(awr.rhoM2015, [df.sza.mean()], [vza], [azi], [ws])
rho99 = awr.get_rho_mobley(awr.rhoM1999, [df.sza.mean()], [vza], [azi], [ws])
Rrs_h = (df.loc[:, 'Lt'] - rho_h * df.loc[:, 'Lsky']) / df.loc[:, 'Ed']
Rrs15 = (df.loc[:, 'Lt'] - rho15 * df.loc[:, 'Lsky']) / df.loc[:, 'Ed']
Rrs99 = (df.loc[:, 'Lt'] - rho99 * df.loc[:, 'Lsky']) / df.loc[:, 'Ed']
# plt.figure()
raw.at[index, 'sza'] = sza
# ------------------
# filtering
# ------------------
ind = awr.filtering(raw.Lt, raw.Lsky, raw.Ed)
clean = raw[ind]
Lt, Lsky, Ed, sza = clean.Lt.values, clean.Lsky.values, clean.Ed.values, clean.sza.values
# -----------------------------
# data processing
# -----------------------------
Rrs99, rho99 = awr.process_wrapper(wl, clean, clean.sza, ws=ws, azi=azi)
Rrs15, rho15 = awr.process_wrapper(wl, clean, clean.sza, ws=ws, azi=azi, method='M15')
Rrs_h, rho_h = awr.process_wrapper(wl, clean, clean.sza, ws=ws, azi=azi, method='osoaa')
Rrs_opt, Rrs_opt_std = awr.process_optimization(wl, Lt, Lsky, Ed, sza, azi=azi)
wl = Rrs99.T.index.get_level_values(1)
date = Rrs99.index.get_level_values(0).date[0].__str__()
# ------------------
# plotting
# ------------------
# rho_h = awr.get_rho_values([df.sza.mean()], [vza], [azi], wl=wl)
# rho15 = awr.get_rho_mobley(awr.rhoM2015, [df.sza.mean()], [vza], [azi], [ws])
# rho99 = awr.get_rho_mobley(awr.rhoM1999, [df.sza.mean()], [vza], [azi], [ws])
#
# Rrs_h = (df.loc[:, 'Lt'] - rho_h * df.loc[:, 'Lsky']) / df.loc[:, 'Ed']
# Rrs15 = (df.loc[:, 'Lt'] - rho15 * df.loc[:, 'Lsky']) / df.loc[:, 'Ed']
# Rrs99 = (df.loc[:, 'Lt'] - rho99 * df.loc[:, 'Lsky']) / df.loc[:, 'Ed']
add_curve(ax, wl, Rrs15.transpose().mean(axis=1), Rrs15.transpose().std(axis=1),
label='M2015 (' + str(rho15) + ')')
label='M2015 (' + str(round(rho15,4)) + ')',c='violet')
add_curve(ax, wl, Rrs99.transpose().mean(axis=1), Rrs99.transpose().std(axis=1), c='orange',
label='M1999(' + str(rho99) + ')')
add_curve(ax, wl, Rrs_h.transpose().mean(axis=1), Rrs_h.transpose().std(axis=1), c='grey',
label='h(' + str(rho_h.mean()) + ')')
ax.set_title('azi=' + str(azi) + ', vza=' + str(vza) + ', sza=' + str(sza))
label='M1999(' + str(round(rho99,4)) + ')')
add_curve(ax, wl, Rrs_h.transpose().mean(axis=1), Rrs_h.transpose().std(axis=1), c='green',
label='h(' + str(round(rho_h.mean(),4)) + ')')
add_curve(ax, wl, Rrs_opt, Rrs_opt_std, c='blue',
label='Optimization')
ax.set_title('azi=' + str(azi) + ', vza=' + str(vza) + ', sza=' + str(round(sza.mean(),2)))
ax.legend(loc='best', frameon=False)
......@@ -140,4 +162,3 @@ for idpr in idprs:
ax.set_xlabel(r'Wavelength (nm)')
fig.savefig(os.path.join(dirfig, 'trios_awr_' + name + '_idpr' + idpr + '.png'))
plt.close()
......@@ -120,7 +120,6 @@ for idpr in idprs:
wl = Rrs.T.index.get_level_values(1)
date = Rrs.index.get_level_values(0).date[0].__str__()
# ------------------
# plotting
# ------------------
......
......@@ -8,6 +8,8 @@ import matplotlib.pyplot as plt
import cmocean
import plotly
import plotly.graph_objs as go
import scipy.optimize as so
import trios.utils.utils as u
import trios.utils.auxdata as ua
from trios.process import *
......@@ -53,7 +55,7 @@ irr = ua.irradiance()
irr.load_F0()
# TODO check noise values (e.g., NEI from Trios), should it be spectral?
noise = 0.1
idpr = '178'
idpr = '177'
# loop over idpr
for idpr in idprs: # [-1:]:
# idpr=idprs[2]
......@@ -69,7 +71,7 @@ for idpr in idprs: # [-1:]:
ofile = os.path.join(odir, 'Rrs_swr_' + date.values[0] + '_idpr' + idpr + '_' + name + '.csv')
header = c.stack()
header.index = header.index.droplevel()
header.to_csv(ofile, header=None)
header.to_csv(ofile, header=False)
dff = pd.DataFrame()
......@@ -78,42 +80,37 @@ for idpr in idprs: # [-1:]:
# -----------------------------------------------
iwr = u.iwr_data(idpr, iwrfiles)
if iwr.file:
df, wl_ = iwr.reader(lat, lon, alt)
reflectance = iwr_process(df, wl_).process()
df = pd.concat([df, reflectance], axis=1)
# df.to_csv(os.path.join(odir, 'trios_iwr_' + name + '_idpr' + idpr + '.csv'))
df, wl_ = iwr.reader(lat, lon, alt) # depth data already corrected for sensor position , delta_Lu_depth=0.07,delta_Edz_depth=-0.28)
else:
continue
# ---------------------------
# Data formatting
# ---------------------------
# set pandas dataframes with data and parameters
df['rounded_depth', ''] = df.prof_Edz.round(1)
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)
N = len(wl_)
x = mean.prof_Edz # - 0.56
x = mean.prof_Edz
depth_ = np.linspace(0, x.max(), 200)
res = fit(N)
res_w, res_rw, res_med = fit(N), fit(N), fit(N)
for idx, wl in enumerate(wl_[:-10]):
aw, bbw = iopw.get_iopw(wl)
F0 = irr.get_F0(wl)
# ---------------------------
# Data processing
# ---------------------------
# load process object
process = iwr_process(wl=wl_)
y = mean.Edz.iloc[:, idx]
sigma = std.Edz.iloc[:, idx]
sigma[sigma < noise] = noise
sigma.fillna(np.inf, inplace=True)
res.popt[idx, :], res.pcov[idx, ...] = curve_fit(iwr_process.f_Edz, x, y, [1.1 * aw, 100],
bounds=([aw, 0], [np.inf, F0]))
res_w.popt[idx, :], res_w.pcov[idx, ...] = curve_fit(iwr_process.f_logEdz, x, np.log(1 + y),
[1.1 * aw, 100], bounds=(
[aw, 0], [np.inf, F0])) # , sigma=sigma, absolute_sigma=True
# res_rw.append( curve_fit(iwr_process.f_Edz, x, y, [1.1 * aw, 100], sigma=sigma, absolute_sigma=False, bounds=([aw, 0], [np.inf, F0])))
y = median.Edz.iloc[:, idx]
res_med.popt[idx, :], res_med.pcov[idx, ...] = curve_fit(iwr_process.f_Edz, x, y, [1.1 * aw, 100],
bounds=([aw, 0], [np.inf, F0]))
# process data
res = process.process(mean, std)
res_w = process.process(mean, std, mode='log')
res_med = process.process(median, std)
# ---------------------------
# Plotting section
# ---------------------------
i = 0
mpl.rcParams.update({'font.size': 12})
......@@ -132,13 +129,13 @@ for idpr in idprs: # [-1:]:
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 = iwr_process().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 = iwr_process().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])
# 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 = iwr_process().f_Edz(depth_, *res_med.popt[idx, :])
ax.plot(depth_, Ed_sim, linestyle='--', c='black', label='median')
ax.semilogy()
......
......@@ -104,7 +104,7 @@ class iwr_data:
return full data frame
:param Edf: file path of irradiance data
:param Edzf: file pat of downward in-water irradiance data
:param Edzf: file path of downward in-water irradiance data
:param Luzf: file path of upward in-water radiance data
:param lat: latitude (decimal)
:param lon: longitude (decimal)
......@@ -116,8 +116,6 @@ class iwr_data:
:return:
'''
df = pd.DataFrame()
# ''' read files with pandas format '''
d = data([1, 0])
......@@ -126,9 +124,9 @@ class iwr_data:
Luz, wl_Luz = d.load_csv(self.Luzf)
#mask negative values TODO save number of discarded data
# Ed.mask(Ed<0,inplace=True)
# Edz.mask(Edz<0,inplace=True)
# Luz.mask(Luz<0,inplace=True)
Ed[Ed < 0] = 0 #.mask(Ed<0,inplace=True)
Edz[Edz < 0] = 0 #.mask(Edz<0,inplace=True)
Luz[Luz < 0] = 0 #.mask(Luz<0,inplace=True)
# copy depth data to Ed frame on date index
# Ed.index = Ed.index.droplevel(level=1)
......@@ -148,10 +146,10 @@ class iwr_data:
newLuz = pd.DataFrame(index=Luz.index, columns=pd.MultiIndex.from_tuples(list(zip(['Luz'] * len(wl), wl)),
names=['param', 'wl']), data=intLuz)
print('read merge ok')
# correct depth data for sensor to sensor distance
newLuz.reset_index(level=1, inplace=True)
newLuz.iloc[:, 0] = Luz.iloc[:, 0] + delta_Lu_depth
newLuz.iloc[:, 0] = newLuz.iloc[:, 0] + delta_Lu_depth
# newEd.reset_index(level=1,inplace=True)
newEdz.reset_index(level=1, inplace=True)
......@@ -259,6 +257,12 @@ class swr_data:
return df, wl
class fit:
def __init__(self, N=0, m=2):
self.popt = np.full([N, m], np.nan)
self.pcov = np.full([N, m, m], np.nan)
class data:
def __init__(self, index_idx=[0]):
# first position should be datetime index
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment