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

test several optimization methods

parent 5f63794c
<?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>
......
<?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
......@@ -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" />
......
......@@ -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)
......
......@@ -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]
......
......@@ -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
......
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