diff --git a/.idea/misc.xml b/.idea/misc.xml index f8f8c44a92f31e0ec4149ad5fb04a57c4e1bf124..a32de7b6f3ee623ad5c65875b13b18844631a8d7 100644 --- a/.idea/misc.xml +++ b/.idea/misc.xml @@ -1,6 +1,6 @@ <?xml version="1.0" encoding="UTF-8"?> <project version="4"> - <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.6" project-jdk-type="Python SDK" /> + <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.5" project-jdk-type="Python SDK" /> <component name="PyCharmProfessionalAdvertiser"> <option name="shown" value="true" /> </component> diff --git a/.idea/vcs.xml b/.idea/vcs.xml index 830674470f8052bb64c64a5b513f5df2815a70c1..94a25f7f4cb416c083d265558da75d457237d671 100644 --- a/.idea/vcs.xml +++ b/.idea/vcs.xml @@ -1,7 +1,6 @@ <?xml version="1.0" encoding="UTF-8"?> <project version="4"> <component name="VcsDirectoryMappings"> - <mapping directory="" vcs="Git" /> <mapping directory="$PROJECT_DIR$" vcs="Git" /> </component> </project> \ No newline at end of file diff --git a/.idea/visu_trios.iml b/.idea/visu_trios.iml index f3cf5399aa748f54217941acbf33f67672ac22f2..64028cfc53dd142c16f1ffc6572373bdfc6826b5 100644 --- a/.idea/visu_trios.iml +++ b/.idea/visu_trios.iml @@ -2,7 +2,7 @@ <module type="PYTHON_MODULE" version="4"> <component name="NewModuleRootManager"> <content url="file://$MODULE_DIR$" /> - <orderEntry type="jdk" jdkName="Python 3.6" jdkType="Python SDK" /> + <orderEntry type="jdk" jdkName="Python 3.5" jdkType="Python SDK" /> <orderEntry type="sourceFolder" forTests="false" /> <orderEntry type="library" name="R User Library" level="project" /> <orderEntry type="library" name="R Skeletons" level="application" /> diff --git a/process/process.py b/process/process.py index 4449d11613d2280006d1aca5a5a9d47277e146cb..3c25eb7f467ca22cfd64fa61f07d7a40a9e8fa8b 100644 --- a/process/process.py +++ b/process/process.py @@ -219,6 +219,11 @@ class iwr_process: '''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) diff --git a/process/run_iwr.py b/process/run_iwr.py index d0dfd488c263b9465582d4016f4ee50599c98fb0..e5d120e7bc4a5a042429b6e1dc28a6a91d34c5d9 100644 --- a/process/run_iwr.py +++ b/process/run_iwr.py @@ -5,6 +5,10 @@ import glob import re import matplotlib as mpl +import matplotlib.pyplot as plt +import cmocean + + import plotly import plotly.graph_objs as go @@ -13,6 +17,14 @@ import utils.auxdata as ua from process.process import * + + +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) + # import aeronet # from config import * @@ -36,7 +48,8 @@ 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 idpr = '178' # loop over idpr for idpr in idprs:#[-1:]: @@ -68,15 +81,6 @@ for idpr in idprs:#[-1:]: q25 = df.groupby('rounded_depth').quantile(0.25) q75 = df.groupby('rounded_depth').quantile(0.75) - import matplotlib.pyplot as plt - import cmocean - - - 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) N = len(wl_) x = mean.prof_Edz #- 0.56 @@ -90,16 +94,18 @@ for idpr in idprs:#[-1:]: 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_Edz, x, y, [1.1 * aw, 100], sigma=sigma, absolute_sigma=True, 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])) - i=0 + + 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 (28, 37, 51, 71, 91, 105, 130, 140, 170): @@ -114,21 +120,21 @@ 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[idx][0]) + 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[idx][0]) - ax.plot(depth_, Ed_sim, linestyle=':', c='black', label='mean, weighted') + 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[idx][0]) + Ed_sim = iwr_process.f_Edz(depth_, *res_med.popt[idx,:]) ax.plot(depth_, Ed_sim, linestyle='--', c='black', label='median') - # ax.semilogy() + 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[idx][0][0], 3)) + '$m^{-1}$, Ed0 =' + str(round(res_w[idx][0][1], 1))) + 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 ' + name + ' idpr' + idpr, fontsize=16) # fig.savefig(os.path.join(dirfig,'trios_iw_idpr' + idpr + '_'+name+'.png'),dpi=600) @@ -142,28 +148,38 @@ for idpr in idprs:#[-1:]: ax.fill_between(x, mean - std, mean + std, alpha=0.35, color=c) - fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(12, 7)) + + mpl.rcParams.update({'font.size': 18}) + fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(15, 6)) + fig.subplots_adjust(left=0.1, right=0.9, hspace=.5, wspace=0.29) iparam = 0 ax = axs.flat[iparam] y, std = res.popt[:,iparam], res.pcov[:,iparam,iparam] add_curve(ax, wl_, y, std, c='red', label='mean') y, std = res_w.popt[:,iparam], res_w.pcov[:,iparam,iparam] - add_curve(ax, wl_, y, std, c='black', label='mean, weighted') + add_curve(ax, wl_, y, std, c='orange', label='log-space') y, std = res_med.popt[:,iparam], res_med.pcov[:,iparam,iparam] - add_curve(ax, wl_, y, std, c='black', label='mean, weighted') + add_curve(ax, wl_, y, std, c='green', label='median') + ax.set_ylabel(r'$K_{d}\ ({m^{-1}})$') + ax.set_xlabel(r'Wavelength (nm)') iparam = 1 ax = axs.flat[iparam] - y, std = res.popt[:,iparam], res.pcov[:,iparam,iparam] + y, std = res.popt[:,iparam], np.sqrt(res.pcov[:,iparam,iparam]) add_curve(ax, wl_, y, std, c='red', label='mean') - y, std = res_w.popt[:,iparam], res_w.pcov[:,iparam,iparam] - add_curve(ax, wl_, y, std, c='black', label='mean, weighted') + y, std = res_w.popt[:,iparam], np.sqrt(res_w.pcov[:,iparam,iparam]) + add_curve(ax, wl_, y, std, c='orange', label='log-space') + y, std = res_med.popt[:,iparam], np.sqrt(res_med.pcov[:,iparam,iparam]) + add_curve(ax, wl_, y, std, c='green', label='median') + add_curve(ax, wl_, mean.Ed.mean(), mean.Ed.std(), c='black', label='Es') - plt.plot(wl_,res.popt[:,iparam],) - plt.fill_between(wl_, - res.popt[:,iparam]- res.pcov[:,iparam,iparam], - res.popt[:,iparam]+ res.pcov[:,iparam,iparam],alpha=0.35, color='grey') + 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)') + fig.suptitle('trios_iwr ' + name + ' idpr' + idpr, fontsize=16) + fig.savefig(os.path.join(dirfig, 'trios_iwr_l2_idpr' + idpr + '_'+name+ '.pdf')) + plt.close() # fig, ax = plt.subplots() # N=df.Edz.shape[1] diff --git a/utils/utils.py b/utils/utils.py index 5ba29766260f2bf24f1b4c6c07d8abcc3828d837..53ecc6918784092d8b477577764b0d22d29b68a7 100644 --- a/utils/utils.py +++ b/utils/utils.py @@ -133,13 +133,13 @@ class iwr_data: #''' interpolate Ed and Lsky data upon Lt wavelength''' wl = wl_Luz - Luz.columns = pd.MultiIndex.from_tuples(zip(['Luz'] * len(wl), wl), names=['param', 'wl']) + Luz.columns = pd.MultiIndex.from_tuples(list(zip(['Luz'] * len(wl), wl)), names=['param', 'wl']) intEd = interp1d(wl_Ed, Ed.values, fill_value='extrapolate')(wl) newEd = pd.DataFrame(index=Ed.index.get_level_values(0), - columns=pd.MultiIndex.from_tuples(zip(['Ed'] * len(wl), wl), names=['param', 'wl']), + columns=pd.MultiIndex.from_tuples(list(zip(['Ed'] * len(wl), wl)), names=['param', 'wl']), data=intEd) intEdz = interp1d(wl_Edz, Edz.values, fill_value='extrapolate')(wl) - newEdz = pd.DataFrame(index=Edz.index, columns=pd.MultiIndex.from_tuples(zip(['Edz'] * len(wl), wl), + newEdz = pd.DataFrame(index=Edz.index, columns=pd.MultiIndex.from_tuples(list(zip(['Edz'] * len(wl), wl)), names=['param', 'wl']), data=intEdz) # correct depth data for sensor to sensor distance