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

implementation of filtering and temporal optimization scheme for above-water radiometry

parent 6078bfa1
include aux/*
\ No newline at end of file
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 as mpl
import plotly
import plotly.graph_objs as go
import utils.utils as u
from trios.utils import utils as u
from trios.process import *
# import aeronet
......
......@@ -5,12 +5,12 @@ setup(
name='trios',
version=__version__,
packages=find_packages(exclude=['build']),
package_data={'': ['*.so']},
package_data={'': ['*.so'],
# # If any package contains *.txt files, include them:
# '': ['*.txt'],
# 'lut': ['data/lut/*.nc'],
# 'aux': ['data/aux/*']
# },
# 'aux': ['aux/*']
},
include_package_data=True,
url='https://gitlab.irstea.fr/ETL-TELQUEL/etl/tree/dev/preprocessing/trios',
......
This diff is collapsed.
Metadata-Version: 1.0
Name: trios
Version: 1.1.1
Version: 1.1.3
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
......
MANIFEST.in
README.md
setup.py
aux/Thuillier_2003_0.3nm.dat
aux/__init__.py
aux/mobley2datasheet.sh
aux/rhoTable_Mobley1999.csv
aux/rhoTable_Mobley1999.txt
aux/rhoTable_Mobley2015.csv
aux/rhoTable_Mobley2015.txt
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
simulation/__init__.py
simulation/rho_snapshot.py
trios/__init__.py
......@@ -19,13 +29,13 @@ 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
trios/utils/DBtrios.py
trios/utils/DButils.py
trios/utils/__init__.py
trios/utils/auxdata.py
trios/utils/format_info4acix.py
trios/utils/sunposition.py
trios/utils/utils.py
visu/__init__.py
visu/data_visu.py
visu/visu.py
\ No newline at end of file
import pandas as pd
from scipy import interpolate, integrate
import scipy.optimize as so
# import plotly.graph_objs as go
from plotly.graph_objs import *
import plotly.graph_objs as go
from utils.utils import reshape as r
import utils.auxdata as ua
from trios.utils.utils import reshape as r
import trios.utils.auxdata as ua
from trios.config import *
......@@ -81,18 +81,71 @@ class awr_process:
rho_mobley = calc().spline_4d(grid, rho_, (ws, sza, vza, azi))
return rho_mobley
def call_process(self, vza=[40], azi=[135], ws=2, aot=0.1):
wl = self.wl
Lt = self.df.loc[:, ("Lt")]
Lsky = self.df.loc[:, ("Lsky")]
Ed = self.df.loc[:, ("Ed")]
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 call_process(self, vza=[40], azi=[135], ws=2, aot=0.1):
# wl = self.wl
# Lt = self.df.loc[:, ("Lt")]
# Lsky = self.df.loc[:, ("Lsky")]
# Ed = self.df.loc[:, ("Ed")]
# 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]):
# '''
#
# :param wl:
# :param Lt:
# :param Lsky:
# :param Ed:
# :param sza:
# :param vza:
# :param azi:
# :param ws:
# :param aot:
# :return:
# '''
#
# # ------------------
# # filtering
# # ------------------
# ind_Ed, notused = calc.spectra_median_filter(Ed)
# ind_sky, notused = calc.spectra_median_filter(Lsky)
# ind = ind_Ed & ind_sky
# Lt, Lsky, Ed, sza = Lt[ind], Lsky[ind], Ed[ind], sza[ind]
#
# rho = self.get_rho_values(np.median(sza), vza, azi, wl=wl, ws=ws, aot=aot)
# self.Rrs = (Lt - rho * Lsky.values) / Ed.values
# self.Rrs.columns = pd.MultiIndex.from_product([['Rrs(awr)'], wl], names=['param', 'wl'])
# return self.Rrs, rho
@staticmethod
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(self, wl, df, sza, vza=[40], azi=[135], ws=[2], aot=[0.1], method='M99'):
print(sza, vza, azi, ws, aot, method)
Rrs, rho = self.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(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], method='M99'):
'''
:param wl:
......@@ -104,13 +157,74 @@ class awr_process:
:param azi:
:param ws:
:param aot:
:param method: 'M99, 'M15, 'osoaa'
:return:
'''
rho = self.get_rho_values([sza], vza, azi, wl=wl, ws=ws, aot=aot)
# -----------------------------
# standard data processing
# -----------------------------
if method == 'osoaa':
rho = self.get_rho_values(np.median(sza), vza, azi, wl=wl, ws=ws, aot=aot)
elif method == 'M99':
rho = self.get_rho_mobley(self.rhoM1999, [np.median(sza)], [vza], [azi], [ws])
elif method == 'M15':
rho = self.get_rho_mobley(self.rhoM2015, [np.median(sza)], [vza], [azi], [ws])
else:
return print('ERROR: no method for rho factor')
self.Rrs = (Lt - rho * Lsky) / Ed
# Rrs.columns = pd.MultiIndex.from_product([['Rrs'], wl], names=['param', 'wl'])
return self.Rrs, rho.mean()
def cost_func(self, x, param, meas, Rrs_bar):
sza, vza, azi = param
Lt, Lsky, Ed = meas
ws = x
rho = self.get_rho_mobley(self.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(self, wl, Lt, Lsky, Ed, sza, vza=[40], azi=[135], ws=[2], aot=[0.1], method='M99'):
return self.Rrs, rho
# ------------------------------
# initialization of mean/median values
# ------------------------------
Rrs, rho = self.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 = []
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(self.cost_func, x0, 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)
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
class swr_process:
......@@ -209,24 +323,25 @@ class iwr_process:
df['rounded_depth', ''] = df.prof_Edz.round(1)
df.groupby('rounded_depth').mean()
# TODO add central part of run_iwr.py here
return self.reflectance
@staticmethod
def f_Edz(depth, Kd, Ed0):
'''simple Edz model for homogeneous water column'''
return Ed0 * np.exp(-Kd*depth)
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
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,x='Luz',y='prof_Luz'):
trace = Scattergl(
def plot_raw(self, x='Luz', y='prof_Luz'):
trace = go.Scattergl(
x=self.df[x].values,
y=self.df[y].values,
......@@ -249,7 +364,7 @@ class iwr_process:
)
layout = Layout(
layout = go.Layout(
height=450,
xaxis=dict(
range=[0, 200],
......@@ -269,14 +384,14 @@ class iwr_process:
zeroline=False,
# nticks=max(6, round(df['Speed'].iloc[-1]/10))
),
margin=Margin(
margin=go.Margin(
t=45,
l=50,
r=50
)
)
return Figure(data=[trace], layout=layout)
return go.Figure(data=[trace], layout=layout)
class self_shading:
......@@ -322,7 +437,7 @@ class calc:
wl = wl[idx_par]
Ed = Ed[idx_par]
par = integrate.trapz(Ed, wl)
par_quanta = integrate.trapz(np.multiply(wl,Ed), wl) / hc
par_quanta = integrate.trapz(np.multiply(wl, Ed), wl) / hc
return par, par_quanta
def earth_sun_correction(self, dayofyear):
......@@ -352,6 +467,19 @@ class calc:
'''
return np.degrees(np.arcsin(np.sin(np.radians(angle_air)) / n))
@staticmethod
def spectra_median_filter(spectra, threshold=0.1):
'''
:param series: pandas object
:param threshold: relative value of median
:return: boolean indices, array of data within interval median +/- threshold
'''
spec = spectra.sum(axis=1)
med = spec.median()
ind = np.abs(1 - spec / med) < 0.1
return ind, spectra[ind]
def spline_2d(self, gin, arr, gout):
'''
Interpolation of a 6D array (arr) with bicubic splines on a 2D grid
......@@ -397,15 +525,18 @@ class calc:
for i in range(N[0]):
for j in range(N[1]):
tmp[i, j, :] = interpolate.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] = interpolate.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] = interpolate.RectBivariateSpline(gin[0], gin[1], tmp[:, :, iband])(gout[0], gout[1],
grid=True)
interp[:, :, iband] = interpolate.RectBivariateSpline(gin[0], gin[1], tmp[:, :, iband])(gout[0],
gout[1],
grid=True)
return interp
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 trios.utils.sunposition import sunpos
from trios.utils import utils as u
from trios.process import *
coordf = glob.glob("/DATA/OBS2CO/data/info/mesures_in_situ.csv")[0]
......
......@@ -3,20 +3,19 @@ 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
import utils.utils as u
from trios.process import *
from utils.sunposition import sunpos
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')
......@@ -26,7 +25,7 @@ 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]'
groupname = dff.columns.values[1] # 'RT$station_area[n]'
def get_param(df, param='Ld', name='Lsky'):
......@@ -41,13 +40,18 @@ for group, df_ in dff.groupby(groupname):
print(group)
ofile = os.path.join(odir, 'awr_' + group + '.csv')
Rrs_df = pd.DataFrame()
Rrs_ = 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]
ws = [df.WIND_SPEED.values.mean() / 2]
# TODO check altitude
alt = 0
......@@ -60,23 +64,186 @@ for group, df_ in dff.groupby(groupname):
# 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.values, Ed.values, sza, ws=ws, azi=azi)
# 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.columns = pd.MultiIndex.from_product([['Rrs(awr)'], wl], names=['param', 'wl'])
# ['Rrs(awr)'] * Rrs.shape[1]
Rrs = pd.concat([df.iloc[:, 0:10], Ed,Lsky,Lt,Rrs], axis=1)
Rrs['sza']=sza
Rrs['azi']=azi[0]
Rrs['rho']=rho.mean()
Rrs_ = pd.concat([Rrs_, Rrs], axis=0)
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_