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

implementation of in-water radiometry (optimization scheme Edz)

parent bfcbf70f
# Default ignored files
/workspace.xml
\ No newline at end of file
import sys, os
import re
import datetime as dt
import julian
dublin_julian_day = 2415020
djd = 42349.44914
jd = djd + dublin_julian_day
dt = julian.from_jd(jd)
print(dt)
\ No newline at end of file
import glob
import re
import matplotlib as mpl
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 18})
from scipy.interpolate import interp1d
import pyodbc
import fiona
import pandas as pd
import pandas_access as pda
import geopandas as gpd
fiona.drvsupport.supported_drivers['kml'] = 'rw' # enable KML support which is disabled by default
fiona.drvsupport.supported_drivers['KML'] = 'rw' # enable KML support which is disabled by default
from trios.utils.sunposition import sunpos
from trios.utils import utils as u
from trios.utils.utils import plot as up
from trios.process import *
project_folder = '/DATA/projet/petit-saut/'
dirfig = os.path.join(project_folder, 'fig')
odir = os.path.join(project_folder, 'data/L2')
coordf = os.path.join(project_folder, 'data/gps_points.kml')
coords = gpd.read_file(coordf)
awrfiles = glob.glob(os.path.join(project_folder, 'data/csv/aw*.csv'))
swrfiles = glob.glob(os.path.join(project_folder, 'data/csv/Lu0*.csv'))
# get idpr
idprs = np.unique([re.split('\.', re.split(r'idpr', x)[1])[0] for x in awrfiles])
idprs = np.unique(np.append(idprs, [re.split('\.', re.split(r'idpr', x)[1])[0] for x in swrfiles]))
idpr = idprs[0]
for idpr in idprs:
name = idpr
for idx, p in enumerate(coords.Name):
if p in idpr:
break
c = coords.iloc[idx]
print(c)
lat = c.geometry.y
lon = c.geometry.x
alt = 35 # c.geometry.z.values[0]
# -----------------------------------------------
# SWR processing
# -----------------------------------------------
uswr = u.swr_data(idpr, swrfiles)
if uswr.file:
df, wl_swr = uswr.reader(lat, lon, alt)
df['sza', ''] = np.nan
for index, row in df.iterrows():
# print index
sza = sunpos(index, lat, lon, alt)[1]
df.at[index, 'sza'] = sza
swr = swr_process(df, wl_swr)
Rrs_swr = swr.call_process()
date = index.date().__str__()
mpl.rcParams.update({'font.size': 18})
fig, ax = plt.subplots( figsize=(7, 6))
up.add_curve(ax, wl_swr, Rrs_swr.transpose().mean(axis=1), Rrs_swr.transpose().std(axis=1), label='swr',
c='black')
ax.set_ylabel(r'$R_{rs}\ (sr^{-1})$')
ax.set_xlabel(r'Wavelength (nm)')
ax.set_title('ID: '+idpr+', '+date+', sza=' + str(round(sza.mean(), 2)))
fig.savefig(os.path.join(dirfig, 'trios_swr_' + date + '_idpr' + idpr + '.png'), bbox_inches='tight')
plt.close()
ofile = os.path.join(odir, 'Rrs_swr_' + date + '_idpr' + idpr + '_PSA_guyane.csv')
Rrs_stat = Rrs_swr.describe()
Rrs_stat.columns=Rrs_stat.columns.droplevel()
Rrs_stat = Rrs_stat.T
Rrs_stat.to_csv(ofile,mode='a')
#
# -----------------------------------------------
# AWR processing
# -----------------------------------------------
azi = 135
vza = 40
awr = u.awr_data(idpr, awrfiles)
if awr.Edf:
index_idx = [0]
d = u.data(index_idx)
Ed, wl_Ed = d.load_csv(awr.Edf)
Lsky, wl_Lsky = d.load_csv(awr.Lskyf)
Lt, wl_Lt = d.load_csv(awr.Ltf)
# ''' interpolate Ed and Lsky data upon Lt wavelength'''
wl = wl_Lt
Lt.columns = pd.MultiIndex.from_tuples(zip(['Lt'] * len(wl), wl), names=['param', 'wl'])
intEd = interp1d(wl_Ed, Ed.values, fill_value='extrapolate')(wl)
newEd = pd.DataFrame(index=Ed.index,
columns=pd.MultiIndex.from_tuples(zip(['Ed'] * len(wl), wl), names=['param', 'wl']),
data=intEd)
intLsky = interp1d(wl_Lsky, Lsky.values, fill_value='extrapolate')(wl)
newLsky = pd.DataFrame(index=Lsky.index, columns=pd.MultiIndex.from_tuples(zip(['Lsky'] * len(wl), wl),
names=['param', 'wl']), data=intLsky)
awr = awr_process()
ws = [2]
print(azi, vza)
Lsky = newLsky # .loc[(newLsky.index.get_level_values(1) == vza) & (newLsky.index.get_level_values(2) == azi)]
Ed = newEd # .loc[(newEd.index.get_level_values(1) == vza) & (newEd.index.get_level_values(2) == azi)]
# Lsky_idx = Lsky.index
# Ed_idx= Ed.index
# Lt_idx = Lt.index
# Lsky.reset_index(level=[1,2],inplace=True)
# Ed.reset_index(level=[1,2],inplace=True)
# Lt.reset_index(level=[1,2],inplace=True)
# merge sensor data on time
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
raw['sza', ''] = np.nan
for index, row in raw.iterrows():
# print index
sza = sunpos(index, lat, lon, alt)[1]
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
# ------------------
Ltm = Lt.mean(axis=0)
Edm = Ed.mean(axis=0)
mpl.rcParams.update({'font.size': 18})
fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(20, 12))
fig.subplots_adjust(left=0.1, right=0.9, hspace=.5, wspace=0.45)
# ---- Ed
ax = axs[0, 0]
up.add_curve(ax, wl, Ed.mean(axis=0),
label=r'$L_{sky}$', c='red') # just to put the two labels
up.add_curve(ax, wl, Ed.mean(axis=0), Ed.std(axis=0),
label=r'$E_s$', c='black')
ax.set_ylabel(r'$E_{d}(0^{+})$')
# ---- Lsky
ax2 = ax.twinx()
up.add_curve(ax2, wl, Lsky.mean(axis=0), Lsky.std(axis=0),
label=r'$L_{sky}$', c='red')
ax2.set_ylabel(r'$L_{sky}$', color='r')
ax2.tick_params('y', colors='r')
ax.set_xlabel(r'Wavelength (nm)')
ax.legend(loc='best', frameon=False)
# ---- Lt vs Lsurf
ax = axs[0, 1]
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) * rho15, Lsky.std(axis=0) * rho15,
label='M2015 (' + str(round(rho15, 4)) + ')', c='violet')
up.add_curve(ax, wl, Lsky.mean(axis=0) * rho99, Lsky.std(axis=0) * rho99, c='orange',
label='M1999(' + str(round(rho99, 4)) + ')')
up.add_curve(ax, wl, Lsky.mean(axis=0) * rho_h, Lsky.std(axis=0) * rho_h, c='green',
label='h(' + str(round(rho_h.mean(), 4)) + ')')
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) * rho15 / Ltm, Lsky.std(axis=0) * rho15,
label='M2015 (' + str(round(rho15, 4)) + ')', c='violet')
up.add_curve(ax, wl, Lsky.mean(axis=0) * rho99 / Ltm, Lsky.std(axis=0) * rho99, c='orange',
label='M1999(' + str(round(rho99, 4)) + ')')
up.add_curve(ax, wl, Lsky.mean(axis=0) * rho_h / Ltm, Lsky.std(axis=0) * rho_h, c='green',
label='h(' + str(round(rho_h.mean(), 4)) + ')')
ax.set_ylabel(r'$L_{surf}/L_t$')
ax.set_xlabel(r'Wavelength (nm)')
# ---- Lw
ax = axs[1, 0]
up.add_curve(ax, wl, Rrs15.mean(axis=0) * Edm, Rrs15.std(axis=0) * Edm,
label='M2015 (' + str(round(rho15, 4)) + ')', c='violet')
up.add_curve(ax, wl, Rrs99.mean(axis=0) * Edm, Rrs99.std(axis=0) * Edm, c='orange',
label='M1999(' + str(round(rho99, 4)) + ')')
up.add_curve(ax, wl, Rrs_h.mean(axis=0) * Edm, Rrs_h.std(axis=0) * Edm, c='green',
label='h(' + str(round(rho_h.mean(), 4)) + ')')
up.add_curve(ax, wl, Rrs_opt * Edm, Rrs_opt_std * Edm, c='blue',
label='Optimization')
ax.set_ylabel(r'$L_{w}\ (sr^{-1})$')
ax.set_xlabel(r'Wavelength (nm)')
# ---- Rrs
ax = axs[1, 1]
up.add_curve(ax, wl_swr, Rrs_swr.transpose().mean(axis=1), Rrs_swr.transpose().std(axis=1), label='swr',
c='black')
up.add_curve(ax, wl, Rrs15.transpose().mean(axis=1), Rrs15.transpose().std(axis=1),
label='M2015 (' + str(round(rho15, 4)) + ')', c='violet')
up.add_curve(ax, wl, Rrs99.transpose().mean(axis=1), Rrs99.transpose().std(axis=1), c='orange',
label='M1999(' + str(round(rho99, 4)) + ')')
up.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)) + ')')
up.add_curve(ax, wl, Rrs_opt, Rrs_opt_std, c='blue',
label='Optimization')
ax.set_ylabel(r'$R_{rs}\ (sr^{-1})$')
ax.set_xlabel(r'Wavelength (nm)')
ax.set_title('azi=' + str(azi) + ', vza=' + str(vza) + ', sza=' + str(round(sza.mean(), 2)))
# ---- delta Rrs
ax = axs[1, 2]
Rrs_swr_ = interp1d(wl_swr, Rrs_swr.transpose().mean(axis=1), fill_value='extrapolate')(wl)
Rrs_swr_[wl > 850] = np.nan
up.add_curve(ax, wl, (Rrs15.mean(axis=0) - Rrs_swr_) / Rrs_swr_,
label='M2015 (' + str(round(rho15, 4)) + ')', c='violet')
up.add_curve(ax, wl, (Rrs99.mean(axis=0) - Rrs_swr_) / Rrs_swr_, c='orange',
label='M1999(' + str(round(rho99, 4)) + ')')
up.add_curve(ax, wl, (Rrs_h.mean(axis=0) - Rrs_swr_) / Rrs_swr_, c='green',
label='h(' + str(round(rho_h.mean(), 4)) + ')')
up.add_curve(ax, wl, (Rrs_opt - Rrs_swr_) / Rrs_swr_, c='blue',
label='Optimization')
ax.set_ylabel(r'$\Delta^{rel} R_{rs} $')
ax.set_xlabel(r'Wavelength (nm)')
ax.legend(loc='best', frameon=False)
fig.suptitle('trios_awr ' + name + ' idpr' + idpr, fontsize=16)
fig.savefig(os.path.join(dirfig, 'trios_awr_' + name + '_idpr' + idpr + '.png'))
plt.close()
#
#
#
# date = c['Date_prel'].dt.strftime('%Y%m%d')
#
# dbf = '/DATA/projet/petit-saut/data/dataTrios_Guyane_20190523.mdb'
# #data = pda.read_table(dbf,'tblData')
#
#
#
#
# driver='DRIVER={Microsoft Access Driver (*.mdb, *.accdb)}'
# driver='DRIVER={'+pyodbc.drivers()[2]+'}'
# # connect to bd TRIOS
# odbc = pyodbc.connect(driver+';DBQ=' + dbf)
#
# query = 'SELECT * FROM tblData WHERE ((tblData.IDDataType LIKE \'SPECTRUM\') )) ' AND ' \
# '((tblData.IDDataTypeSub1 LIKE \'CALIBRATED\') OR (tblData.IDDataTypeSub1 LIKE \'CALCULATED\')))'
#
# ramses_df = pd.read_sql(query, odbc)
#
#
......@@ -16,10 +16,11 @@ setup(
url='https://gitlab.irstea.fr/ETL-TELQUEL/etl/tree/dev/preprocessing/trios',
license='MIT',
author='T. Harmel',
author_email='tristan.harmel@ntymail.com',
author_email='tristan.harmel@gmail.com',
description='Package to help trios TriOS radiometer data for various above-water or in-water setups',
# TODO update Dependent packages (distributions)
install_requires=['dash','dash_core_components','dash_html_components','pandas', 'scipy', 'numpy', 'netCDF4', 'matplotlib', 'docopt', 'GDAL', 'python-dateutil'],
install_requires=['dash','dash_core_components','dash_html_components','pandas', 'scipy', 'numpy',
'pyodbc', 'netCDF4', 'matplotlib', 'docopt', 'GDAL', 'python-dateutil'],
entry_points={
'console_scripts': [
......
Metadata-Version: 1.0
Name: trios
Version: 1.1.3
Version: 1.1.4
Summary: Package to help trios TriOS radiometer data for various above-water or in-water setups
Home-page: https://gitlab.irstea.fr/ETL-TELQUEL/etl/tree/dev/preprocessing/trios
Author: T. Harmel
Author-email: tristan.harmel@ntymail.com
Author-email: tristan.harmel@gmail.com
License: MIT
Description: UNKNOWN
Platform: UNKNOWN
......@@ -12,6 +12,8 @@ aux/sensors_spec.py
aux/surface_reflectance_factor_rho_coarse_aerosol_rg0.60_sig0.60.csv
aux/surface_reflectance_factor_rho_fine_aerosol_rg0.06_sig0.46.csv
aux/water_coef.txt
exe/__init__.py
exe/db_martinez.py
simulation/__init__.py
simulation/rho_snapshot.py
trios/__init__.py
......
......@@ -4,6 +4,7 @@ dash_html_components
pandas
scipy
numpy
pyodbc
netCDF4
matplotlib
docopt
......
......@@ -10,7 +10,6 @@ import trios.utils.auxdata as ua
from trios.config import *
class awr_process:
def __init__(self, df=None, wl=None):
self.df = df
......@@ -217,7 +216,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, bounds=(0,15), 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)
......@@ -314,7 +313,7 @@ class iwr_process:
self.df = df
self.wl = wl
def process(self, meas, std, mode ='linear'):
def process(self, meas, std, mode='linear'):
wl_ = self.wl
################
......@@ -329,6 +328,8 @@ class iwr_process:
N = len(wl_)
x = meas.prof_Edz # - 0.56
res = uu.fit(N)
if mode == 'lsq':
res = uu.fit(N, 4)
for idx, wl in enumerate(wl_[:-10]):
aw, bbw = iopw.get_iopw(wl)
......@@ -340,47 +341,74 @@ class iwr_process:
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]))
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':
[1.1 * aw, 100], bounds=(
[aw, 0], [np.inf, F0])) # , sigma=sigma, absolute_sigma=True
elif mode == 'lsq':
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))
sig_Edz = self.format_sigma(std.Edz.iloc[:, idx],meas.Edz.iloc[:, idx], 0.1)
sig_Luz = self.format_sigma(std.Luz.iloc[:, idx],meas.Luz.iloc[:, idx], 1e-3)
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]]
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!
res.popt[idx, :], res.pcov[idx, ...] = lsq.x, calc().cov_from_jac(lsq.jac, cost)
if mode == 'lsq':
# TODO formalize and do more clever things for Quality Control
# discard retrieval if error covariance > threshold error covariance median
QC_idx = res.pcov[:, 3,3] > 20 * np.nanmedian(res.pcov[:, 3,3])
res.popt[QC_idx, 3] = np.nan
return res
#@staticmethod
def format_sigma(self, sigma, rescale = 1, noise=0.1):
'''
:param sigma:
:return:
'''
sigma = (sigma + noise) / rescale
return sigma.fillna(np.inf)
# @staticmethod
def f_Edz(self, depth, Kd, Ed0):
'''simple Edz model for homogeneous water column'''
return Ed0 * np.exp(-Kd * depth)
#@staticmethod
# @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
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):
return Lw0minus * np.exp(-KLu * depth)
def cost_func(self, x, z, mes, sigma):
z_Edz = z[0]
z_Lu = z[1]
Edz = mes[0]
Lu = mes[1]
sig_Edz = sigma[0]
sig_Luz = sigma[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)
cost_f1 = (Edz - self.f_Edz(z_Edz, x[0], x[1])) / sig_Edz
cost_f2 = (Lu - self.f_Lu(z_Lu, x[2], x[3])) / sig_Luz
return np.append(cost_f1, cost_f2)
def Kd(self, depth, Edz):
Kd = np.diff(Edz) / np.diff(depth)
......@@ -585,3 +613,30 @@ class calc:
grid=True)
return interp
def cov_from_jac(self, jac, cost):
'''
Compute covariance from jacobian matrix computed in optimization processes
Use Moore-Penrose inverse discarding zero singular values.
:param jac: jacobian
:param cost: cost residual
:return: pcov
'''
from scipy.linalg import svd
M, N = jac.shape
_, s, VT = svd(jac, full_matrices=False)
threshold = np.finfo(float).eps * max(jac.shape) * s[0]
s = s[s > threshold]
VT = VT[:s.size]
pcov = np.dot(VT.T / s ** 2, VT)
if pcov is None:
# indeterminate covariance
pcov = np.zeros((N, N), dtype=float)
pcov.fill(np.inf)
else:
s_sq = cost / (M - N)
pcov = pcov * s_sq
return pcov
import glob
import re
import matplotlib as mpl
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 18})
from scipy.interpolate import interp1d
......@@ -12,20 +14,21 @@ coordf = glob.glob("/DATA/OBS2CO/data/info/mesures_in_situ.csv")[0]
coords = pd.read_csv(coordf, sep=';')
dirfig = os.path.abspath('/DATA/OBS2CO/data/trios/fig')
awrfiles = glob.glob("/DATA/OBS2CO/data/trios/raw/aw*idpr*.csv")
awrfiles = glob.glob("/DATA/OBS2CO/data/trios/raw/2018/aw*idpr*.csv")
# awrfiles = glob.glob("/DATA/OBS2CO/data/trios/test_setup/raw/aw*idpr*.csv")
swrfiles = glob.glob("/DATA/OBS2CO/data/trios/raw/Lu0*idpr*.csv")
swrfiles = glob.glob("/DATA/OBS2CO/data/trios/raw/2018/Lu0*idpr*.csv")
iopw = ua.iopw()
iopw.load_iopw()
def add_curve(ax, x, mean, std, c='red', label=''):
def add_curve(ax, x, mean, std=None, c='red', label='',**kwargs):
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)
alpha=0.8, label=label,*kwargs)
if np.any(std):
ax.fill_between(x,
mean - std,
mean + std, alpha=0.35, color=c)
idpr = '167'
......@@ -41,8 +44,7 @@ for idpr in idprs:
alt = 0 # c['Altitude']
name = c['ID_lac'].values[0]
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(10, 8))
fig.subplots_adjust(left=0.1, right=0.9, hspace=.5, wspace=0.65)
# -----------------------------------------------
# SWR processing
......@@ -58,9 +60,6 @@ for idpr in idprs:
df.at[index, 'sza'] = sza
swr = swr_process(df, wl_swr)
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')
# -----------------------------------------------
# AWR processing
......@@ -137,15 +136,75 @@ for idpr in idprs:
# ------------------
# plotting
# ------------------
Ltm = Lt.mean(axis=0)
Edm = Ed.mean(axis=0)
mpl.rcParams.update({'font.size': 18})
fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(20, 12))
fig.subplots_adjust(left=0.1, right=0.9, hspace=.5, wspace=0.45)
# ---- Ed
ax = axs[0,0]
add_curve(ax, wl, Ed.mean(axis=0),
label=r'$L_{sky}$',c='red') # just to put the two labels
add_curve(ax, wl, Ed.mean(axis=0), Ed.std(axis=0),
label=r'$E_s$',c='black')
ax.set_ylabel(r'$E_{d}(0^{+})$')