This repository has been archived by the owner on Nov 10, 2020. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Hyperspectral_Calibration_20191212.txt
305 lines (252 loc) · 14 KB
/
Hyperspectral_Calibration_20191212.txt
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
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
# -*- coding: utf-8 -*-
"""
*** Radiometric calibration of hyperspectral imagery ***
This framework will process the raw hyperspectral VNIR and SWIR imagery using the pre-computed
calibration models and real-time downwelling irradiance data, convert raw image digital number(DN)
to reflectance and save as netCDF files.
Input: 1) VNIR and SWIR raw hyperspectral imagery
2) Spectral data from downwelling irradiance sensor
3) Pre-computed calibration models (includes four models: vnir_new, vnir_middle, vnir_old, swir_new; models are not avaialbe for swir_old and swir_middle
because no downwelling irridiance data/sensor was avaiable for that time period)
Output: Hyperspectral reflectance imagery with netCDF format
@author: Remote Sensing Lab at Saint Louis University
"""
import json
import numpy as np
import os
import spectral.io.envi as envi
from netCDF4 import Dataset
raw_root = "/home/extractor/sites/ua-mac/raw_data"
# raw_root = "/home/extractor/hs_calib"
# extract spectral profiles from environmentlogger.json
def irradiance_time_extractor(camera_type,envlog_file):
# For the environmental logger records after 04/26/2016, there would be 24 files per day (1 file per hour, 5 seconds per record)
# Convert json fiel to dictionary format file
with open(envlog_file, "r") as fp:
lines = fp.readlines()
slines = "".join(lines)
js = json.loads(slines)
# assume that time stamp follows in 5 second increments across records since 5 sec/record
num_readings = len(js["environment_sensor_readings"])
if "spectrometers" in js["environment_sensor_readings"][0]:
if camera_type == "swir_new":
num_bands = len(js["environment_sensor_readings"][0]["spectrometers"]["NIRQuest-512"]["spectrum"])
else:
num_bands = len(js["environment_sensor_readings"][0]["spectrometers"]["FLAME-T"]["spectrum"])
else:
num_bands = len(js["environment_sensor_readings"][0]["spectrometer"]["spectrum"])
spectra = np.zeros((num_readings, num_bands))
times = []
for idx in range(num_readings):
# read time stamp
time_current = js["environment_sensor_readings"][idx]["timestamp"]
C = time_current.replace("."," ").replace("-"," ").replace(":","")
ArrayTime=C.split(" ")
time_current_r = int(ArrayTime[3])
times.append(time_current_r)
# read spectrum from irridiance sensors
if "spectrometers" in js["environment_sensor_readings"][idx]:
if camera_type == "swir_new":
spectrum = js["environment_sensor_readings"][0]["spectrometers"]["NIRQuest-512"]["spectrum"]
else:
spectrum = js["environment_sensor_readings"][idx]["spectrometers"]["FLAME-T"]["spectrum"]
else:
spectrum = js["environment_sensor_readings"][idx]["spectrometer"]["spectrum"]
spectra[idx,:] = spectrum
return times, spectra
# replace rfl_img variable in netcdf with given matrix
def update_netcdf(inp, rfl_data, camera_type):
print("Updating %s" % inp)
out = inp.replace(".nc", "_newrfl.nc")
with Dataset(inp) as src, Dataset(out, "w") as dst:
# copy global attributes all at once via dictionary
dst.setncatts(src.__dict__)
# copy dimensions
for name, dimension in src.dimensions.items():
dst.createDimension(name, (len(dimension) if not dimension.isunlimited() else None))
# copy all file data except for the excluded
for name, variable in src.variables.items():
if name == "Google_Map_View":
continue
# Create variables
var_dict = (src[name].__dict__)
if '_FillValue' in var_dict.keys():
x = dst.createVariable(name, variable.datatype, variable.dimensions, fill_value=var_dict['_FillValue'])
del var_dict['_FillValue']
else:
x = dst.createVariable(name, variable.datatype, variable.dimensions)
# Set variables to values
if name != "rfl_img":
print("...%s" % name)
dst[name][:] = src[name][:]
else:
if camera_type=='vnir_old':
print("...%s (subset)" % name)
dst[name][:679,:,:] = rfl_data
# 679-955 set to NaN
print("...NaNs")
dst[name][679:,:,:] = np.nan
elif camera_type == "vnir_middle":
print("...%s (subset)" % name)
dst[name][:662,:,:] = rfl_data
# 679-955 set to NaN
print("...NaNs")
dst[name][662:,:,:] = np.nan
else:
print("...%s" % name)
dst[name][:] = rfl_data
# copy variable attributes all at once via dictionary
dst[name].setncatts(var_dict)
if 'rfl_img' not in src.variables:
print("...adding rfl_img")
dst.createVariable("rfl_img", "f4")
dst.variables['rfl_img'] = rfl_data
# apply calibration algorithm to the raw data
def apply_calibration(raw_filepath):
print("Calibrating %s" % raw_filepath)
# get necessary paths from path to _raw file
raw_dir = os.path.dirname(raw_filepath)
raw_file = os.path.basename(raw_filepath)
md_file = os.path.join(raw_dir, "%s_metadata.json" % raw_file[:-4])
date = raw_filepath.split("/")[-3]
timestamp = raw_filepath.split("/")[-2]
envlog_dir = os.path.join(raw_root, "EnvironmentLogger/%s" % date)
# determine type of sensor and age of camera
if raw_filepath.find("VNIR") > -1:
if date < "2018-08-18":
camera_type = "vnir_old"
num_spectral_bands = 955
num_bands_irradiance = 1024
image_scanning_time = 540
elif "2018-08-18" <= date < "2019-02-26":
camera_type = "vnir_middle"
num_spectral_bands = 939
num_bands_irradiance = 1024
image_scanning_time = 540
else:
camera_type = "vnir_new"
num_spectral_bands = 939
num_bands_irradiance = 3648
# it is obverved that it takes an average of 3.5 mins/scan = 210 seconds
image_scanning_time = 210
else:
if date < "2019-02-26": # Note that no calibration models are available for old&middle swir data
camera_type = "swir_old_middle"
else:
camera_type = "swir_new"
num_spectral_bands = 275
num_bands_irradiance = 512
image_scanning_time = 210
print("MODE: ---------- %s ----------" % camera_type)
# load the raw data set
print("Loading %s.hdr" % raw_filepath)
try:
raw = envi.open(raw_filepath +'.hdr')
# img_DN = raw.load()
img_DN = raw.open_memmap()
#head_file = envi.read_envi_header(data_fullpath +'.hdr')
except IOError:
print('No such file named %s' % raw_filepath)
# Apply calibration procedure if camera_type == vnir_old, vnir_middle, vnir_new or swir_new. Since no calibration models are available for
# swir_old and swir_middle, so directly convert old&middel SWIR raw data to netcdf format
if camera_type =="swir_old_middle":
# Convert the raw swir_old and swir_middle data to netCDF
img_DN = np.rollaxis(img_DN, 2, 0)
# Generate output path and call the netCDF conversion function, convert the raw old&middle swir data to netcdf
out_path = os.path.dirname(raw_filepath.replace("raw_data", "Level_1").replace("SWIR", "swir_netcdf").replace("VNIR", "vnir_netcdf"))
out_file = os.path.join(out_path, "%s_netcdf_L1_ua-mac_%s.nc" % (camera_type.split("_")[0], timestamp))
update_netcdf(out_file, img_DN, camera_type)
# free up memory
del img_DN
else: # when camera_type == vnir_old, vnir_middle, vnir_new or swir_new, apply pre-computed calibration models
# Load the previously created calibration models based on the camera_type
best_matched = os.path.join(raw_root + "/" + "calibration_new", camera_type, 'best_matched_index.npy')
bias = os.path.join(raw_root + "/" + "calibration_new", camera_type, 'bias_coeff.npy')
gain = os.path.join(raw_root + "/" + "calibration_new", camera_type, 'gain_coeff.npy')
# read EnvLog data
print("Reading EnvLog files in %s" % envlog_dir)
envlog_tot_time = []
envlog_spectra = np.array([], dtype=np.int64).reshape(0, num_bands_irradiance)
for ef in os.listdir(envlog_dir):
if ef.endswith("environmentlogger.json"):
print(ef)
time, spectrum = irradiance_time_extractor(camera_type,os.path.join(envlog_dir, ef))
envlog_tot_time += time
# print("concatenating %s onto %s" % (spectrum.shape, envlog_spectra.shape))
envlog_spectra = np.vstack([envlog_spectra, spectrum])
# Find the best match time range between image time stamp and EnvLog time stamp
num_irridiance_record = int(image_scanning_time/5) # 210/5=4.2 ----> 5 seconds per record
# concatenation of hour mins and seconds of the image time stamp (eg., 12-38-49 to 123849)
with open(md_file) as json_file:
img_meta_data = json.load(json_file)
meta_data_time = img_meta_data['lemnatec_measurement_metadata']['gantry_system_variable_metadata']['time']
image_time = meta_data_time[-8:]
image_time = int(image_time.replace(":",""))
# compute the absolute difference between
print("Computing mean spectrum")
abs_diff_time = np.zeros((len(envlog_tot_time)))
for k in range(len(envlog_tot_time)):
abs_diff_time[k] = abs(image_time - envlog_tot_time[k])
ind_closet_time = np.argmin(abs_diff_time) # closest time index
mean_spectrum = np.mean(envlog_spectra[ind_closet_time : ind_closet_time + num_irridiance_record-1, :], axis=0)
# load pre-computed the best matched index between image and irradiance sensor spectral bands
best_matched_index = np.load(best_matched)
test_irridance = mean_spectrum[best_matched_index.astype(int).tolist()]
test_irridance_re = np.resize(test_irridance, (1, num_spectral_bands))
# load and apply precomputed coefficient to convert irradiance to DN
b = np.load(bias)
g = np.load(gain)
if camera_type == "vnir_old":
test_irridance_re = test_irridance_re[:,0:679]
img_DN = img_DN[:,:,0:679]
if camera_type == "vnir_middle":
test_irridance_re = test_irridance_re[:,0:662]
img_DN = img_DN[:,:,0:662]
irrad2DN = (g * test_irridance_re) + b
# reflectance computation
print("Computing reflectance")
rfl_data = img_DN/irrad2DN
rfl_data = np.rollaxis(rfl_data, 2, 0)
# free up memory
del img_DN
del irrad2DN
# prepare output paths
print("Generating output")
out_path = os.path.dirname(raw_filepath.replace("raw_data", "Level_1").replace("SWIR", "swir_netcdf").replace("VNIR", "vnir_netcdf"))
if not os.path.isdir(out_path):
os.makedirs(out_path)
# save as ENVI file (RGB bands: 392, 252, 127)
#out_file = os.path.join('ref_%s.hdr' % raw_file)
#envi.save_image(out_file, Ref, dtype=np.float32, interleave='bil', force = 'True', metadata=head_file)
# Save Ref as a .npy file
#out_file = os.path.join(out_path, 'ref_%s.npy' % raw_file)
#np.save(out_file, rfl_data)
# Write to nc file
out_file = os.path.join(out_path, "%s_netcdf_L1_ua-mac_%s.nc" % (camera_type.split("_")[0], timestamp))
update_netcdf(out_file, rfl_data, camera_type)
# free up memory
del rfl_data
# TODO: This will come from the extractor message
input_paths = [
# swir_old
# NA os.path.join(raw_root, "SWIR/2017-04-16/2017-04-16__11-50-46-707/c6079666-b686-4481-9a4f-0663f5f43a6a_raw"),
# swir_new
# OK os.path.join(raw_root, "SWIR/2018-09-22/2018-09-22__13-21-35-977/05b7ad1a-a2d7-4dfc-bcec-b1a394ec0892_raw"),
# OK os.path.join(raw_root, "SWIR/2018-10-11/2018-10-11__12-11-43-420/dc60c7d5-24bc-432b-a7cb-98ac1da73154_raw"),
# vnir_old
# OK os.path.join(raw_root, "VNIR/2017-04-15/2017-04-15__11-33-42-265/76efd15f-928a-49f7-a008-4877b8842129_raw"),
# OK os.path.join(raw_root, "VNIR/2017-04-17/2017-04-17__16-39-35-738/a5096c7a-c052-4728-a29b-a5a6119c366c_raw"),
# OK os.path.join(raw_root, "VNIR/2017-04-17/2017-07-08__06-30-15-622/5154c9fc-0a51-4a4d-9c79-242539f057ad_raw"),
# os.path.join(raw_root, "VNIR/2017-05-13/2017-05-13__12-00-39-756/1bcc7cd0-1205-45a3-b4b3-cbcac6236754_raw"), # Killed
# os.path.join(raw_root, "VNIR/2017-06-18/2017-06-18__14-34-24-390/41a0b327-83ff-4131-b1fd-5ee5254760b6_raw"), # Killed
# os.path.join(raw_root, "VNIR/2017-07-27/2017-07-27__15-05-11-667/d1643679-bef3-4179-9912-d63bf4cd53c6_raw"),
#os.path.join(raw_root, "VNIR/2017-08-23/2017-08-23__09-21-43-959/ea0e3408-ed1c-412d-aa68-e75ce2e902b1_raw"),
#os.path.join(raw_root, "VNIR/2018-09-26/2018-09-26__13-36-22-843/73e4642f-1ac7-430a-a00e-a04224bde9db_raw"),
# vnir_middle
# NA os.path.join(raw_root, "VNIR/2018-08-18/2018-08-18__11-11-41-890/c5f4d50f-44ad-4e23-9d92-f10e62110ac7_raw"),
# NA os.path.join(raw_root, "VNIR/2018-10-08/2018-10-08__11-41-01-365/5e39a30f-d343-405e-a140-db26dc72eb59_raw"),
# vnir_new
# GEN os.path.join(raw_root, "VNIR/2019-06-17/2019-06-17__14-03-29-760/d51b6f4c-9246-4da8-9a6c-786bb1dc21bf_raw"),
]
for p in input_paths:
apply_calibration(p)