forked from bishun945/pyOWT
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathOpticalVariables.py
183 lines (132 loc) · 6.56 KB
/
OpticalVariables.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
from pandas import read_csv
import numpy as np
import os
import yaml
class OpticalVariables():
def __init__(self, Rrs, band, sensor=None):
if not isinstance(Rrs, np.ndarray):
raise TypeError("Input 'Rrs' should be np.ndarray type.")
# manipulate dimension of input Rrs as we always assume it is 3d nparray (raster-like)
# the wavelength should be on the shape[2] dim
if np.ndim(Rrs) == 1:
# here assume a signle spectrum
self.Rrs = Rrs.reshape(1, 1, Rrs.shape[0])
elif np.ndim(Rrs) == 2:
# here assume shape[0] is sample and shape[1] is wavelength
self.Rrs = Rrs.reshape(Rrs.shape[0], 1, Rrs.shape[1])
elif np.ndim(Rrs) == 3:
# here assume shape[2] is wavelength
self.Rrs = Rrs
else:
raise ValueError("Input 'Rrs' should only have 1, 2, or 3 dims!")
self.band = np.array(band)
self.sensor = sensor
self.AVW = None
self.Area = None
self.NDI = None
# TODO: check the input band fits the selected sensor range
# dont_TODO: if AVW ends by 700 nm, this list has to be modified
with open('data/sensor_band_library.yaml', 'r') as file:
sensor_lib = yaml.load(file, Loader=yaml.FullLoader)
self.sensor_AVW_bands_library = sensor_lib['lib_800']['sensor_AVW_bands_library']
self.sensor_RGB_bands_library = sensor_lib['lib_800']['sensor_RGB_bands_library']
# dont_TODO: if AVW ends by 700 nm, this list has to be modified
self.sensor_RGB_min_max = sensor_lib['lib_800']['sensor_RGB_min_max']
self.AVW_regression_coef = sensor_lib['lib_800']['AVW_regression_coef']
self.available_sensors = ', '.join(self.sensor_AVW_bands_library.keys())
if self.sensor is None:
# the input Rrs are assumed to be hyperspectral
self.spectral_attr = "hyper"
ref_sensor_RGB_bands = [443, 560, 665]
self.sensor_RGB_bands = [self.band[np.argmin(abs(self.band - v))] for v in ref_sensor_RGB_bands]
self.sensor_band_min = 400
self.sensor_band_max = 800
self.AVW_conver_coef = [0, 1, 0, 0, 0, 0]
else:
if self.sensor not in self.sensor_AVW_bands_library.keys():
raise ValueError(f"The input `sensor` couldn't be found in the library: {self.available_sensors}")
# if `sensor` specifized, trigger `conver_AVW_multi_to_hyper`
self.spectral_attr = "multi"
# define RGB bands, by B, R, G order
ref_sensor_RGB_bands = self.sensor_RGB_bands_library[self.sensor]
self.sensor_RGB_bands = [self.band[np.argmin(abs(self.band - v))] for v in ref_sensor_RGB_bands]
# define min max ranges
self.sensor_band_min = self.sensor_RGB_min_max[self.sensor]["min"]
self.sensor_band_max = self.sensor_RGB_min_max[self.sensor]["max"]
# read coefficients to convert AVW_multi to AVW_hyper
proj_root = os.path.dirname(os.path.abspath(__file__))
fn = os.path.join(proj_root, self.AVW_regression_coef)
d = read_csv(fn)
self.AVW_convert_coef = d[d["variable"] == sensor][["0", "1", "2", "3", "4", "5"]].values.tolist()[0]
# run calculation
self.calculate_AVW()
self.calculate_Area()
self.calculate_NDI()
class ArrayWithAttributes:
'''This subclass add attributes to np.array
Outputs of `OpticalVariables` (AVW, Area, NDI) can have attributes indicating
some key features during the calculation. For example, we should know which
bands were used to calcualte AVW, Area, and NDI.
Examples:
A = ArrayWithAttributes(np.array([1, 2, 3]), author='Shun')
print(A) # Output: [1 2 3]
print(A.author) # Output: Shun
'''
def __init__(self, array, **attributes):
self.array = np.asarray(array)
self.__dict__.update(attributes)
def __getitem__(self, item):
return self.array[item]
def __setitem__(self, key, value):
self.array[key] = value
def __repr__(self):
return repr(self.array)
def __str__(self):
return str(self.array)
def __getattr__(self, name):
return getattr(self.array, name)
def __array__(self):
return self.array
def convert_AVW_multi_to_hyper(self):
self.AVW_hyper = np.zeros(self.AVW_multi.shape)
for i in range(len(self.AVW_convert_coef)):
self.AVW_hyper += self.AVW_convert_coef[i] * (self.AVW_multi**i)
def calculate_AVW(self):
idx_for_AVW = (self.band >= self.sensor_band_min) & (self.band <= self.sensor_band_max)
bands_for_AVW = self.band[idx_for_AVW]
Rrs_for_AVW = self.Rrs[:, :, idx_for_AVW]
self.AVW_init = np.sum(Rrs_for_AVW, axis=-1) / np.sum(Rrs_for_AVW / bands_for_AVW[None, None, :], axis=-1)
if self.spectral_attr == "hyper":
self.AVW_hyper = self.AVW_init
else:
self.AVW_multi = self.AVW_init
self.convert_AVW_multi_to_hyper()
self.AVW = self.AVW_hyper
def calculate_Area(self):
bands_for_Area = np.array(self.sensor_RGB_bands)
Rrs_for_Area = self.Rrs[:, :, np.where(np.isin(self.band, bands_for_Area))[0]]
self.Area = np.trapz(x=bands_for_Area, y=Rrs_for_Area, axis=-1)
def calculate_NDI(self):
Rrs_for_NDI = self.Rrs[:, :, np.where(np.isin(self.band, self.sensor_RGB_bands[1:]))[0]] # [1:] for G and R
self.NDI = -np.diff(Rrs_for_NDI, axis=-1)[:, :, 0] / np.sum(Rrs_for_NDI, axis=-1)
def run(self):
# TODO: deprecate this func in the future
import warnings
warnings.warn(
"No need to calculate via 'ov.run()'. "
"The optical variables will be calculated directly once you create this instance. "
"This function is deprecated and will be removed in the future versions.",
DeprecationWarning,
stacklevel=2
)
self.calculate_AVW()
self.calculate_Area()
self.calculate_NDI()
if __name__ == "__main__":
# ov = OpticalVariables(Rrs=1, band=1, sensor="OLCI_S3B")
# print(ov.sensor_RGB_bands)
band = np.arange(400, 801, step = 2)
# Rrs = np.random.normal(loc = 0, scale = 1, size = len(band))
Rrs = np.full(len(band), 1)
ov = OpticalVariables(Rrs=Rrs, band=band)
print(ov.AVW)