run_iwr.py 13.6 KB
Newer Older
1

2
3
4
5
6
7
import os, sys
import pandas as pd
import numpy as np
import glob
import re
import matplotlib as mpl
8
import matplotlib.pyplot as plt
9
10
11
from matplotlib.backends.backend_pdf import PdfPages
from scipy.interpolate import interp1d

12
import cmocean
13
14
import plotly
import plotly.graph_objs as go
15
16
import scipy.optimize as so

17
18
import trios.utils.utils as u
import trios.utils.auxdata as ua
19
from trios.process import *
20

21
22

class fit:
23
24
25
    def __init__(self, N=0, m=2):
        self.popt = np.full([N, m], np.nan)
        self.pcov = np.full([N, m, m], np.nan)
26
27


28
29
30
31
32
33
34
# import aeronet
# from config import *


# ------------------------------------------------
# above-water data files
dirfig = os.path.abspath('/DATA/OBS2CO/data/trios/fig')
35
odir = os.path.abspath('/DATA/OBS2CO/data/trios/in_water')
36

37
iwrfiles = glob.glob("/DATA/OBS2CO/data/trios/raw/2018/uw*idpr*.csv")
38
39
40

coordf = glob.glob("/DATA/OBS2CO/data/info/mesures_in_situ.csv")[0]
coords = pd.read_csv(coordf, sep=';')
41
coords['Date_prel'] = pd.to_datetime(coords['Date_prel'])
42
43
44
# get idpr numbers
idprs = np.unique([re.findall(r'idpr(\d+)', x)[0] for x in iwrfiles])

45

46
def add_curve(ax, x, mean, std=None, c='red', label='',**kwargs):
47
    ax.plot(x, mean, linestyle='solid', c=c, lw=2.5,
48
49
50
51
52
            alpha=0.8, label=label,*kwargs)
    if np.any(std):
        ax.fill_between(x,
                        mean - std,
                        mean + std, alpha=0.35, color=c)
53

54

55
56
57
58
59
60
################
# load aux data
iopw = ua.iopw()
iopw.load_iopw()
irr = ua.irradiance()
irr.load_F0()
61
# TODO check noise values (e.g., NEI from Trios), should it be spectral?
62
noise = 0.1
63
64
65
66
idpr = '169'

idx_list_for_plot=(28, 37, 51, 71, 91, 105, 130, 140, 170)
#idx_list_for_plot=(165,170,175,180,185,190,195,200,205)
67
# loop over idpr
68
for idpr in idprs[0:]:
69
70
71
72
    #    idpr=idprs[2]
    print(idpr)
    try:
        c = coords[coords.ID_prel == int(idpr)]  # .values[0]
73
        date = c['Date_prel'].dt.strftime('%Y%m%d')
74
75
76
77
78
        lat = c['Lat'].values[0]
        lon = c['Lon'].values[0]
        alt = 0  # c['Altitude']
        name = c['ID_lac'].values[0]

79
80
81
82
        # -----------------------------------------------
        #   write data header
        # -----------------------------------------------
        ofile = os.path.join(odir, 'iwr_' + date.values[0] + '_idpr' + idpr + '_' + name + '.csv')
83
84
        header = c.stack()
        header.index = header.index.droplevel()
85
        header.to_csv(ofile, header=False)
86

87
88
89
90
91
92
        dff = pd.DataFrame()

        # -----------------------------------------------
        #   IWR processing
        # -----------------------------------------------
        iwr = u.iwr_data(idpr, iwrfiles)
93
94
95
96
97
        # if iwr.file:
        df, wl_ = iwr.reader(lat, lon,
                             alt)  # depth data already corrected for sensor position  , delta_Lu_depth=0.07,delta_Edz_depth=-0.28)
        # else:
        #    continue
98
99
100
101
102
103

        # ---------------------------
        # Data formatting
        # ---------------------------
        # set pandas dataframes with data and parameters
        df['rounded_depth', ''] = df.prof_Edz.round(1)
104
105
106
107
108

        # data filtering
        df[df.rounded_depth < -0.05] = np.nan
        #df[df == 0] = np.nan

