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

inclusion of iwr process in main, fix for first guess x0 in iwr.process

parent 9c5a17f4
......@@ -75,10 +75,26 @@ Usage:
## Running the tests
- For `awr` data:
```
trios_processing ./test/data/ 150 awr --lat 42.30351823 --lon 9.462897398 --odir ./test/results --no_clobber
trios_processing ./test/data/ 150 awr --lat 42.30351823 --lon 9.462897398 --odir ./test/results --method M99 --name _M99 --plot --figdir ./test/fig
trios_processing ./test/data/ 150 awr --lat 42.30351823 --lon 9.462897398 --odir ./test/results --method osoaa --name _osoaa --plot
trios_processing ./test/data/ 150 awr --lat 42.30351823 --lon 9.462897398 --odir ./test/results --method osoaa --name _osoaa --plot --figdir ./test/fig
trios_processing ./test/data/ 150 awr --lat 42.30351823 --lon 9.462897398 --odir ./test/results --method osoaa --ws 5 --aot=0.2 --name _osoaaws5_aot0.1 --plot --figdir ./test/fig
trios_processing ./test/data/ 150 awr --lat 42.30351823 --lon 9.462897398 --odir ./test/results --method temp_opt --name _optimization --plot --figdir ./test/fig
```
- For `iwr` data:
```
trios_processing ./test/data/ 150 iwr --lat 42.30351823 --lon 9.462897398 --odir ./test/results --no_clobber --plot --figdir ./test/fig
```
- For `swr` data:
```
trios_processing ./test/data/ 150 swr --lat 42.30351823 --lon 9.462897398 --odir ./test/results --no_clobber --plot --figdir ./test/fig
```
......
This diff is collapsed.
......@@ -16,3 +16,8 @@ F0_file = os.path.join(root, '../aux/Thuillier_2003_0.3nm.dat')
wl_step = 3
wl_common = np.arange(320,950+wl_step,wl_step)
# wavelength indexes for iwr plotting
# for close red/IR
# idx_list_for_plot=(165,170,175,180,185,190,195,200,205)
# through visible/IR
idx_list_for_plot=(28, 37, 51, 71, 91, 105, 130, 140, 170)
\ No newline at end of file
......@@ -3,7 +3,8 @@
Usage:
trios_processing <input_dir> <IDpr> <measurement_type> --lat <lat> --lon <lon> \
[--altitude=alt] [--ofile <ofile>] [--odir <odir>] [--plot] [--figdir <figdir>] \
[--name <name>] [--method <method>] [--no_clobber]
[--name <name>] [--method <method>] [--no_clobber] \
[--vza <vza>] [--azi <azi>] [--ws <ws>] [--aot <aot>]
trios_processing -h | --help
trios_processing -v | --version
......@@ -26,6 +27,14 @@ Options:
--name name Keyword to append to file name and figures [default: ]
--method method Keyword for the method to apply for data processing.
For awr: M99, M15, osoaa, temp_opt [default: M99]
--vza vza Viewing zenith angle (in deg) of the awr radiometer (applicable for all awr methods)
[default: 40]
--azi azi Relative azimuth angle (in deg) of the awr radiometer (applicable for all awr methods)
[default: 135]
--ws ws Wind speed (in m/s) (applicable for awr methods: M99, M15, osoaa)
[default: 2]
--aot aot Aerosol optical thickness at 550 nm (applicable for awr methods: osoaa)
[default: 0.1]
--no_clobber Do not process <input_dir> <IDpr> files if <output_file> already exists.
'''
......@@ -65,6 +74,10 @@ def main():
plot = args['--plot']
figdir = os.path.abspath(args['--figdir'])
noclobber = args['--no_clobber']
azi = float(args['--azi'])
vza = float(args['--vza'])
ws = float(args['--ws'])
aot = float(args['--aot'])
try:
type_ = type_list[meas_type]
......@@ -103,17 +116,15 @@ def main():
ax.set_ylabel(r'$R_{rs}\ (sr^{-1})$')
ax.set_xlabel(r'Wavelength (nm)')
ax.set_title('ID: ' + idpr + ', ' + date + ', sza=' + str(round(df.sza.mean(), 2)))
fig.savefig(os.path.join(figdir, 'trios_swr_' + date + '_idpr' + idpr + name+'.png'), bbox_inches='tight')
fig.savefig(os.path.join(figdir, 'trios_swr_' + date + '_idpr' + idpr + name + '.png'),
bbox_inches='tight')
plt.close()
elif meas_type == 'awr':
# -----------------------------------------------
# AWR processing
# -----------------------------------------------
azi = 135
vza = 40
ws = 2
method = 'M99'
uawr = u.awr_data(idpr, files)
if uawr.Edf:
df, wl = uawr.reader(lat, lon, alt)
......@@ -127,17 +138,39 @@ def main():
print('Skip processing: data already processed with "--no_clobber" set')
return
awr = awr_process(df, wl, name, idpr)
if plot:
figfile = os.path.join(figdir, 'trios_awr_' + date + '_idpr' + idpr + name +'.png')
figfile = os.path.join(figdir, 'trios_awr_' + date + '_idpr' + idpr + name + '.png')
else:
figfile = ""
Rrs = awr.call_process(method, ofile, vza=vza, azi=azi, plot_file=figfile)
awr = awr_process(df, wl, name, idpr)
Rrs = awr.call_process(method, ofile, vza=vza, azi=azi, ws=ws, aot=aot, plot_file=figfile)
elif meas_type == 'iwr':
pass
# -----------------------------------------------
# IWR processing
# -----------------------------------------------
uiwr = u.iwr_data(idpr, files)
if uiwr.file:
df, wl = uiwr.reader(lat, lon, alt)
date = df.index[0].date().__str__()
if ofile:
ofile = os.path.join(odir, ofile)
else:
ofile = os.path.join(odir, 'Rrs_iwr_' + date + '_idpr' + idpr + name + '.csv')
if noclobber and os.path.exists(ofile):
print('Skip processing: data already processed with "--no_clobber" set')
return
if plot:
figfile = os.path.join(figdir, 'trios_iwr_' + date + '_idpr' + idpr + name + '.pdf')
else:
figfile = ""
iwr = iwr_process(df, wl, name, idpr)
results = iwr.call_process(ofile, plot_file=figfile)
if __name__ == "__main__":
......
# standard lib
import pandas as pd
from scipy import interpolate, integrate
import scipy.optimize as so
# for plotting pruposes
import cmocean
import plotly.graph_objs as go
from matplotlib.backends.backend_pdf import PdfPages
import matplotlib as mpl
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 18})
# package
from trios.utils.utils import plot as up
import trios.utils.utils as uu
from trios.utils.utils import reshape as r
......@@ -92,7 +96,7 @@ class awr_process:
def call_process(self, method='M99', ofile="", vza=40, azi=135, ws=2, aot=0.1, plot_file=""):
wl = self.wl
vza, azi, ws = [vza], [azi], [ws] # formatting for interpolation functions
vza, azi, ws, aot = [vza], [azi], [ws], [aot] # formatting for interpolation functions
# ------------------
# filtering
......@@ -111,7 +115,7 @@ class awr_process:
elif method == 'osoaa':
Rrs, rho = self.process_wrapper(wl, clean, clean.sza, vza=vza, azi=azi, ws=ws, aot=aot, method=method)
elif method == 'temp_opt':
Rrs, Rrs_opt_std = self.process_optimization(wl, Lt, Lsky, Ed, sza, vza=vza, azi=azi)
Rrs, rho, Rrs_opt, Rrs_opt_std = self.process_optimization(wl, Lt, Lsky, Ed, sza, vza=vza, azi=azi)
self.Rrs = Rrs
......@@ -121,7 +125,6 @@ class awr_process:
Rrs_stat = Rrs_stat.T
Rrs_stat.to_csv(ofile)
if plot_file:
# ------------------
# plotting
......@@ -155,21 +158,21 @@ class awr_process:
up.add_curve(ax, wl, Lt.mean(axis=0), Lt.std(axis=0),
label=r'$L_t$', c='black')
up.add_curve(ax, wl, Lsky.mean(axis=0) * rho, Lsky.std(axis=0) * rho,
label=method+' (' + str(round(rho, 4)) + ')', c='violet')
label=method + ' (' + str(round(rho, 4)) + ')', c='violet')
ax.set_ylabel(r'$L_t\ or L_{surf}$')
ax.set_xlabel(r'Wavelength (nm)')
# ---- Proportion o(Lt - Lsurf ) /Lt
ax = axs[0, 2]
up.add_curve(ax, wl, Lsky.mean(axis=0) * rho / Ltm, Lsky.std(axis=0) * rho,
label=method+' (' + str(round(rho, 4)) + ')', c='violet')
label=method + ' (' + str(round(rho, 4)) + ')', c='violet')
ax.set_ylabel(r'$L_{surf}/L_t$')
ax.set_xlabel(r'Wavelength (nm)')
# ---- Lw
ax = axs[1, 0]
up.add_curve(ax, wl, Rrs.mean(axis=0) * Edm, Rrs.std(axis=0) * Edm,
label=method+' (' + str(round(rho, 4)) + ')', c='violet')
label=method + ' (' + str(round(rho, 4)) + ')', c='violet')
ax.set_ylabel(r'$L_{w}\ (sr^{-1})$')
ax.set_xlabel(r'Wavelength (nm)')
......@@ -177,7 +180,7 @@ class awr_process:
# ---- Rrs
ax = axs[1, 1]
up.add_curve(ax, wl, Rrs.transpose().mean(axis=1), Rrs.transpose().std(axis=1),
label=method+' (' + str(round(rho, 4)) + ')', c='violet')
label=method + ' (' + str(round(rho, 4)) + ')', c='violet')
ax.set_ylabel(r'$R_{rs}\ (sr^{-1})$')
ax.set_xlabel(r'Wavelength (nm)')
......@@ -263,7 +266,7 @@ class awr_process:
# initialization of mean/median values
# ------------------------------
Rrs, rho = self.process(wl, Lt, Lsky, Ed, sza, ws=ws, azi=azi)
Rrs_bar = Rrs.mean(axis=0)
Rrs_bar = np.nanmedian(Rrs, axis=0)
# -----------------------------
# non-linear optimization
......@@ -288,7 +291,9 @@ class awr_process:
print(res_lsq.x, res_lsq.cost)
Rrs_bar = np.mean(Rrs_est, axis=0)
Rrs_std = np.std(Rrs_est, axis=0)
return Rrs_bar, Rrs_std
Rrs_est = pd.DataFrame(Rrs_est)
Rrs_est.columns = pd.MultiIndex.from_product([['Rrs'], wl], names=['param', 'wl'])
return Rrs_est, np.mean(rho_est, axis=0), Rrs_bar, Rrs_std
@classmethod
def filtering(cls, Lt, Lsky, Ed, **kargs):
......@@ -400,9 +405,256 @@ class swr_process:
class iwr_process:
def __init__(self, df=None, wl=None, ):
def __init__(self, df=None, wl=None, name="", idpr=""):
self.df = df
self.wl = wl
self.name = name
self.idpr = idpr
def call_process(self, ofile="", plot_file=""):
# ---------------------------
# Data formatting
# ---------------------------
df = self.df
wl = self.wl
# set pandas dataframes with data and parameters
df['rounded_depth', ''] = df.prof_Edz.round(1)
# data filtering
df[df.rounded_depth < -0.05] = np.nan
# df[df == 0] = np.nan
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)
# ---------------------------
# 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')
Kd = res_lsq.popt[:, 0]
Edm = res_lsq.popt[:, 1]
Es = mean.Ed.mean()
Lwm = res_lsq.popt[:, 3]
rrs = Lwm / Edm
Es_from_Edm = 0.96 * Edm
Lwp = 0.541 * Lwm
Rrs_Edp = Lwp / Es
Rrs_Edm = Lwp / Es_from_Edm
Rrs_df = pd.DataFrame(Rrs_Edp)
Rrs_df.columns = ['Rrs']
Rrs_df['Kd'] = Kd
if ofile:
# -----------------------------------------------
# write result data
# -----------------------------------------------
Rrs_df.to_csv(ofile)
if plot_file:
# ---------------------------
# Plotting section
# ---------------------------
pdf = PdfPages(plot_file)
# -------------------------------------
# Profile plotting
# -------------------------------------
# Ed profile
i = 0
x = mean.prof_Edz
depth_ = np.linspace(0, x.max(), 200)
mpl.rcParams.update({'font.size': 12})
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)
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,
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, :])
ax.plot(depth_, Ed_sim, linestyle='-', c='black', label='mean')
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, :])
ax.plot(depth_, Ed_sim, linestyle='--', c='black', label='median')
ax.semilogy()
# ax.colorbar()
ax.set_ylabel(r'$E_d\ ({mW\cdot m^{-2}\cdot nm^{-1}})$')
ax.set_xlabel('Depth (m)')
ax.set_title(r'$\lambda = $' + str(round(wl[idx], 1)) + ' nm, Kd = ' +
str(round(res_w.popt[idx, 0], 3)) + '$m^{-1}$, Ed0 =' + str(round(res_w.popt[idx, 1], 1)),
fontsize=12)
ax.legend(loc='best', frameon=False)
fig.suptitle('trios_iwr ' + self.name + ' idpr' + self.idpr, fontsize=16)
# fig.savefig(os.path.join(dirfig,'trios_iw_idpr' + idpr + '_'+name+'.png'),dpi=600)
# fig.savefig(os.path.join(dirfig, 'trios_iwr_Ed_profile_idpr' + idpr + '_' + name + '.pdf'))
pdf.savefig()
plt.close()
# Lu profile
i = 0
x = mean.prof_Luz
depth_ = np.linspace(0, x.max(), 200)
mpl.rcParams.update({'font.size': 12})
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)
for idx in idx_list_for_plot:
ax = axs.flat[i]
i += 1
y = mean.Luz.iloc[:, idx]
sigma = std.Luz.iloc[:, idx]
yerr = [y - q25.Luz.iloc[:, idx], q75.Luz.iloc[:, idx] - y]
yerr = sigma
ax.errorbar(x, y, yerr=yerr, ecolor='lightgray', alpha=0.43, fmt='none', elinewidth=3, capsize=0,
label='std')
ax.scatter(x, y,
# s = std.Luz.iloc[:,50]/std.Luz.iloc[:,50].mean()+100,
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:])
ax.plot(depth_, Lu_sim, linestyle='-', c='black', label='mean')
# Lu_sim = iwr_process().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])
# # ax.plot(depth_, Lu_sim, linestyle=':', c='red', label='mean, relativ weighted')
# Lu_sim = iwr_process().f_Luz(depth_, *res_med.popt[idx, :])
# ax.plot(depth_, Lu_sim, linestyle='--', c='black', label='median')
ax.semilogy()
# ax.colorbar()
ax.set_ylabel(r'$L_u\ ({mW\cdot m^{-2}\cdot sr^{-1}\cdot nm^{-1}})$')
ax.set_xlabel('Depth (m)')
ax.set_title(r'$\lambda = $' + str(round(wl[idx], 1)) + r'$ nm, K_{Lu} = $' +
str(round(res_lsq.popt[idx, 2], 3)) + '$m^{-1}$, Lu0 =' + str(
round(res_lsq.popt[idx, 3], 1)),
fontsize=12)
ax.legend(loc='best', frameon=False)
fig.suptitle('trios_iwr ' + self.name + ' idpr' + self.idpr, fontsize=16)
# fig.savefig(os.path.join(dirfig,'trios_iw_idpr' + idpr + '_'+name+'.png'),dpi=600)
# fig.savefig(os.path.join(dirfig, 'trios_iwr_Lu_profile_idpr' + idpr + '_' + name + '.pdf'))
pdf.savefig()
plt.close()
# -------------------------------------
# Spectra plotting
# -------------------------------------
mpl.rcParams.update({'font.size': 18})
fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(15, 12))
fig.subplots_adjust(left=0.1, right=0.9, hspace=.5, wspace=0.29)
iparam = 0
ax = axs.flat[iparam]
y, sigma = res.popt[:, iparam], res.pcov[:, iparam, iparam]
up.add_curve(ax, wl, y, sigma, c='red', label='mean')
y, sigma = res_w.popt[:, iparam], res_w.pcov[:, iparam, iparam]
up.add_curve(ax, wl, y, sigma, c='orange', label='log-space')
y, sigma = res_med.popt[:, iparam], res_med.pcov[:, iparam, iparam]
up.add_curve(ax, wl, y, sigma, c='green', label='median')
y, sigma = res_lsq.popt[:, iparam], np.sqrt(res_lsq.pcov[:, iparam, iparam])
up.add_curve(ax, wl, y, sigma, c='blue', label='lsq')
ax.set_ylabel(r'$K_{d}\ ({m^{-1}})$')
ax.set_xlabel(r'Wavelength (nm)')
iparam = 1
ax = axs.flat[iparam]
y, sigma = res.popt[:, iparam], np.sqrt(res.pcov[:, iparam, iparam])
up.add_curve(ax, wl, y, sigma, c='red', label='mean')
y, sigma = res_w.popt[:, iparam], np.sqrt(res_w.pcov[:, iparam, iparam])
up.add_curve(ax, wl, y, sigma, c='orange', label='log-space')
y, sigma = res_med.popt[:, iparam], np.sqrt(res_med.pcov[:, iparam, iparam])
up.add_curve(ax, wl, y, sigma, c='green', label='median')
y, sigma = res_lsq.popt[:, iparam], np.sqrt(res_lsq.pcov[:, iparam, iparam])
up.add_curve(ax, wl, y, sigma, c='blue', label='lsq')
up.add_curve(ax, wl, mean.Ed.mean(), mean.Ed.std(), c='black', label='Es')
ax.legend(loc='best', frameon=False)
ax.set_ylabel(r'$E_{d}(0^{-})\ ({mW\cdot m^{-2}\cdot nm^{-1}})$')
ax.set_xlabel(r'Wavelength (nm)')
iparam = 2
ax = axs.flat[iparam]
y, sigma = res_lsq.popt[:, iparam], np.sqrt(res_lsq.pcov[:, iparam, iparam]) * 0
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}})$')
ax.set_xlabel(r'Wavelength (nm)')
iparam = 3
ax = axs.flat[iparam]
y, sigma = res_lsq.popt[:, iparam], np.sqrt(res_lsq.pcov[:, iparam, iparam]) * 0
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})$')
ax.set_xlabel(r'Wavelength (nm)')
fig.suptitle('trios_iwr ' + self.name + ' idpr' + self.idpr, fontsize=16)
# fig.savefig(os.path.join(dirfig, 'trios_iwr_l2_idpr' + idpr + '_' + name + '.pdf'))
pdf.savefig()
plt.close()
# Edm = res_lsq.popt[:,1]
# Es = mean.Ed.mean()
# Lwm = res_lsq.popt[:,3]
# rrs = Lwm / Edm
# Es_from_Edm = 0.96 * Edm
# Lwp = 0.541 * Lwm
# Rrs_Edp = Lwp / Es
# Rrs_Edm = Lwp / Es_from_Edm
mpl.rcParams.update({'font.size': 18})
fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(15, 12))
fig.subplots_adjust(left=0.1, right=0.9, hspace=.5, wspace=0.29)
ax = axs.flat[0]
up.add_curve(ax, wl, Es / Edm, c='red', label='Es over Ed0-')
ax.set_xlabel(r'Wavelength (nm)')
ax.legend(loc='best', frameon=False)
ax = axs.flat[1]
up.add_curve(ax, wl, rrs, c='red', label='rrs')
ax.set_xlabel(r'Wavelength (nm)')
ax.legend(loc='best', frameon=False)
ax = axs.flat[2]
up.add_curve(ax, wl, Rrs_Edm, c='red', label='Rrs (from Ed0-)')
up.add_curve(ax, wl, Rrs_Edp, c='blue', label='Rrs (from Ed0+)')
ax.set_xlabel(r'Wavelength (nm)')
ax.legend(loc='best', frameon=False)
ax = axs.flat[3]
fig.suptitle('trios_iwr ' + self.name + ' idpr' + self.idpr, fontsize=16)
pdf.savefig()
plt.close()
pdf.close()
def process(self, meas, std, mode='linear'):
wl_ = self.wl
......@@ -446,8 +698,8 @@ class iwr_process:
sigma = (sig_Edz, sig_Luz)
sigma = (1, 1)
x0 = [1.1 * aw, meas.Ed.iloc[:, idx].mean(), 1.1 * aw, meas.Luz.iloc[0, idx]]
x0 = [1.1 * aw, min(F0,meas.Ed.iloc[:, idx].mean()), 1.1 * aw, meas.Luz.iloc[0, idx]]
print('x0',wl,F0,aw, x0)
lsq = so.least_squares(self.cost_func, x0, args=(z, y, sigma),
bounds=([aw, 0, aw / 2, 0], [np.inf, F0, np.inf, np.inf]))
cost = 2 * lsq.cost # res.cost is half sum of squares!
......
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