import os, sys import pandas as pd import numpy as np import glob import io import matplotlib as mpl import matplotlib.pyplot as plt import plotly import plotly.graph_objs as go import cmocean import scipy.optimize as so from trios.process import * from trios.utils.sunposition import sunpos dir = '/DATA/OBS2CO/data/sabine/data/raw' odir = '/DATA/OBS2CO/data/sabine/data/awr' dirfig = '/DATA/OBS2CO/data/sabine/fig' file = os.path.join(dir, 'Ed_Ld_Lu_all_highroc_cruise_oslofjord.dat') dateparse = lambda x: pd.datetime.strptime(x, '%Y-%m-%d %H:%M:%S') dff = pd.read_csv(file, sep="\t", date_parser=dateparse, index_col=2) awr = awr_process() awr.get_rho_values([30], [40], [135]) groupname = dff.columns.values[1] # 'RT$station_area[n]' def get_param(df, param='Ld', name='Lsky'): ind = df.columns.str.contains(param) df = df.loc[:, ind] wl = df.columns.str.replace(param + '_', '').astype('float') df.columns = pd.MultiIndex.from_product([[name], wl], names=['param', 'wl']) return df for group, df_ in dff.groupby(groupname): print(group) ofile = os.path.join(odir, 'awr_' + group + '.csv') Rrs_df = pd.DataFrame() # ----------------------------- # data formatting # ----------------------------- raw_ = pd.DataFrame() for time, df in df_.groupby(df_.index): azi = [int(df.iloc[:, 0].unique()[0].split('_')[-1])] lat = df.LAT.values.mean() lon = df.LON.values.mean() ws = [df.WIND_SPEED.values.mean() / 2] # TODO check altitude alt = 0 Lsky = get_param(df, 'Ld', 'Lsky') Lt = get_param(df, 'Lu', 'Lt') Ed = get_param(df, 'Ed', 'Ed') wl = Lsky.columns.levels[1].values # compute solar angle (mean between first and last acquisition time sza = sunpos(time, lat, lon, alt)[1] sza = np.repeat(sza, Lt.__len__()) # Rrs, rho = awr.process(wl, Lt, Lsky, Ed, sza, ws=ws, azi=azi) # keep metadata from de and reshape dataframe df.columns = pd.MultiIndex.from_product([df.columns.values, ['']], names=['param', 'wl']) # ['Rrs(awr)'] * Rrs.shape[1] raw = pd.concat([df.iloc[:, 0:10], Ed, Lsky, Lt], axis=1) raw['sza'] = sza.mean() raw['azi'] = azi[0] # raw['rho']=rho.mean() raw_ = pd.concat([raw_, raw], axis=0) def filtering(Lt, Lsky, Ed, **kargs): ''' :param Lt: :param Lsky: :param Ed: :param kargs: :return: ''' ind_Ed, notused = calc.spectra_median_filter(Ed, kargs) ind_sky, notused = calc.spectra_median_filter(Lsky, kargs) ind = ind_Ed & ind_sky return ind def process_wrapper(wl, df, sza, vza=[40], azi=[135], ws=[2], aot=[0.1], method='M99'): print(sza, vza, azi, ws, aot, method) Rrs, rho = process(wl, df.Lt, df.Lsky.values, df.Ed.values, sza, vza, azi, ws, aot, method) Rrs.columns = pd.MultiIndex.from_product([['Rrs'], wl], names=['param', 'wl']) return Rrs, rho def process(wl, Lt, Lsky, Ed, sza, vza=[40], azi=[135], ws=[2], aot=[0.1], method='M99'): ''' :param wl: :param Lt: :param Lsky: :param Ed: :param sza: :param vza: :param azi: :param ws: :param aot: :param method: 'M99, 'M15, 'osoaa' :return: ''' # ----------------------------- # standard data processing # ----------------------------- if method == 'osoaa': rho = awr.get_rho_values(np.median(sza), vza, azi, wl=wl, ws=ws, aot=aot) elif method == 'M99': rho = awr.get_rho_mobley(awr.rhoM1999, [np.median(sza)], [vza], [azi], [ws]) elif method == 'M15': rho = awr.get_rho_mobley(awr.rhoM2015, [np.median(sza)], [vza], [azi], [ws]) else: return print('ERROR: no method for rho factor') Rrs = (Lt - rho * Lsky) / Ed # Rrs.columns = pd.MultiIndex.from_product([['Rrs'], wl], names=['param', 'wl']) return Rrs, rho.mean() def cost_func(x, param, meas, Rrs_bar): sza, vza, azi = param Lt, Lsky, Ed = meas ws = x rho = awr.get_rho_mobley(awr.rhoM1999, [sza], [vza], [azi], [ws]) Rrs = (Lt - rho * Lsky) / Ed return Rrs - Rrs_bar # x_ave = x.mean() # return np.sum(np.abs(x - x_ave)) def process_optimization(wl, Lt, Lsky, Ed, sza, vza=[40], azi=[135], ws=[2], aot=[0.1], method='M99'): plt.figure() # ------------------------------ # initialization of mean/median values # ------------------------------ Rrs, rho = process(wl, Lt, Lsky, Ed, sza, ws=ws, azi=azi) Rrs_bar = Rrs.mean(axis=0) # ----------------------------- # non-linear optimization # ----------------------------- for j in range(10): x_est = [] res = [] Rrs_est = [] rho_est = [] plt.plot(wl, Rrs_bar) for i in range(len(Lt)): geom = [sza[i], vza, azi] meas = [Lt[i], Lsky[i], Ed[i]] x0 = ws res_lsq = so.least_squares(cost_func, x0, args=(geom, meas, Rrs_bar)) res.append(res_lsq) x_est.append(res_lsq.x[0]) Rrs, rho = process(wl, Lt[i], Lsky[i], Ed[i], sza[i], ws=res_lsq.x[0], azi=azi) Rrs_est.append(Rrs) rho_est.append(rho) 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 # ------------------ # filtering # ------------------ ind = 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 # ----------------------------- Rrs, rho = process_wrapper(wl, raw_, raw_.sza, ws=ws, azi=azi) Rrs_15, rho = process_wrapper(wl, raw_, raw_.sza, ws=ws, azi=azi, method='M15') Rrs_h, rho = process_wrapper(wl, raw_, raw_.sza, ws=ws, azi=azi, method='osoaa') Rrs_opt, Rrs_opt_std = process_optimization(wl, Lt, Lsky, Ed, sza, azi=azi) wl = Rrs.T.index.get_level_values(1) date = Rrs.index.get_level_values(0).date[0].__str__() 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) ax.fill_between(x, mean - std, mean + std, alpha=0.35, color=c) plt.rcParams.update({'font.size': 18}) fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(10, 8)) fig.subplots_adjust(left=0.16, right=0.9, hspace=.5, wspace=0.65) add_curve(ax, wl, Rrs[~ind].transpose().mean(axis=1), Rrs[~ind].transpose().std(axis=1), label='QC-excluded', c='red') add_curve(ax, wl, Rrs[ind].transpose().mean(axis=1), Rrs[ind].transpose().std(axis=1), label='M99, QC-passed', c='black') add_curve(ax, wl, Rrs_h[ind].transpose().mean(axis=1), Rrs_h[ind].transpose().std(axis=1), label='H, QC-passed', c='green') add_curve(ax, wl, Rrs_15[ind].transpose().mean(axis=1), Rrs_15[ind].transpose().std(axis=1), label='M15, QC-passed', c='violet') add_curve(ax, wl, Rrs_opt, Rrs_opt_std, label='Optimization', c='blue') ax.legend(loc='best', frameon=False) ax.set_ylabel(r'$R_{rs}\ (sr^{-1})$') ax.set_xlabel(r'$Wavelength\ (nm)$') fig.savefig(os.path.join(dirfig, 'trios_awr_' + date + '_' + group + '.png')) plt.close() raw_df = pd.concat([raw_df, raw_], axis=0) Rrs_df = pd.concat([Rrs_df, Rrs], axis=1) Rrs_df.to_csv(ofile) if False: dff = Rrs # _df # dff.columns = pd.MultiIndex.from_product([['Rrs(awr)'], dff.columns], names=['param', 'wl']) # wl = dff.columns.droplevel().values trace = [] i = 0 for station, df in dff.groupby('RT$station_area[n]'): i = i + 1 for time, df_ in df.iterrows(): print(time) trace.append(go.Scattergl( x=wl, # spectrum, y=df_.loc['Rrs(awr)'].values, text='', name=station, legendgroup=station, mode='lines', line={ 'color': 'rgba' + str(cmocean.cm.thermal(i * 25)[:-1])} # 'line': go.scatter.Line(color='rgba({}, {}, {}, {})'.format(*s_m.to_rgba(parameters[i]).flatten()), width=2), , )) plotly.offline.plot(trace)