109
110
111
112
113
114
        mean = df.groupby('rounded_depth').mean()
        median = df.groupby('rounded_depth').median()
        std = df.groupby('rounded_depth').std()
        q25 = df.groupby('rounded_depth').quantile(0.25)
        q75 = df.groupby('rounded_depth').quantile(0.75)

115
116
117
118
119
        # ---------------------------
        # Data processing
        # ---------------------------
        # load process object
        process = iwr_process(wl=wl_)
120

121
122
123
124
        # process data
        res = process.process(mean, std)
        res_w = process.process(mean, std, mode='log')
        res_med = process.process(median, std)
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
        res_lsq = process.process(mean, std, mode='lsq')

        Kd = res_lsq.popt[:, 0]
        Edm = res_lsq.popt[:, 1]
        Es = mean.Ed.mean()
        Lwm = res_lsq.popt[:, 3]
        rrs = Lwm / Edm
        Es_from_Edm = 0.96 * Edm
        Lwp = 0.541 * Lwm
        Rrs_Edp = Lwp / Es
        Rrs_Edm = Lwp / Es_from_Edm

        Rrs_df = pd.DataFrame(Rrs_Edp)
        Rrs_df.columns = ['Rrs']
        Rrs_df['Kd']=Kd

        # -----------------------------------------------
        #   write data header
        # -----------------------------------------------
        Rrs_df.to_csv(ofile,mode='a')


        # Open SWR files for comparison
        swr = pd.read_csv(
            '/DATA/OBS2CO/data/trios/surface_water/Rrs_swr_' + date.values[0] + '_idpr' + idpr + '_' + name + '.csv',
            skiprows=57)
151

152
153
154
        # ---------------------------
        # Plotting section
        # ---------------------------
155
156
157
158
159
160
161

        pdf = PdfPages(os.path.join(dirfig, 'trios_iwr_all_idpr' + idpr + '_' + name + '.pdf'))

        # -------------------------------------
        # Profile plotting
        # -------------------------------------
        # Ed profile
162
        i = 0
163
164
        x = mean.prof_Edz
        depth_ = np.linspace(0, x.max(), 200)
165
166
167

        mpl.rcParams.update({'font.size': 12})

168
169
        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)
170
171

        for idx in idx_list_for_plot:
172
173
174
175
176
            ax = axs.flat[i]
            i += 1
            y = mean.Edz.iloc[:, idx]
            sigma = std.Edz.iloc[:, idx]
            yerr = [median.Edz.iloc[:, idx] - q25.Edz.iloc[:, idx], q75.Edz.iloc[:, idx] - median.Edz.iloc[:, idx]]
177
            ax.errorbar(x, y, yerr=yerr, ecolor='lightgray', alpha=0.43, fmt='none', elinewidth=3, capsize=0, label='std')
178
179
180
181
182
            ax.scatter(x, y,
                       # s = std.Edz.iloc[:,50]/std.Edz.iloc[:,50].mean()+100,
                       c=mean.Ed.iloc[:, idx],
                       alpha=0.6, cmap=cmocean.cm.thermal, label=None
                       )
183
            Ed_sim = iwr_process().f_Edz(depth_, *res.popt[idx, :])
184
            ax.plot(depth_, Ed_sim, linestyle='-', c='black', label='mean')
185
            Ed_sim = iwr_process().f_Edz(depth_, *res_w.popt[idx, :])
186
            ax.plot(depth_, Ed_sim, linestyle=':', c='black', label='log-space')
187
188
            # Ed_sim = iwr_process.f_Edz(depth_, *res_rw[idx][0])
            # ax.plot(depth_, Ed_sim, linestyle=':', c='red', label='mean, relativ weighted')
189
            Ed_sim = iwr_process().f_Edz(depth_, *res_med.popt[idx, :])
190
191
            ax.plot(depth_, Ed_sim, linestyle='--', c='black', label='median')

192
            ax.semilogy()
193
194
195
196
            # 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 = ' +
197
198
                         str(round(res_w.popt[idx, 0], 3)) + '$m^{-1}$, Ed0 =' + str(round(res_w.popt[idx, 1], 1)),
                         fontsize=12)
