run_iwr.py 7.83 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
9
10
11
import matplotlib.pyplot as plt
import cmocean


12
13
14
15
16
17
import plotly
import plotly.graph_objs as go

import utils.utils as u
import utils.auxdata as ua

18
from trios.process import *
19

20
21
22
23
24
25
26
27


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)

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
38
39
40

iwrfiles = glob.glob("/DATA/OBS2CO/data/trios/raw/uw*idpr*.csv")

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
47
48
49
50
51
52

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)

53
54
55
56
57
58
################
# load aux data
iopw = ua.iopw()
iopw.load_iopw()
irr = ua.irradiance()
irr.load_F0()
59
60
#TODO check noise values (e.g., NEI from Trios), should it be spectral?
noise = 0.1
61
62
63
64
65
66
67
idpr = '178'
# loop over idpr
for idpr in idprs:#[-1:]:
    #    idpr=idprs[2]
    print(idpr)
    try:
        c = coords[coords.ID_prel == int(idpr)]  # .values[0]
68
        date=c['Date_prel'].dt.strftime('%Y%m%d')
69
70
71
72
73
        lat = c['Lat'].values[0]
        lon = c['Lon'].values[0]
        alt = 0  # c['Altitude']
        name = c['ID_lac'].values[0]

74
75
76
77
78
        ofile=os.path.join(odir,'Rrs_swr_'+date.values[0]+'_idpr'+idpr+'_'+name+'.csv')
        header = c.stack()
        header.index = header.index.droplevel()
        header.to_csv(ofile, header=None)

79
80
81
82
83
84
85
86
87
88
89
        dff = pd.DataFrame()

        # -----------------------------------------------
        #   IWR processing
        # -----------------------------------------------
        iwr = u.iwr_data(idpr, iwrfiles)
        if iwr.file:
            df, wl_ = iwr.reader(lat, lon, alt)
            reflectance = iwr_process(df, wl_).process()
            df = pd.concat([df, reflectance], axis=1)

90
            #df.to_csv(os.path.join(odir, 'trios_iwr_' + name + '_idpr' + idpr + '.csv'))
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110

        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)


        N = len(wl_)
        x = mean.prof_Edz #- 0.56
        depth_ = np.linspace(0, x.max(), 200)
        res=fit(N)
        res_w,res_rw,res_med=fit(N),fit(N),fit(N)

        for idx, wl in enumerate(wl_[:-10]):
            aw, bbw = iopw.get_iopw(wl)
            F0 = irr.get_F0(wl)

            y = mean.Edz.iloc[:, idx]
            sigma = std.Edz.iloc[:, idx]
111
            sigma[sigma<noise]=noise
112
113
            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]))
114
            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
115
116
117
118
119
            #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
120
121
122

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

123
124
125
126
127
128
129
130
131
132
133
134
135
136
        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):
            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]]
            ax.errorbar(x, y, yerr=yerr, ecolor='lightgray', fmt='none', elinewidth=3, capsize=0)
            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
                       )
137
            Ed_sim = iwr_process.f_Edz(depth_, *res.popt[idx,:])
138
            ax.plot(depth_, Ed_sim, linestyle='-', c='black', label='mean')
139
140
            Ed_sim = iwr_process.f_Edz(depth_, *res_w.popt[idx,:])
            ax.plot(depth_, Ed_sim, linestyle=':', c='black', label='log-space')
141
142
            #Ed_sim = iwr_process.f_Edz(depth_, *res_rw[idx][0])
            #ax.plot(depth_, Ed_sim, linestyle=':', c='red', label='mean, relativ weighted')
143
            Ed_sim = iwr_process.f_Edz(depth_, *res_med.popt[idx,:])
144
145
            ax.plot(depth_, Ed_sim, linestyle='--', c='black', label='median')

146
            ax.semilogy()
147
148
149
150
            # 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 = ' +
151
                         str(round(res_w.popt[idx,0], 3)) + '$m^{-1}$, Ed0 =' + str(round(res_w.popt[idx,1], 1)), fontsize=12)
152
153
154
155
156
157
158
        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_iw_idpr' + idpr + '_' + name + '.pdf'))
        plt.close()


159
160
161
162

        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)
163
164
165
166
167
168

        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]
169
        add_curve(ax, wl_, y, std, c='orange', label='log-space')
170
        y, std = res_med.popt[:,iparam], res_med.pcov[:,iparam,iparam]
171
172
173
        add_curve(ax, wl_, y, std, c='green', label='median')
        ax.set_ylabel(r'$K_{d}\ ({m^{-1}})$')
        ax.set_xlabel(r'Wavelength (nm)')
174
175
176

        iparam = 1
        ax = axs.flat[iparam]
177
        y, std = res.popt[:,iparam], np.sqrt(res.pcov[:,iparam,iparam])
178
        add_curve(ax, wl_, y, std, c='red', label='mean')
179
180
181
182
183
        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')
184

185
186
187
188
189
190
        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()
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206

        # fig, ax = plt.subplots()
        # N=df.Edz.shape[1]
        # for wl, group in df.Edz.iloc[:,range(0,N,10)].iteritems():
        #
        #     ax.scatter(df.prof_Edz, group, marker='o',  c=cmocean.cm.thermal(wl/850*255))
        #
        # ax.set_ylim((0, 500000))

        #
        # writing output file
        # dff.to_csv(os.path.join(dirout, 'trios_awr_' + name + '_idpr' + idpr + '.csv'))

    except:
        plt.close()
        continue