process_test_setup.py 4.55 KB
Newer Older
1
2
3
4
import glob
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

5
6
from trios.utils.sunposition import sunpos
from trios.utils import utils as u
7
from trios.process import *
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126


coordf = glob.glob("/DATA/OBS2CO/data/info/mesures_in_situ.csv")[0]
coords = pd.read_csv(coordf, sep=';')
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")

idpr='167'

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
# -----------------------------------------------

swr = u.swr_data(idpr, swrfiles)
if swr.file:
    df, wl = swr.reader(lat, lon, alt)
    Rrs_swr = swr_process(df, wl).process()

# -----------------------------------------------
#   AWR processing
# -----------------------------------------------
awr = u.awr_data(idpr, awrfiles)

index_idx=[2,0,1]

d=u.data(index_idx)
Ed, wl_Ed = d.load_csv(awr.Edf)
Lsky, wl_Lsky = d.load_csv(awr.Lskyf)
Lt0, wl_Lt = d.load_csv(awr.Ltf)

# ''' interpolate Ed and Lsky data upon Lt wavelength'''
wl = wl_Lt
Lt0.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]
fig, axs = plt.subplots(nrows=3, ncols=4, figsize=(16, 10))
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]
        i=i+1
        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
        df = pd.merge_asof(Lt, Ed, left_index=True, right_index=True, tolerance=pd.Timedelta("2 seconds"),
                           direction="nearest")
        df = pd.merge_asof(df, Lsky, left_index=True, right_index=True, tolerance=pd.Timedelta("2 seconds"),
                           direction="nearest")

        # add solar angle data and idpr
        # compute solar angle (mean between fisrt and last aqcuisition time
        df['sza', ''] = np.nan
        for index, row in df.iterrows():
            # print index
            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])

        Rrs_h =(df.loc[:,'Lt'] -rho_h*df.loc[:,'Lsky'])/ df.loc[:,'Ed']
        Rrs15 = (df.loc[:,'Lt'] -rho15*df.loc[:,'Lsky'])/ df.loc[:,'Ed']

        Rrs99 = (df.loc[:,'Lt'] -rho99*df.loc[:,'Lsky'])/ df.loc[:,'Ed']
        #plt.figure()


        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)
        add_curve(ax,wl,Rrs_swr.transpose().mean(axis=1),Rrs_swr.transpose().std(axis=1),label='swr',c='black')
        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')

        ax.set_title('azi='+str(azi)+', vza='+str(vza))


        ax.legend(loc='best', frameon=False)

        ax.set_ylabel(r'$R_{rs}\  (sr^{-1})$')
        ax.set_xlabel(r'Wavelength (nm)')

Lt.index.names