/
RISA_WC_Experiment_DRIVER.py
186 lines (134 loc) · 8 KB
/
RISA_WC_Experiment_DRIVER.py
1
2
3
4
5
6
7
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
127
128
129
130
131
132
133
134
135
136
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
"""
This is the driver script for running weathering crust radiative transfer experiments
for the RISA project.
This requires imports from the BiOSNICAR_GO_PY repository and must therefore be run in
the BiOSNICAR_GO_PY directory. This can be downloaded from www.github.com/jmcook1186/BiOSNICAR_GO_PY/
This script along with RISA_WC_Experiment_FUNCS.py must then be copied into the
top level directory and run in the BioSNICAR_py Python 3 environment.
This driver script calls functions stored in RISA_WC_Experiment_FUNCS.py.
There are three functions which themselves call local functions or functions from BioSNICAR_GO_PY.
These three are:
find_best_params(): This is used to do multiple runs of SNICAR and find the parameter set that
best simulates a given field-measured spectrum
match_field_spectrum(): for a given set of input params this function plots measured vs simulated spectra
build_LUT(): constructs look up tables o be used by inverse model
inverse_model(): Uses lookup tables to predict cell concentrations from spectra
isolate_biological_effect(): calculate the proportion of albedo change attributed to biological or physical
processes by spectral differencing.
compare_predicted_and_measured(): provide a comparison between predictions made by the inverse model,
2BDA model and measured cell concentrations. Saves a figure to savepath.
All calls to SNICAR assume solar zenith = 45 degrees, underlying surface has BBA = 0.1, Picard (2016) ice refractive
index is used, grains are spherical.
"""
import numpy as np
import dask
import pandas as pd
from RISA_WC_Experiment_FUNCS import *
# FIRST determine whether we work in albedo or reflectance -
# apply ARF for reflectance, do not apply ARF for albedo
apply_ARF = False
plot_ARF = False
# define path to save figures
savepath = '/home/tothepoles/Desktop/'
# provide list of samples to include as examples of each surface type
CIsites = ['21_7_S4', '13_7_SB3', '14_7_SB6', '15_7_S4', '15_7_SB1',
'15_7_SB5', '21_7_S2', '21_7_SB3', '22_7_S2',
'22_7_S4', '23_7_SB1', '23_7_SB2', '23_7_S4', 'WI_3',
'WI_8', 'WI_10', 'WI_11']
DISPsites = ['DISP1','DISP2','DISP3','DISP4','DISP5','DISP6',
'DISP7','DISP8', 'DISP9','DISP10','DISP11','DISP12',
'DISP13','DISP14', '27_7_16_DISP1','27_7_16_DISP3']
LAsites = ['13_7_SB1','13_7_SB3','14_7_SB2','14_7_SB3','22_7_SB3']
HAsites = ['14_7_SB9','13_7_SB2','14_7_SB1','14_7_SB5','14_7_SB9',
'14_7_SB10','15_7_SB3','20_7_SB1','21_7_SB2','21_7_SB4',
'22_7_SB2','22_7_SB5','23_7_SB4', '23_7_SB5']
# generate ARF data
ARF_CI_DF = pd.read_csv('/datadrive2/BigIceSurfClassifier/Training_Data/ARF_master.csv')
ARF_CI = np.mean(np.array(ARF_CI_DF[ARF_CI_DF.columns.intersection(CIsites)]),axis=1)
ARF_CI = ARF_CI[::10]
ARF_HA_DF = pd.read_csv('/datadrive2/BigIceSurfClassifier/Training_Data/ARF_master.csv')
ARF_HA = np.mean(np.array(ARF_HA_DF[ARF_HA_DF.columns.intersection(HAsites)]),axis=1)
ARF_HA = ARF_HA[::10]
# set path to field spectral database - reflectance if ARF is applied, albedo if not
# The ARF adjusts the SNICAR simulated albedo to give nadir reflectance
if apply_ARF:
field_data_fname = '/datadrive2/BigIceSurfClassifier/Training_Data/HCRF_master_16171819.csv'
else:
field_data_fname = '/datadrive2/BigIceSurfClassifier/Training_Data/Albedo_master.csv'
path_to_metadata = '/datadrive2/BigIceSurfClassifier/Spectra_Metadata.csv'
# set parallel arrays of filenames plus density, radius, thickness, model algal concn (ppb)
# and measured algal concn (cells/mL) for each sample to plot.
# These value are likely generated by repeated calls to find_best_params()
if apply_ARF:
fnames= ['23_7_SB1','14_7_SB6','14_7_SB9','14_7_SB1','22_7_SB5','14_7_SB2', '22_7_SB3', 'RAIN2']
rho = [[650,650],[650,650],[750,750],[800,800],[700,700],[800,800],[750,750],[900,900]]
rds = [[650,650],[650,650],[800,800],[800,800],[650,650],[750,750],[750,750],[900,900]]
dz = [[0.001,0.05],[0.001,0.08],[0.001,0.04],[0.001,0.03],[0.001,0.025],[0.001,0.06],[0.001,0.04],[0.001,0.04]]
alg = [[0,0],[0,0],[14000,0],[18000,0],[14000,0],[1500,0],[4000,0],[0,0]] #cells in ppb
measured_cells=[0,0,21875,30312,22813,3062,6687,0] #cells in cells/mL
else:
fnames= ['23_7_SB1','14_7_SB6','14_7_SB9','14_7_SB1','22_7_SB5','14_7_SB2', '22_7_SB3', 'RAIN2']
rho = [[600,600],[650,650],[800,800],[800,800],[700,700],[800,800],[700,700],[900,900]]
rds = [[600,600],[650,650],[800,800],[800,800],[700,700],[800,800],[700,700],[900,900]]
dz = [[0.001,0.05],[0.001,0.08],[0.001,0.04],[0.001,0.03],[0.001,0.02],[0.001,0.06],[0.001,0.04],[0.001,0.04]]
alg = [[0,0],[0,0],[14000,0],[18000,0],[14000,0],[1500,0],[4000,0],[0,0]] #cells in ppb
measured_cells=[0,0,21875,30312,22813,3062,6687,0] #cells in cells/mL
#################
# CALL FUNCTIONS
#################
# 1) if aiming to find best param set for simulating a specific field spectrum,
# provide sample ID or one of "CImean", "LAmean" or "HAmean" to find_best_params().
# If no LAPs, set clean=True
#namelist=['13_7_SB1','13_7_SB2','13_7_SB3','13_7_SB5','14_7_SB1','14_7_SB2'],\
# '14_7_SB3','14_7_SB5','14_7_SB6','14_7_SB8','14_7_SB9','14_7_SB10',\
# '15_7_SB1','15_7_SB2','15_7_SB3','15_7_SB4','15_7_SB5','17_7_SB1',\
# '17_7_SB2','20_7_SB1','20_7_SB2','20_7_SB3','21_7_SB2','21_7_SB3',\
# '21_7_SB4','21_7_SB9','22_7_SB2','22_7_SB3','22_7_SB5','22_7_SB6',\
# '22_7_SB7','23_7_SB1','23_7_SB2','23_7_SB4','23_7_SB5','24_7_SB1']
# ResultArray = np.zeros(shape=(len(namelist),6))
# for i in np.arange(0,len(namelist),1):
# fn = namelist[i]
# Results = find_best_params(field_data_fname, fn, CIsites, LAsites, HAsites, weight = 1, clean=False)
# Results = np.array(Results)
# best_idx = Results[:,1].argmin()
# best_params = Results[best_idx,0]
# best_error = Results[best_idx,1]
# best_dens, best_rds, best_dz, best_alg, best_zen = zip(best_params)
# ResultArray[i,:] = np.array(best_dens),np.array(best_rds),np.array(best_dz),np.array(best_alg),np.array(best_zen),np.array(best_error)
# import pandas as pd
# Out = pd.DataFrame(data=ResultArray,columns=['dens','rds','dz','alg','zen','spec_err'])
# Out.index=namelist
# Out.to_csv('/home/tothepoles/Desktop/Out.csv')
# 2) once best param set has been determined using find_best_params for 8 spectra and
# they are collatd into parallel arrays, pass them to this function to display
# them in a multipanel figure
# match_field_spectra(field_data_fname, fnames, rho, rds, dz, alg, measured_cells,\
# CIsites, LAsites, HAsites, apply_ARF, plot_ARF, ARF_CI, ARF_HA, savepath)
# 3) to estimate biological and physical contributions to albedo reduction
# delAbioLA, delAphysLA, delAbioHA, delAphysHA = isolate_biological_effect(field_data_fname,CIsites,LAsites,HAsites, savepath)
# print("**LIGHT ALGAE**")
# print("Total BBA change = {}".format(delAbioLA+delAphysLA))
# print("Change in BBA due to glacier algae = {}".format(delAbioLA))
# print("Change in BBA due to weathering crust degradation = {}".format(delAphysLA))
# print("**HEAVY ALGAE**")
# print("Total BBA change = {}".format(delAbioHA+delAphysHA))
# print("Change in BBA due to glacier algae = {}".format(delAbioHA))
# print("Change in BBA due to weathering crust degradation = {}".format(delAphysHA))
# # BUILD LUT
# dz = [0.02, 0.03, 0.04, 0.05, 0.075 , 0.1]
# densities = [600, 650, 700, 750, 800, 850]
# radii = [600, 700, 800, 900, 1000]
# algae = [0, 5000, 7500, 10000, 12500, 15000, 17500, 20000, 25000, 30000]
# solzens = [40, 45, 50, 55]
# wavelengths = np.arange(0.2,5,0.01)
# save_LUT = True
# savepath = '/home/tothepoles/Desktop/'
# LUT = build_LUT(solzens, dz, densities, radii, algae, wavelengths, save_LUT,\
# apply_ARF, ARF_CI, ARF_HA, savepath)
# 5) run inverse model
path_to_LUTs = '/home/tothepoles/Desktop/'
field_data_fname = '/datadrive2/BigIceSurfClassifier/Training_Data/HCRF_master_16171819.csv'
Out = inverse_model(field_data_fname,path_to_LUTs)
Out.to_csv('/home/tothepoles/Desktop/Out.csv')
# 6) Field vs model comparison
#compare_predicted_and_measured(savepath)