199
200
201
        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)
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
        # fig.savefig(os.path.join(dirfig, 'trios_iwr_Ed_profile_idpr' + idpr + '_' + name + '.pdf'))
        pdf.savefig()
        plt.close()

        # Lu profile
        i = 0
        x = mean.prof_Luz
        depth_ = np.linspace(0, x.max(), 200)
        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 idx_list_for_plot:
            ax = axs.flat[i]
            i += 1
            y = mean.Luz.iloc[:, idx]
            sigma = std.Luz.iloc[:, idx]
            yerr = [y - q25.Luz.iloc[:, idx], q75.Luz.iloc[:, idx] - y]
            yerr = sigma
            ax.errorbar(x, y, yerr=yerr, ecolor='lightgray', alpha=0.43, fmt='none', elinewidth=3, capsize=0, label='std')
            ax.scatter(x, y,
                       # s = std.Luz.iloc[:,50]/std.Luz.iloc[:,50].mean()+100,
                       c=mean.Ed.iloc[:, idx],
                       alpha=0.6, cmap=cmocean.cm.thermal, label=None
                       )
            Lu_sim = iwr_process().f_Lu(depth_, *res_lsq.popt[idx, -2:])
            ax.plot(depth_, Lu_sim, linestyle='-', c='black', label='mean')
            # Lu_sim = iwr_process().f_Luz(depth_, *res_w.popt[idx, :])
            # ax.plot(depth_, Lu_sim, linestyle=':', c='black', label='log-space')
            # # Lu_sim = iwr_process.f_Luz(depth_, *res_rw[idx][0])
            # # ax.plot(depth_, Lu_sim, linestyle=':', c='red', label='mean, relativ weighted')
            # Lu_sim = iwr_process().f_Luz(depth_, *res_med.popt[idx, :])
            # ax.plot(depth_, Lu_sim, linestyle='--', c='black', label='median')

            ax.semilogy()
            # ax.colorbar()
            ax.set_ylabel(r'$L_u\ ({mW\cdot m^{-2}\cdot sr^{-1}\cdot nm^{-1}})$')
            ax.set_xlabel('Depth (m)')
            ax.set_title(r'$\lambda = $' + str(round(wl_[idx], 1)) + r'$ nm, K_{Lu} = $' +
                         str(round(res_lsq.popt[idx, 2], 3)) + '$m^{-1}$, Lu0 =' + str(round(res_lsq.popt[idx, 3], 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)
        # fig.savefig(os.path.join(dirfig, 'trios_iwr_Lu_profile_idpr' + idpr + '_' + name + '.pdf'))
        pdf.savefig()

250
251
        plt.close()

252
253
254
        # -------------------------------------
        # Sprectra plotting
        # -------------------------------------
255
        mpl.rcParams.update({'font.size': 18})
256
        fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(15, 12))
257
        fig.subplots_adjust(left=0.1, right=0.9, hspace=.5, wspace=0.29)
258
259
260

        iparam = 0
        ax = axs.flat[iparam]
261
262
263
264
265
266
267
268
        y, sigma = res.popt[:, iparam], res.pcov[:, iparam, iparam]
        add_curve(ax, wl_, y, sigma, c='red', label='mean')
        y, sigma = res_w.popt[:, iparam], res_w.pcov[:, iparam, iparam]
        add_curve(ax, wl_, y, sigma, c='orange', label='log-space')
        y, sigma = res_med.popt[:, iparam], res_med.pcov[:, iparam, iparam]
        add_curve(ax, wl_, y, sigma, c='green', label='median')
        y, sigma = res_lsq.popt[:, iparam], np.sqrt(res_lsq.pcov[:, iparam, iparam])
        add_curve(ax, wl_, y, sigma, c='blue', label='lsq')
269
270
        ax.set_ylabel(r'$K_{d}\ ({m^{-1}})$')
        ax.set_xlabel(r'Wavelength (nm)')
271
272
273

        iparam = 1
        ax = axs.flat[iparam]
274
275
276
277
278
279
280
281
        y, sigma = res.popt[:, iparam], np.sqrt(res.pcov[:, iparam, iparam])
        add_curve(ax, wl_, y, sigma, c='red', label='mean')
        y, sigma = res_w.popt[:, iparam], np.sqrt(res_w.pcov[:, iparam, iparam])
        add_curve(ax, wl_, y, sigma, c='orange', label='log-space')
        y, sigma = res_med.popt[:, iparam], np.sqrt(res_med.pcov[:, iparam, iparam])
        add_curve(ax, wl_, y, sigma, c='green', label='median')
        y, sigma = res_lsq.popt[:, iparam], np.sqrt(res_lsq.pcov[:, iparam, iparam])
        add_curve(ax, wl_, y, sigma, c='blue', label='lsq')
282
283
284
285
        add_curve(ax, wl_, mean.Ed.mean(), mean.Ed.std(), c='black', label='Es')
        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)')
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303

        iparam = 2
        ax = axs.flat[iparam]
        y, sigma = res_lsq.popt[:, iparam], np.sqrt(res_lsq.pcov[:, iparam, iparam]) * 0
        sigma[sigma > 10 * np.nanmedian(sigma)] = np.nan
        add_curve(ax, wl_, y, sigma, c='blue', label='lsq')
        ax.set_ylabel(r'$K_{Lu}\ ({m^{-1}})$')
        ax.set_xlabel(r'Wavelength (nm)')

        iparam = 3
        ax = axs.flat[iparam]

        y, sigma = res_lsq.popt[:, iparam], np.sqrt(res_lsq.pcov[:, iparam, iparam]) * 0
        sigma[sigma > 10 * np.nanmedian(sigma)] = np.nan
        add_curve(ax, wl_, y, sigma, c='blue', label='lsq')
        ax.set_ylabel(r'$L_{w}(0^{-})\ ({mW\cdot m^{-2}\cdot sr^{-1}\cdot nm^{-1})$')
        ax.set_xlabel(r'Wavelength (nm)')

