process_compar_awr.py 9.72 KB
Newer Older
1
2
import glob
import re
3
4

import matplotlib as mpl
5
import matplotlib.pyplot as plt
6
plt.rcParams.update({'font.size': 18})
7
8
from scipy.interpolate import interp1d

9
10
from trios.utils.sunposition import sunpos
from trios.utils import utils as u
11
from trios.process import *
12
13
14
15
16

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')

17
awrfiles = glob.glob("/DATA/OBS2CO/data/trios/raw/2018/aw*idpr*.csv")
18
# awrfiles = glob.glob("/DATA/OBS2CO/data/trios/test_setup/raw/aw*idpr*.csv")
19
swrfiles = glob.glob("/DATA/OBS2CO/data/trios/raw/2018/Lu0*idpr*.csv")
20
21
22
23

iopw = ua.iopw()
iopw.load_iopw()

24

25
def add_curve(ax, x, mean, std=None, c='red', label='',**kwargs):
26
    ax.plot(x, mean, linestyle='solid', c=c, lw=2.5,
27
28
29
30
31
            alpha=0.8, label=label,*kwargs)
    if np.any(std):
        ax.fill_between(x,
                        mean - std,
                        mean + std, alpha=0.35, color=c)
32
33
34
35
36
37


idpr = '167'

# get idpr numbers
idprs = np.unique([re.findall(r'idpr(\d+)', x)[0] for x in swrfiles])
38
# idprs = np.array(['170'])
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
# loop over idpr
for idpr in idprs:
    c = coords[coords.ID_prel == int(idpr)]  # .values[0]
    lat = c['Lat'].values[0]
    lon = c['Lon'].values[0]
    alt = 0  # c['Altitude']
    name = c['ID_lac'].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()

    # -----------------------------------------------
    #   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
104
105
106
107
        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")
108
109
110

        # add solar angle data and idpr
        # compute solar angle (mean between fisrt and last aqcuisition time
111
112
        raw['sza', ''] = np.nan
        for index, row in raw.iterrows():
113
114
            # print index
            sza = sunpos(index, lat, lon, alt)[1]
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
            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
        # ------------------
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
        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^{+})$')

        # ---- Lsky
        ax2 = ax.twinx()
        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]
        add_curve(ax, wl, Lt.mean(axis=0), Lt.std(axis=0),
                  label=r'$L_t$',c='black')
        add_curve(ax, wl, Lsky.mean(axis=0)*rho15, Lsky.std(axis=0)*rho15,
                  label='M2015 (' + str(round(rho15,4)) + ')',c='violet')
        add_curve(ax, wl, Lsky.mean(axis=0)*rho99, Lsky.std(axis=0)*rho99, c='orange',
                  label='M1999(' + str(round(rho99,4)) + ')')
        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]


        add_curve(ax, wl, Lsky.mean(axis=0)*rho15/Ltm, Lsky.std(axis=0)*rho15,
                  label='M2015 (' + str(round(rho15,4)) + ')',c='violet')
        add_curve(ax, wl, Lsky.mean(axis=0)*rho99/Ltm, Lsky.std(axis=0)*rho99, c='orange',
                  label='M1999(' + str(round(rho99,4)) + ')')
        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)')
189

190
191
192
193
194
195
196
197
198
199
200
201
        # ---- Lw
        ax = axs[1,0]
        add_curve(ax, wl, Rrs15.mean(axis=0)*Edm, Rrs15.std(axis=0)*Edm,
                  label='M2015 (' + str(round(rho15,4)) + ')',c='violet')
        add_curve(ax, wl, Rrs99.mean(axis=0)*Edm, Rrs99.std(axis=0)*Edm, c='orange',
                  label='M1999(' + str(round(rho99,4)) + ')')
        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)) + ')')
        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)')
202

203
204
205
        # ---- Rrs
        ax = axs[1,1]
        add_curve(ax, wl_swr, Rrs_swr.transpose().mean(axis=1), Rrs_swr.transpose().std(axis=1), label='swr', c='black')
206
        add_curve(ax, wl, Rrs15.transpose().mean(axis=1), Rrs15.transpose().std(axis=1),
207
                  label='M2015 (' + str(round(rho15,4)) + ')',c='violet')
208
        add_curve(ax, wl, Rrs99.transpose().mean(axis=1), Rrs99.transpose().std(axis=1), c='orange',
209
210
211
212
213
                  label='M1999(' + str(round(rho99,4)) + ')')
        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)) + ')')
        add_curve(ax, wl, Rrs_opt, Rrs_opt_std, c='blue',
                  label='Optimization')
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
        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
        add_curve(ax, wl, (Rrs15.mean(axis=0)-Rrs_swr_)/Rrs_swr_,
                  label='M2015 (' + str(round(rho15,4)) + ')',c='violet')
        add_curve(ax, wl, (Rrs99.mean(axis=0)-Rrs_swr_)/Rrs_swr_, c='orange',
                  label='M1999(' + str(round(rho99,4)) + ')')
        add_curve(ax, wl, (Rrs_h.mean(axis=0)-Rrs_swr_)/Rrs_swr_, c='green',
                  label='h(' + str(round(rho_h.mean(),4)) + ')')
        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)
233

234
235
236
        fig.suptitle('trios_awr ' + name + ' idpr' + idpr, fontsize=16)
        fig.savefig(os.path.join(dirfig, 'trios_awr_' + name + '_idpr' + idpr + '.png'))
        plt.close()