From 6b4880fc14441934c6a4f9df8fac36f84dcec4d1 Mon Sep 17 00:00:00 2001
From: Tristan Harmel <tristan.harmel@get.omp.eu>
Date: Tue, 12 Mar 2019 11:06:30 +0100
Subject: [PATCH] add 'run' file for swr, define output format
---
build/lib/process/process.py | 225 ++++++++++++++++++++----
build/lib/process/process_compar_awr.py | 49 ++++--
build/lib/process/process_test_setup.py | 2 +-
build/lib/utils/utils.py | 100 ++++++-----
build/lib/visu/data_visu.py | 73 ++++----
process/run_iwr.py | 16 +-
process/run_swr.py | 86 +++++++++
setup.py | 2 +-
trios.egg-info/PKG-INFO | 2 +-
trios.egg-info/SOURCES.txt | 6 +
trios.egg-info/requires.txt | 3 +
trios.egg-info/top_level.txt | 2 +
12 files changed, 419 insertions(+), 147 deletions(-)
create mode 100644 process/run_swr.py
diff --git a/build/lib/process/process.py b/build/lib/process/process.py
index 6e6f4a7..3c25eb7 100644
--- a/build/lib/process/process.py
+++ b/build/lib/process/process.py
@@ -1,11 +1,14 @@
import numpy as np
import pandas as pd
-import scipy.interpolate as si
+from scipy import interpolate, integrate
+from scipy.optimize import curve_fit
+
import plotly.plotly as py
# import plotly.graph_objs as go
from plotly.graph_objs import *
from utils.utils import reshape as r
+import utils.auxdata as ua
from config import *
@@ -30,25 +33,33 @@ class awr_process:
self.rhoM1999.dropna(inplace=True)
self.rhoM2015.dropna(inplace=True)
- def get_rho_values(self, sza, vza, azi, ws=2, aot=0.1, wl=[550]):
+ def get_rho_values(self, sza, vza, azi, ws=[2], aot=[0.1], wl=[550], sunglint=False):
'''
-
- :param sza:
- :param vza:
- :param azi:
- :param ws:
- :param aot:
- :param wl:
+ Interpolate the rho factor values from tabulated data
+
+ :param sza: solar zenith angle in deg, array-like
+ :param vza: view zenith angle in deg, array-like
+ :param azi: relative azimuth in deg (=0 when looking at Sun), array-like
+ :param ws: wind speed, m/s,(based on Cox-Munk parametrization of surface roughness) array-like
+ :param aot: aerosol optical thickness at 550 nm, array-like
+ :param wl: wavelength in nm, array-like
+ :param sunglint: add sunglint component in rho calculation if True
:return:
'''
+ grid = self.rho.rho.index.levels
-
- grid = self.rho.rho.index.levels[2:]
# convert pandas dataframe into 6D array of the tabulated rho values for interpolation
- rho_ = r().df2ndarray(self.rho, 'rho')
+ rhoname = 'rho'
+ if sunglint:
+ rhoname = 'rho_g'
+
+ rho_6d = r().df2ndarray(self.rho, rhoname)
+
+ rho_ = calc().spline_2d(grid[-2:], rho_6d, (azi, vza))
+
+ rho_wl = calc().spline_4d(grid[:-2], rho_, (ws, aot, wl, sza))
- rho_wl = calc().spline_4d(grid, rho_[1, 1, ...], (wl, sza, azi, vza))
return rho_wl.squeeze()
def get_rho_mobley(self, rhodf, sza, vza, azi, ws):
@@ -78,10 +89,13 @@ class awr_process:
Lt = self.df.loc[:, ("Lt")]
Lsky = self.df.loc[:, ("Lsky")]
Ed = self.df.loc[:, ("Ed")]
- sza = self.df.loc[:, ("sza")]
- return self.process(wl, Lt, Lsky, Ed, sza, vza, azi, ws, aot)
+ sza = self.df.loc[:, ("sza")].values.mean()
+ Rrs = self.process(wl, Lt, Lsky, Ed, sza, vza, azi, ws, aot)
+ Rrs.columns = pd.MultiIndex.from_product([['Rrs(awr)'], self.Rrs.columns], names=['param', 'wl'])
+ self.Rrs = Rrs
+ return self.Rrs
- def process(self, wl, Lt, Lsky, Ed, sza, vza=[40], azi=[135], ws=2, aot=0.1):
+ def process(self, wl, Lt, Lsky, Ed, sza, vza=[40], azi=[135], ws=[2], aot=[0.1]):
'''
:param wl:
@@ -96,9 +110,8 @@ class awr_process:
:return:
'''
- rho = self.get_rho_values([sza.values.mean()],vza,azi,wl=wl )
+ rho = self.get_rho_values([sza], vza, azi, wl=wl, ws=ws, aot=aot)
self.Rrs = (Lt - rho * Lsky) / Ed
- self.Rrs.columns = pd.MultiIndex.from_product([['Rrs(awr)'], self.Rrs.columns], names=['param', 'wl'])
return self.Rrs, rho
@@ -108,13 +121,78 @@ class swr_process:
self.df = df
self.wl = wl
- def process(self):
+ def call_process(self, shade_corr=False):
wl = self.wl
- df = self.df
-
- Rrs = df.loc[:, ("Lu0+")] / df.loc[:, ("Ed")]
+ Lu = self.df.loc[:, ("Lu0+")]
+ Ed = self.df.loc[:, ("Ed")]
+ sza = self.df.loc[:, ("sza")].values.mean()
+ Rrs = self.process(Lu, Ed, sza, wl, shade_corr=shade_corr)
Rrs.columns = pd.MultiIndex.from_product([['Rrs(swr)'], Rrs.columns], names=['param', 'wl'])
self.Rrs = Rrs
+ return Rrs
+
+ def process(self, Lu, Ed, sza, wl, R=0.05, shade_corr=False):
+ Rrs = Lu / Ed
+ ang_w = calc().angle_w(sza)
+
+ iopw = ua.iopw()
+ iopw.load_iopw()
+ iopw.get_iopw(wl)
+ a, bb = iopw.aw, iopw.bbw
+ # TODO add particulate and dissolved component to a and bb values
+ # a,bb = aux.get_iop(..., withwater=True)
+ acdom = ua.cdom(0.5, wl).get_acdom()
+ a = a + acdom + 0.4
+ bb = bb + 0.05
+ if shade_corr:
+ Rrs = self.shade_corr(Rrs, R, ang_w, a, bb, wl)
+ # Rrs.columns = pd.MultiIndex.from_product([['Rrs(swr)'], Rrs.columns], names=['param', 'wl'])
+ self.Rrs = Rrs
+ self.a = a
+ self.bb = bb
+ self.acdom = acdom
+ return self.Rrs
+
+ def epsilon(self, K, R, ang_w):
+ '''
+ epsilon from Shang et al, 2017, Applied Optics
+ :param K:
+ :param R:
+ :param ang_w: Sun zenith angle below surface (in deg)
+ :return:
+ '''
+
+ self.eps = np.array(1 - np.exp(-K * R / np.tan(np.radians(ang_w))))
+ return self.eps
+
+ def K(self, a, bb, ang_w):
+ '''
+ K (sum attenuation coef. of Lu in and outside the shade) from Shang et al, 2017, Applied Optics
+ :param a: total absorption coefficient (m-1)
+ :param bb: total backscattering coefficient (m-1)
+ :param ang_w: Sun zenith angle below surface (in deg)
+ :return:
+ '''
+ sin_ang_w = np.sin(np.radians(ang_w))
+ self.K_ = (3.15 * sin_ang_w + 1.15) * a * np.exp(-1.57 * bb) \
+ + (5.62 * sin_ang_w - 0.23) * bb * np.exp(-0.5 * a)
+ return self.K_
+
+ def shade_corr(self, Rrs, R, ang_w, a, bb, wl, wl_cutoff=900):
+ '''
+ Correction of shading error from Shang et al, 2017, Applied Optics
+ :param Rrs:
+ :param R:
+ :param ang_w:
+ :param a:
+ :param bb:
+ :return:
+ '''
+
+ K = self.K(a, bb, ang_w)
+ eps = self.epsilon(K, R, ang_w)
+ eps[wl > wl_cutoff] = 0
+ self.Rrs = Rrs / (1 - eps)
return self.Rrs
@@ -127,18 +205,33 @@ class iwr_process:
wl = self.wl
df = self.df
- Rrs = df.loc[:, ("Lu0+")] / df.loc[:, ("Ed")]
- Rrs.columns = pd.MultiIndex.from_product([['Rrs(swr)'], Rrs.columns], names=['param', 'wl'])
- self.Rrs = Rrs
- return self.Rrs
+ reflectance = df.loc[:, ("Luz")] / df.loc[:, ("Edz")]
+ reflectance.columns = pd.MultiIndex.from_product([['reflectance'], reflectance.columns], names=['param', 'wl'])
+ self.reflectance = reflectance
+
+ df['rounded_depth', ''] = df.prof_Edz.round(1)
+ df.groupby('rounded_depth').mean()
+
+ return self.reflectance
+
+ @staticmethod
+ def f_Edz(depth, Kd, Ed0):
+ '''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)
- def plot_raw(self):
+ def plot_raw(self,x='Luz',y='prof_Luz'):
trace = Scattergl(
- x=self.df['Luz'].values,
- y=self.df['depth_Luz'].values,
+ x=self.df[x].values,
+ y=self.df[y].values,
text=self.df.index.get_level_values(0),
hoverinfo="text",
@@ -172,8 +265,8 @@ class iwr_process:
title=''
),
yaxis=dict(
- range=[min(-5, min(self.df['depth_Luz'])),
- max(0, max(self.df['depth_Luz']))],
+ range=[min(-5, min(self.df[y])),
+ max(0, max(self.df[y]))],
showline=False,
fixedrange=True,
zeroline=False,
@@ -200,7 +293,7 @@ class self_shading:
self.eps_dif_EuZ = 2.70
def epsilon(self, sza):
- eps = si.interp1d(self.ang, self.eps_dif_EuZ)(sza)
+ eps = interpolate.interp1d(self.ang, self.eps_dif_EuZ)(sza)
return eps
@@ -208,6 +301,33 @@ class calc:
def __init__(self):
pass
+ def PAR(self, wl, Ed):
+ '''
+ Compute instantaneous PAR from Ed spectrum.
+ PAR in mW m-2
+ PAR_quanta in µmol photon m-2 s-1
+ :param wl:
+ :param Ed:
+ :return:
+ '''
+ # ---------------------------------------------
+ # PARAMETERS
+ # Planck constant in J s or W s2
+ h = 6.6260695729e-3 # d-34
+ # light speed in m s-1
+ c = 2.99792458e0 # d8
+ # Avogadro Number in mol-1
+ Avogadro = 6.0221412927e0 # d23
+ hc = Avogadro * h * c
+ # ---------------------------------------------
+
+ idx_par = (wl >= 400) & (wl <= 700)
+ wl = wl[idx_par]
+ Ed = Ed[idx_par]
+ par = integrate.trapz(Ed, wl)
+ par_quanta = integrate.trapz(np.multiply(wl,Ed), wl) / hc
+ return par, par_quanta
+
def earth_sun_correction(self, dayofyear):
'''
Earth-Sun distance correction factor for adjustment of mean solar irradiance
@@ -226,6 +346,40 @@ class calc:
return bidir
+ def angle_w(self, angle_air, n=1.334):
+ '''
+ convert above surface angle (angle_air) into sub-surface angle
+ :param angle_air in deg
+ :param n: refractive index of water
+ :return: sub-surface angle in deg
+ '''
+ return np.degrees(np.arcsin(np.sin(np.radians(angle_air)) / n))
+
+ def spline_2d(self, gin, arr, gout):
+ '''
+ Interpolation of a 6D array (arr) with bicubic splines on a 2D grid
+ corresponding to the 5th and 6th dimensions of arr.
+ Return 4D array interpolated on gout.
+
+ :param gin: regular 2D grid of the tabulated data (tuple/array/list of arrays)
+ :param arr: tabulated data (N dimensions, interpolation on N-1 and N)
+ :param gout: new 2D grid on which data are interpolated (with dims 2 and 3 of the same length);
+ (tuple/array/list of arrays)
+ :return: Interpolated data (1D or 3D array depending on the dimension shapes of gout
+ '''
+
+ N = arr.shape
+ interp = np.zeros(N[:-2])
+
+ for i in range(N[0]):
+ for j in range(N[1]):
+ for k in range(N[2]):
+ for l in range(N[3]):
+ interp[i, j, k, l] = interpolate.RectBivariateSpline(gin[0], gin[1], arr[i, j, k, l, ...])(
+ gout[0], gout[1], grid=False)
+
+ return interp
+
def spline_4d(self, gin, lut, gout):
'''
Interpolation with two successive bicubic splines on a regular 4D grid.
@@ -239,7 +393,6 @@ class calc:
(tuple/array/list of arrays)
:return: Interpolated data (1D or 3D array depending on the dimension shapes of gout
'''
- import scipy.interpolate as si
N = gin[0].__len__(), gin[1].__len__(), gin[2].__len__(), gin[3].__len__()
Nout = gout[0].__len__(), gout[1].__len__(), gout[2].__len__()
@@ -247,15 +400,15 @@ class calc:
for i in range(N[0]):
for j in range(N[1]):
- tmp[i, j, :] = si.RectBivariateSpline(gin[2], gin[3], lut[i, j, :, :])(gout[2], gout[3], grid=False)
+ tmp[i, j, :] = interpolate.RectBivariateSpline(gin[2], gin[3], lut[i, j, :, :])(gout[2], gout[3], grid=False)
if Nout[0] == Nout[1] == 1:
interp = np.ndarray(Nout[2])
for iband in range(Nout[2]):
- interp[iband] = si.RectBivariateSpline(gin[0], gin[1], tmp[:, :, iband])(gout[0], gout[1], grid=False)
+ interp[iband] = interpolate.RectBivariateSpline(gin[0], gin[1], tmp[:, :, iband])(gout[0], gout[1], grid=False)
else:
interp = np.ndarray([Nout[0], Nout[1], Nout[2]])
for iband in range(Nout[2]):
- interp[:, :, iband] = si.RectBivariateSpline(gin[0], gin[1], tmp[:, :, iband])(gout[0], gout[1],
+ interp[:, :, iband] = interpolate.RectBivariateSpline(gin[0], gin[1], tmp[:, :, iband])(gout[0], gout[1],
grid=True)
return interp
diff --git a/build/lib/process/process_compar_awr.py b/build/lib/process/process_compar_awr.py
index 9ccfa92..6fae7a2 100644
--- a/build/lib/process/process_compar_awr.py
+++ b/build/lib/process/process_compar_awr.py
@@ -13,6 +13,7 @@ from scipy.interpolate import interp1d
from utils.sunposition import sunpos
import utils.utils as u
+import utils.auxdata as ua
from process.process import *
coordf = glob.glob("/DATA/OBS2CO/data/info/mesures_in_situ.csv")[0]
@@ -23,6 +24,8 @@ awrfiles = glob.glob("/DATA/OBS2CO/data/trios/raw/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")
+iopw = ua.iopw()
+iopw.load_iopw()
def add_curve(ax, x, mean, std, c='red', label=''):
ax.plot(x, mean, linestyle='solid', c=c, lw=2.5,
@@ -36,7 +39,7 @@ idpr = '167'
# get idpr numbers
idprs = np.unique([re.findall(r'idpr(\d+)', x)[0] for x in swrfiles])
-idprs = np.array(['156'])
+#idprs = np.array(['170'])
# loop over idpr
for idpr in idprs:
c = coords[coords.ID_prel == int(idpr)] # .values[0]
@@ -52,23 +55,29 @@ for idpr in idprs:
# SWR processing
# -----------------------------------------------
- swr = u.swr_data(idpr, swrfiles)
- if swr.file:
- df, wl_swr = swr.reader(lat, lon, alt)
- Rrs_swr = swr_process(df, wl_swr).process()
- add_curve(ax, wl_swr, Rrs_swr.transpose().mean(axis=1), Rrs_swr.transpose().std(axis=1), label='swr', c='black')
-
+ 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()
+ 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
# -----------------------------------------------
+ azi = 135
+ vza = 40
awr = u.awr_data(idpr, awrfiles)
if awr.Edf:
-
index_idx = [0]
- azi = 135
- vza = 40
d = u.data(index_idx)
Ed, wl_Ed = d.load_csv(awr.Edf)
@@ -115,9 +124,9 @@ for idpr in idprs:
sza = sunpos(index, lat, lon, alt)[1]
df.at[index, 'sza'] = sza
- rho_h = awr.get_rho_values([df.sza.min()], [vza], [azi], wl=wl)
- rho15 = awr.get_rho_mobley(awr.rhoM2015, [df.sza.min()], [vza], [azi], [ws])
- rho99 = awr.get_rho_mobley(awr.rhoM1999, [df.sza.min()], [vza], [azi], [ws])
+ rho_h = awr.get_rho_values([df.sza.mean()], [vza], [azi], wl=wl)
+ rho15 = awr.get_rho_mobley(awr.rhoM2015, [df.sza.mean()], [vza], [azi], [ws])
+ rho99 = awr.get_rho_mobley(awr.rhoM1999, [df.sza.mean()], [vza], [azi], [ws])
Rrs_h = (df.loc[:, 'Lt'] - rho_h * df.loc[:, 'Lsky']) / df.loc[:, 'Ed']
Rrs15 = (df.loc[:, 'Lt'] - rho15 * df.loc[:, 'Lsky']) / df.loc[:, 'Ed']
@@ -125,12 +134,14 @@ for idpr in idprs:
Rrs99 = (df.loc[:, 'Lt'] - rho99 * df.loc[:, 'Lsky']) / df.loc[:, 'Ed']
# plt.figure()
- add_curve(ax, wl, Rrs15.transpose().mean(axis=1), Rrs15.transpose().std(axis=1), label='M2015')
- add_curve(ax, wl, Rrs99.transpose().mean(axis=1), Rrs99.transpose().std(axis=1), c='orange', label='M1999')
- add_curve(ax, wl, Rrs_h.transpose().mean(axis=1), Rrs_h.transpose().std(axis=1), c='grey', label='h')
+ add_curve(ax, wl, Rrs15.transpose().mean(axis=1), Rrs15.transpose().std(axis=1),
+ label='M2015 (' + str(rho15) + ')')
+ add_curve(ax, wl, Rrs99.transpose().mean(axis=1), Rrs99.transpose().std(axis=1), c='orange',
+ label='M1999(' + str(rho99) + ')')
+ add_curve(ax, wl, Rrs_h.transpose().mean(axis=1), Rrs_h.transpose().std(axis=1), c='grey',
+ label='h(' + str(rho_h.mean()) + ')')
-
- ax.set_title('azi=' + str(azi) + ', vza=' + str(vza)+', sza='+str(sza))
+ ax.set_title('azi=' + str(azi) + ', vza=' + str(vza) + ', sza=' + str(sza))
ax.legend(loc='best', frameon=False)
@@ -138,4 +149,4 @@ for idpr in idprs:
ax.set_xlabel(r'Wavelength (nm)')
fig.savefig(os.path.join(dirfig, 'trios_awr_' + name + '_idpr' + idpr + '.png'))
plt.close()
- Lt.index.names
+
diff --git a/build/lib/process/process_test_setup.py b/build/lib/process/process_test_setup.py
index 55d6cae..b3bf378 100644
--- a/build/lib/process/process_test_setup.py
+++ b/build/lib/process/process_test_setup.py
@@ -71,7 +71,7 @@ fig.subplots_adjust(left=0.1, right=0.9, hspace=.5, wspace=0.65)
i=0
for azi, Lt1 in Lt0.groupby(level=2):
for vza,Lt in Lt1.groupby(level=1):
- ax=axs.flat[i]
+ ax = axs.flat[i]
i=i+1
print(azi,vza)
diff --git a/build/lib/utils/utils.py b/build/lib/utils/utils.py
index a2b7bca..53ecc69 100644
--- a/build/lib/utils/utils.py
+++ b/build/lib/utils/utils.py
@@ -10,7 +10,7 @@ class awr_data:
Above-water radiometry
'''
- def __init__(self, idpr, files):
+ def __init__(self, idpr=None, files=None):
# ''' get file names for Ed, Lsky and Lt data'''
self.file = list(filter(lambda x: 'idpr' + idpr in x, files))
file = self.file
@@ -41,7 +41,7 @@ class awr_data:
df = pd.DataFrame()
# ''' read files with pandas format '''
- d=data(index_idx)
+ d = data(index_idx)
Ed, wl_Ed = d.load_csv(self.Edf)
Lsky, wl_Lsky = d.load_csv(self.Lskyf)
Lt, wl_Lt = d.load_csv(self.Ltf)
@@ -90,7 +90,7 @@ class iwr_data:
self.Luzf = list(filter(lambda x: 'Luz' in x, file))
self.idpr = idpr
- def reader(self, lat=None, lon=None, alt=0, name=''):
+ def reader(self, lat, lon, alt=0, name='', delta_Lu_depth=0, delta_Edz_depth=0):
'''
Read above-water data files for a given acquisition series (idpr),
merge the different data types:
@@ -105,6 +105,9 @@ class iwr_data:
:param lat: latitude (decimal)
:param lon: longitude (decimal)
:param alt: altitude (m)
+ :param delta_Lu_depth: adjustment of actual depth for Lu sensor (distance from depth sensor);
+ in meters for depth counted positively
+ :param delta_Edz_depth: similar to delta_Lu_depth for Edz sensor
:param idpr: ID of the acquisition series
:return:
'''
@@ -112,29 +115,42 @@ class iwr_data:
df = pd.DataFrame()
# ''' read files with pandas format '''
- Ed, wl_Ed = self.load_csv(self.Edf)
- Edz, wl_Edz = self.load_csv(self.Edzf)
- Luz, wl_Luz = self.load_csv(self.Luzf)
+ d = data([1, 0])
+
+ Ed, wl_Ed = d.load_csv(self.Edf)
+ Edz, wl_Edz = d.load_csv(self.Edzf)
+ Luz, wl_Luz = d.load_csv(self.Luzf)
+
+ #mask negative values TODO save number of discarded data
+ # Ed.mask(Ed<0,inplace=True)
+ # Edz.mask(Edz<0,inplace=True)
+ # Luz.mask(Luz<0,inplace=True)
+
+
# copy depth data to Ed frame on date index
- Ed.index = Ed.index.droplevel(level=1)
+ # Ed.index = Ed.index.droplevel(level=1)
- # ''' interpolate Ed and Lsky data upon Lt wavelength'''
+ #''' 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,
- columns=pd.MultiIndex.from_tuples(zip(['Ed'] * len(wl), wl), names=['param', 'wl']),
+ newEd = pd.DataFrame(index=Ed.index.get_level_values(0),
+ 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)
- # merge sensor data on time
+ # correct depth data for sensor to sensor distance
Luz.reset_index(level=1, inplace=True)
+ Luz.iloc[:, 0] = Luz.iloc[:, 0] + delta_Lu_depth
# newEd.reset_index(level=1,inplace=True)
newEdz.reset_index(level=1, inplace=True)
+ newEdz.iloc[:, 0] = newEdz.iloc[:, 0] + delta_Edz_depth
+
+ # merge sensor data on time
df = pd.merge_asof(Luz, newEd, left_index=True, right_index=True, tolerance=pd.Timedelta("2 seconds"),
direction="nearest")
df = pd.merge_asof(df, newEdz, left_index=True, right_index=True, suffixes=('_Luz', '_Edz'),
@@ -150,26 +166,26 @@ class iwr_data:
df.at[index, 'sza'] = sza
df['idpr', ''] = self.idpr
- df['name', ''] = self.name
+ df['name', ''] = name
return df, wl
- def load_csv(self, file):
-
- dateparse = lambda x: pd.datetime.strptime(x, '%Y-%m-%d %H:%M:%S')
- if len(file) > 1:
- print('Warning! Multiple files found but only one expected, process first file of the list:')
- print(file)
- file = file[0]
- df = pd.read_csv(file, sep=';', index_col=[1, 0], na_values=['-NAN'])
- df = df.dropna(axis=1, how='all').dropna(axis=0, how='all')
- df.index = df.index.set_levels([pd.to_datetime(df.index.levels[0]), df.index.levels[1]])
- df.columns = df.columns.astype('float') # str.extract('(\d+)',expand=False).astype('float')
- # resort to get data in increasing time order
- df.sort_index(inplace=True, level=0)
- wl = df.columns
-
- return df, wl
+ # def load_csv(self, file):
+ #
+ # dateparse = lambda x: pd.datetime.strptime(x, '%Y-%m-%d %H:%M:%S')
+ # if len(file) > 1:
+ # print('Warning! Multiple files found but only one expected, process first file of the list:')
+ # print(file)
+ # file = file[0]
+ # df = pd.read_csv(file, sep=';', index_col=[1, 0], na_values=['-NAN'])
+ # df = df.dropna(axis=1, how='all').dropna(axis=0, how='all')
+ # df.index = df.index.set_levels([pd.to_datetime(df.index.levels[0]), df.index.levels[1]])
+ # df.columns = df.columns.astype('float') # str.extract('(\d+)',expand=False).astype('float')
+ # # resort to get data in increasing time order
+ # df.sort_index(inplace=True, level=0)
+ # wl = df.columns
+ #
+ # return df, wl
class swr_data:
@@ -235,28 +251,28 @@ class swr_data:
class data:
- def __init__(self,index_idx=[0]):
+ def __init__(self, index_idx=[0]):
# first position should be datetime index
# followed by the other parameters used for indexing (e.g. azimuth, view angle)
- self.index_idx=index_idx
+ self.index_idx = index_idx
pass
def load_csv(self, file):
+ print(file)
dateparse = lambda x: pd.datetime.strptime(x, '%Y-%m-%d %H:%M:%S')
if len(file) > 1:
- print( 'Warning! Multiple files found but only one expected, process first file of the list:')
+ print('Warning! Multiple files found but only one expected, process first file of the list:')
print(file)
- file = file[0]
- #df = pd.read_csv(file, date_parser=dateparse, sep=';', index_col=0, na_values=['-NAN'])
- df=pd.read_csv(file, sep=';', na_values=['-NAN'])
- df = df.dropna(axis=1, how='all').dropna(axis=0, how='all')
-
- #get list of indexes
- col=df.columns.values[self.index_idx]
- df[col[0]]=pd.to_datetime(df[col[0]])
+ file_ = file[0]
+ # df = pd.read_csv(file, date_parser=dateparse, sep=';', index_col=0, na_values=['-NAN'])
+ df = pd.read_csv(file_, sep=';', na_values=['-NAN'])
- df.set_index(col.tolist(),inplace=True)
+ # get list of indexes
+ col = df.columns.values[self.index_idx]
+ df[col[0]] = pd.to_datetime(df[col[0]])
+ df.set_index(col.tolist(), inplace=True)
+ df = df.dropna(axis=1, how='all').dropna(axis=0, how='all')
df.columns = df.columns.astype('float') # str.extract('(\d+)',expand=False).astype('float')
# resort to get data in increasing time order
df.sort_index(inplace=True)
diff --git a/build/lib/visu/data_visu.py b/build/lib/visu/data_visu.py
index 72a35f8..539a7b0 100644
--- a/build/lib/visu/data_visu.py
+++ b/build/lib/visu/data_visu.py
@@ -13,6 +13,7 @@ import pandas as pd
import plotly.graph_objs as go
from dash.dependencies import Input, Output
+
def main():
app = dash.Dash()
# app.css.append_css('data/aeronet_layout.css')
@@ -58,9 +59,9 @@ def main():
'width': '100%', 'display': 'inline-block',
}),
html.Div([
- html.H4('File...', id='filename',style={'float': 'left','width': '60%'}),
+ html.H4('File...', id='filename', style={'float': 'left', 'width': '60%'}),
html.Div([
- html.H4('Color variable:',style={'margin-bottom': '0','width':'50%'}),
+ html.H4('Color variable:', style={'margin-bottom': '0', 'width': '50%'}),
html.Div([
dcc.Dropdown(
id='color-column',
@@ -73,7 +74,7 @@ def main():
# value='Linear',
# labelStyle={'width': '25%','display': 'inline-block'})],
],
- style={'width': '50%','display': 'inline-block'})],
+ style={'width': '50%', 'display': 'inline-block'})],
style={'width': '40%',
'margin-top': '0',
'display': 'inline-block',
@@ -101,8 +102,8 @@ def main():
value='Lt'),
],
style={'width': '48.9%',
- 'float': 'right', }),],
- style={'width':'100%','margin-block-start': '1%'}),
+ 'float': 'right', }), ],
+ style={'width': '100%', 'margin-block-start': '1%'}),
html.Div([
html.Div([
@@ -123,7 +124,7 @@ def main():
}),
],
- style={'height': '20%'},
+ style={'height': '20%'},
# className='row'
),
@@ -191,27 +192,26 @@ def main():
'margin-left': 'auto',
'margin-right': 'auto'})
-
def figure_spectrum(df, column_name, color_column_name):
dff = df
- layout = go.Layout( xaxis={'title': 'Wavelength (nm)'},
- yaxis={'title': column_name},
- margin={'l': 50, 'b': 40, 't': 20, 'r': 50},
- hovermode='closest',
- height=300,
- font=dict(color='#CCCCCC'),
- titlefont=dict(color='#CCCCCC', size=14),
- plot_bgcolor="#191A1A",
- paper_bgcolor="#020202")
-
- if not(column_name in dff.columns.get_level_values(0)):
- return {'data': [],
- 'layout': layout}
+ layout = go.Layout(xaxis={'title': 'Wavelength (nm)'},
+ yaxis={'title': column_name},
+ margin={'l': 50, 'b': 40, 't': 20, 'r': 50},
+ hovermode='closest',
+ height=300,
+ font=dict(color='#CCCCCC'),
+ titlefont=dict(color='#CCCCCC', size=14),
+ plot_bgcolor="#191A1A",
+ paper_bgcolor="#020202")
+
+ if not (column_name in dff.columns.get_level_values(0)):
+ return {'data': [],
+ 'layout': layout}
parameters = df.loc[:, (color_column_name)].values
dff = dff.loc[:, ([column_name])]
wl = dff.columns.droplevel().values
- dff = dff.stack(level=['wl'])
+ #dff = dff.stack(level=['wl'])
norm = mpl.colors.Normalize(vmin=np.nanmin(parameters), vmax=np.nanmax(parameters))
# create a ScalarMappable and initialize a data structure
cmap = cm.Spectral
@@ -219,11 +219,14 @@ def main():
s_m.set_array([])
i = 0
trace = []
- for date, x in dff.groupby(level=0):
+ #for date, x_ in dff.groupby(level=0):
+ for idx, x in dff.iterrows():
+ date = x.name.__str__()
+ #print(idx,x)
trace.append(go.Scattergl(
x=wl, # spectrum,
- y=x[column_name].values.flatten(),
- text=dff.index.get_level_values(0),
+ y=x[column_name],
+ text=date,#str(parameters[i]),#dff.index.get_level_values(0),
name=str(parameters[i]),
mode='lines',
marker={
@@ -244,7 +247,6 @@ def main():
'layout': layout
}
-
def list_data(contents, level=0):
# content_type, content_string = contents.split(',')
# decoded = base64.b64decode(content_string)
@@ -253,6 +255,7 @@ def main():
# remove useless variables
c = c.drop(filter(lambda s: re.match('.*(Wave|Tri|[sS]ite|Dat)', s), c))
return [{'label': i, 'value': i} for i in c]
+
#
def parse_contents(contents):
global df
@@ -275,7 +278,6 @@ def main():
parse_contents(contents)
return contents
-
@app.callback(Output('filename', 'children'),
[Input('upload-data', 'filename')])
def show_filename(filename):
@@ -289,68 +291,59 @@ def main():
def update_dropdown(contents):
return list_data(contents)
-
@app.callback(Output('spectrum1', 'options'),
[Input('dataset', 'children')])
def update_dropdown(contents):
return list_data(contents, level=0)
-
@app.callback(Output('spectrum2', 'options'),
[Input('dataset', 'children')])
def update_dropdown(contents):
return list_data(contents, level=0)
-
@app.callback(Output('spectrum3', 'options'),
[Input('dataset', 'children')])
def update_dropdown(contents):
return list_data(contents, level=0)
-
@app.callback(Output('spectrum4', 'options'),
[Input('dataset', 'children')])
def update_dropdown(contents):
return list_data(contents, level=0)
-
# selection from time series graph -> spectrum graph
@app.callback(Output('graph1', 'figure'),
[Input('spectrum1', 'value'),
Input('color-column', 'value'),
Input('dataset', 'children')])
- def spectrum_figure(column_name, color_column_name,void):
+ def spectrum_figure(column_name, color_column_name, void):
return figure_spectrum(df, column_name, color_column_name)
-
# selection from time series graph -> spectrum graph
@app.callback(Output('graph2', 'figure'),
[Input('spectrum2', 'value'),
Input('color-column', 'value'),
Input('dataset', 'children')])
- def spectrum_figure(column_name, color_column_name,void):
+ def spectrum_figure(column_name, color_column_name, void):
return figure_spectrum(df, column_name, color_column_name)
-
@app.callback(Output('graph3', 'figure'),
[Input('spectrum3', 'value'),
Input('color-column', 'value'),
Input('dataset', 'children')])
- def spectrum_figure(column_name, color_column_name,void):
+ def spectrum_figure(column_name, color_column_name, void):
return figure_spectrum(df, column_name, color_column_name)
-
# selection from time series graph -> spectrum graph
@app.callback(Output('graph4', 'figure'),
[Input('spectrum4', 'value'),
Input('color-column', 'value'),
Input('dataset', 'children')])
- def spectrum_figure(column_name, color_column_name,void):
+ def spectrum_figure(column_name, color_column_name, void):
return figure_spectrum(df, column_name, color_column_name)
-
app.run_server(port=8060)
+
if __name__ == "__main__":
main()
-
diff --git a/process/run_iwr.py b/process/run_iwr.py
index e5d120e..357695a 100644
--- a/process/run_iwr.py
+++ b/process/run_iwr.py
@@ -32,7 +32,7 @@ class fit:
# ------------------------------------------------
# above-water data files
dirfig = os.path.abspath('/DATA/OBS2CO/data/trios/fig')
-dirout = os.path.abspath('/DATA/OBS2CO/data/trios/above_water')
+dirout = os.path.abspath('/DATA/OBS2CO/data/trios/in_water')
iwrfiles = glob.glob("/DATA/OBS2CO/data/trios/raw/uw*idpr*.csv")
@@ -42,6 +42,14 @@ coords
# get idpr numbers
idprs = np.unique([re.findall(r'idpr(\d+)', x)[0] for x in iwrfiles])
+
+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)
+
################
# load aux data
iopw = ua.iopw()
@@ -142,12 +150,6 @@ for idpr in idprs:#[-1:]:
plt.close()
- 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)
mpl.rcParams.update({'font.size': 18})
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(15, 6))
diff --git a/process/run_swr.py b/process/run_swr.py
new file mode 100644
index 0000000..4d6f25c
--- /dev/null
+++ b/process/run_swr.py
@@ -0,0 +1,86 @@
+import base64
+import pandas as pd
+import numpy as np
+import glob
+import io
+import os
+from textwrap import dedent as d
+import re
+import matplotlib.pyplot as plt
+import plotly
+import plotly.graph_objs as go
+from scipy.interpolate import interp1d
+
+from utils.sunposition import sunpos
+import utils.utils as u
+import utils.auxdata as ua
+from process.process import *
+
+plot=False
+odir = os.path.abspath('/DATA/OBS2CO/data/trios/surface_water')
+coordf = glob.glob("/DATA/OBS2CO/data/info/mesures_in_situ.csv")[0]
+coords = pd.read_csv(coordf, sep=';')
+coords['Date_prel']=pd.to_datetime(coords['Date_prel'])
+
+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/test_setup/raw/aw*idpr*.csv")
+swrfiles = glob.glob("/DATA/OBS2CO/data/trios/raw/Lu0*idpr*.csv")
+
+iopw = ua.iopw()
+iopw.load_iopw()
+
+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)
+
+
+idpr = '167'
+
+# get idpr numbers
+idprs = np.unique([re.findall(r'idpr(\d+)', x)[0] for x in swrfiles])
+#idprs = np.array(['170'])
+# loop over idpr
+for idpr in idprs:
+ c = coords[coords.ID_prel == int(idpr)] # .values[0]
+ date=c['Date_prel'].dt.strftime('%Y%m%d')
+ lat = c['Lat'].values[0]
+ lon = c['Lon'].values[0]
+ alt = 0 # c['Altitude']
+ name = c['ID_lac'].values[0]
+
+ ofile=os.path.join(odir,'Rrs_swr_idpr'+idpr+'_'+name+'_'+date.values[0]+'.csv')
+ header = c.stack()
+ header.index = header.index.droplevel()
+ header.to_csv(ofile, header=None)
+
+
+
+ # -----------------------------------------------
+ # 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()
+ Rrs_stat = Rrs_swr.describe()
+ Rrs_stat.columns=Rrs_stat.columns.droplevel()
+ Rrs_stat = Rrs_stat.T
+ Rrs_stat.to_csv(ofile,mode='a')
+ if plot:
+ 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)
+ 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')
diff --git a/setup.py b/setup.py
index 422fe17..e450da8 100644
--- a/setup.py
+++ b/setup.py
@@ -19,7 +19,7 @@ setup(
author_email='tristan.harmel@ntymail.com',
description='Package to help process TriOS radiometer data for various above-water or in-water setups',
# TODO update Dependent packages (distributions)
- install_requires=['pandas', 'scipy', 'numpy', 'netCDF4', 'matplotlib', 'docopt', 'GDAL', 'python-dateutil'],
+ install_requires=['dash','dash_core_components','dash_html_components','pandas', 'scipy', 'numpy', 'netCDF4', 'matplotlib', 'docopt', 'GDAL', 'python-dateutil'],
entry_points={
'console_scripts': [
diff --git a/trios.egg-info/PKG-INFO b/trios.egg-info/PKG-INFO
index de2f375..2dce818 100644
--- a/trios.egg-info/PKG-INFO
+++ b/trios.egg-info/PKG-INFO
@@ -1,6 +1,6 @@
Metadata-Version: 1.0
Name: trios
-Version: 1.1.0
+Version: 1.1.1
Summary: Package to help process 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
diff --git a/trios.egg-info/SOURCES.txt b/trios.egg-info/SOURCES.txt
index 9bd76c5..bd300bd 100644
--- a/trios.egg-info/SOURCES.txt
+++ b/trios.egg-info/SOURCES.txt
@@ -1,9 +1,13 @@
README.md
setup.py
+aux/__init__.py
+aux/sensors_spec.py
process/__init__.py
process/process.py
process/process_compar_awr.py
+process/process_sabine.py
process/process_test_setup.py
+process/run_iwr.py
simulation/__init__.py
simulation/rho_snapshot.py
trios.egg-info/PKG-INFO
@@ -12,8 +16,10 @@ trios.egg-info/dependency_links.txt
trios.egg-info/entry_points.txt
trios.egg-info/requires.txt
trios.egg-info/top_level.txt
+utils/DBtrios.py
utils/DButils.py
utils/__init__.py
+utils/auxdata.py
utils/format_info4acix.py
utils/sunposition.py
utils/utils.py
diff --git a/trios.egg-info/requires.txt b/trios.egg-info/requires.txt
index 25ab4e9..a424c23 100644
--- a/trios.egg-info/requires.txt
+++ b/trios.egg-info/requires.txt
@@ -1,3 +1,6 @@
+dash
+dash_core_components
+dash_html_components
pandas
scipy
numpy
diff --git a/trios.egg-info/top_level.txt b/trios.egg-info/top_level.txt
index 14f0538..60bfdba 100644
--- a/trios.egg-info/top_level.txt
+++ b/trios.egg-info/top_level.txt
@@ -1,3 +1,5 @@
+aux
+build
process
simulation
utils
--
GitLab