304
        fig.suptitle('trios_iwr ' + name + ' idpr' + idpr, fontsize=16)
305
306
307
        # fig.savefig(os.path.join(dirfig, 'trios_iwr_l2_idpr' + idpr + '_' + name + '.pdf'))
        pdf.savefig()

308
        plt.close()
309

310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
        # Edm = res_lsq.popt[:,1]
        # Es = mean.Ed.mean()
        # Lwm = res_lsq.popt[:,3]
        # rrs = Lwm / Edm
        # Es_from_Edm = 0.96 * Edm
        # Lwp = 0.541 * Lwm
        # Rrs_Edp = Lwp / Es
        # Rrs_Edm = Lwp / Es_from_Edm

        mpl.rcParams.update({'font.size': 18})
        fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(15, 12))
        fig.subplots_adjust(left=0.1, right=0.9, hspace=.5, wspace=0.29)

        ax = axs.flat[0]
        add_curve(ax, wl_, Es / Edm, c='red', label='Es over Ed0-')
        ax.set_xlabel(r'Wavelength (nm)')
        ax.legend(loc='best', frameon=False)
        ax = axs.flat[1]
        add_curve(ax, wl_, rrs, c='red', label='rrs')
        ax.set_xlabel(r'Wavelength (nm)')
        ax.legend(loc='best', frameon=False)
        ax = axs.flat[2]
        add_curve(ax, wl_, Rrs_Edm, c='red', label='Rrs (from Ed0-)')
        add_curve(ax, wl_, Rrs_Edp, c='blue', label='Rrs (from Ed0+)')
        add_curve(ax, swr.wl, swr['mean'], swr['std'], c='black', label='Rrs (from SWR)')
        ax.set_xlabel(r'Wavelength (nm)')
        ax.legend(loc='best', frameon=False)
        ax = axs.flat[3]
        Rrs_swr = interp1d(swr.wl, swr['mean'], fill_value='extrapolate')(wl_)
        diff = Rrs_Edp - Rrs_swr
        add_curve(ax, wl_, diff, c='black', label='Rrs(IWR) - Rrs(SWR)')

        fig.suptitle('trios_iwr ' + name + ' idpr' + idpr, fontsize=16)
        pdf.savefig()
        plt.close()
        pdf.close()
346
347
348
349

    except:
        plt.close()
        continue