Newer
Older
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
306
307
308
pro metis_l1_prep
; start the log
journal,'output/metis_l1_prep_log.txt'
; some definitions
metis_datatype = list( $
'vl-image', $ ; 0
'uv-image', $ ; 1
'pcu-accumul', $ ; 2
'vl-temp-matrix', $ ; 3
'uv-temp-matrix', $ ; 4
'vl-c-ray-mat', $ ; 5
'uv-c-ray-mat', $ ; 6
'pcu-event-list', $ ; 7
'pcu-test-event-list', $ ; 8
'light-curve' $ ; 9
)
metis_obj_size = list( $
2048L, $ ; 0
1024L, $ ; 1
1024L, $ ; 2
2048L, $ ; 3
1024L, $ ; 4
2048L, $ ; 5
1024L $ ; 6
)
; read the auxiliary file - here we have all the inputs we need
input = json_parse('input/contents.json', /toarray, /tostruct)
; prepare the spice kernels
kernels = [ $
input.tls_file_name, $
input.tsc_file_name $
]
; load the spice kernels
cspice_furnsh, kernels
; read l0 fits structure
fits_info, input.l0_file_name, n_ext = n_ext, extname = extname, /silent
; read the primary hdu
ext_no = 0
image = mrdfits(input.l0_file_name, ext_no, primary_header, /silent)
; if the data is not an image, read the data binary-table extension
if fxpar(primary_header, 'NAXIS') eq 0 then begin
ext_no += 1
data_bin_table = mrdfits(input.l0_file_name, ext_no, data_extension_header, /silent)
endif else data_bin_table = !null
; read the metadata extension
ext_no += 1
metadata_bin_table = mrdfits(input.l0_file_name, ext_no, metadata_extension_header, /silent)
; identify the data product and its nominal size
datatype = fxpar(metadata_extension_header, 'DATATYPE')
if datatype le 6 then begin
naxis = 2
naxis1 = metis_obj_size[datatype]
naxis2 = naxis1
bitpix = 16
data_size = naxis1 * naxis2
data_volume = (data_size * bitpix) / 8.
endif else begin
naxis = 0
bitpix = 8
data_size = 0
data_volume = 0
endelse
; if the data product is an accumulation matrix, look for the accumulation vector extension and read it if it exists
if datatype eq 2 then begin
if fxpar(metadata_extension_header, 'M') eq 256 then begin
ext_no += 1
accumul_vector = mrdfits(input.l0_file_name, ext_no, vector_extension_header, /silent)
endif else accumul_vector = !null
endif
; NOTE - the radialization extension is ignored
; read the planning data
planning_data = json_parse(input.planning_file_name, /toarray, /tostruct)
; definitions for the primary header
; version of the fits file
version = string(input.l1_version + 1, format = '(I02)')
; creation and acquisition times in utc
date = date_conv(systime(/julian, /utc), 'FITS') + 'Z'
obt_beg = fxpar(primary_header, 'OBT_BEG')
obt_end = fxpar(primary_header, 'OBT_END')
obt_avg = (obt_beg + obt_end) / 2.0D
date_beg = solo_obt2utc(decode_obt(obt_beg, /from_double))
date_end = solo_obt2utc(decode_obt(obt_end, /from_double))
date_avg = solo_obt2utc(decode_obt(obt_avg, /from_double))
date_beg_string = strreplace(strreplace(strmid(date_beg, 0, 19), '-', ''), ':', '')
; name of the fits file
file_name_fields = strsplit(fxpar(primary_header, 'FILENAME'), '_', /extract)
file_name = 'solo_l1_' + file_name_fields[2] + '_' + date_beg_string + '_v' + version + '.fits'
out_file_name = 'output/' + file_name
; instrument keywords
; TODO - complete with filter information
if datatype eq 0 or datatype eq 3 or datatype eq 5 or datatype eq 9 then begin
filter = 'VL'
wavelnth = 0.0
wavemin = 0.0
wavemax = 0.0
waveband = 'Visible light 580-640 nm'
endif else begin
filter = 'UV'
wavelnth = 1215.67
wavemin = 0.0
wavemax = 0.0
waveband = 'H I Lyman-alpha 121.6 nm'
endelse
detector = filter + 'DA'
telescope = 'SOLO/Metis/' + detector
; campaign keywords
obs_id_fields = strsplit(planning_data.obs_id, '_', /extract)
soop_type = obs_id_fields[2]
obs_type = obs_id_fields[4]
; build the fits file extensions
; join the metadata extension header to the primary header removing unwanted keywords
i = where(strmid(primary_header, 0, 8) eq 'COMP_RAT')
j = where(strmid(metadata_extension_header, 0, 8) eq 'EXTNAME ')
n = n_elements(metadata_extension_header)
primary_header = [ $
primary_header[0 : i], $
metadata_extension_header[j + 1 : n - 2], $
primary_header[i + 1 : *] $
]
; adjust the primary header (it is almost the same for all data product types)
fxaddpar, primary_header, 'FILENAME', file_name, 'FITS filename'
fxaddpar, primary_header, 'PARENT', file_basename(input.l0_file_name), 'Name of the parent file that got processed to the current one', before = 'APID'
fxaddpar, primary_header, 'DATE', date, 'Date and time of FITS file creation', before = 'OBT_BEG'
fxaddpar, primary_header, 'DATE-OBS', date_beg, 'Same as DATE-BEG', before = 'OBT_BEG'
fxaddpar, primary_header, 'DATE-BEG', date_beg, 'Start time of observation', before = 'OBT_BEG'
fxaddpar, primary_header, 'DATE-AVG', date_avg, 'Average time of observation', before = 'OBT_BEG'
fxaddpar, primary_header, 'DATE-END', date_end, 'End time of observation', before = 'OBT_BEG'
fxaddpar, primary_header, 'TIMESYS', 'UTC', 'System used for time keywords', before = 'OBT_BEG'
fxaddpar, primary_header, 'TIMRDER', 0.0, 'Estimated random error in time values', before = 'OBT_BEG'
fxaddpar, primary_header, 'TIMSYER', 0.0, 'Estimated systematic error in time values', before = 'OBT_BEG'
fxaddpar, primary_header, 'LEVEL', 'L1', 'Data processing level'
fxaddpar, primary_header, 'CREATOR', 'metis_l1_prep.pro', 'FITS creation software'
fxaddpar, primary_header, 'OBSRVTRY', 'Solar Orbiter', 'Satellite name', before = 'INSTRUME'
fxaddpar, primary_header, 'TELESCOP', telescope, 'Telescope that took the measurement', before = 'INSTRUME'
fxaddpar, primary_header, 'DETECTOR', detector, 'Subunit/sensor', before = 'DATAMIN'
fxaddpar, primary_header, 'OBJECT', 'TBD', 'The use of the keyword OBJECT is [TBD]', before = 'DATAMIN'
fxaddpar, primary_header, 'OBS_MODE', planning_data.obs_mode, 'Observation mode or study that has been used to acquire this image', before = 'DATAMIN'
fxaddpar, primary_header, 'OBS_TYPE', obs_type, 'Subfield of OBS_ID that only contains an encoded version of OBS_MODE', before = 'DATAMIN'
fxaddpar, primary_header, 'FILTER', filter, 'Filter used to acquire this image', before = 'DATAMIN'
fxaddpar, primary_header, 'WAVELNTH', wavelnth, 'Characteristic wavelength at which the observation was taken', before = 'DATAMIN'
fxaddpar, primary_header, 'WAVEMIN', wavemin, 'The shortest wavelength at which the net (approximate) response function becomes 0.05 times the maximum response. ', before = 'DATAMIN'
fxaddpar, primary_header, 'WAVEMAX', wavemax, 'The longest wavelength at which the net (approximate) response function becomes 0.05 times the maximum response', before = 'DATAMIN'
fxaddpar, primary_header, 'WAVEBAND', waveband, 'Description of the wavelength band', before = 'DATAMIN'
fxaddpar, primary_header, 'XPOSURE', 0, 'Total effective exposure time of the observation, in seconds', before = 'DATAMIN'
fxaddpar, primary_header, 'NSUMEXP', 0, 'Number of images summed together to form the observation', before = 'DATAMIN'
fxaddpar, primary_header, 'TELAPSE', 0, 'Total elapsed time between the beginning and end of the complete observation in seconds, including any dead times between exposures', before = 'DATAMIN'
fxaddpar, primary_header, 'SOOPNAME', planning_data.soop_name, 'SOOP(s) that this observation belongs to', before = 'DATAMIN'
fxaddpar, primary_header, 'SOOPTYPE', soop_type, before = 'DATAMIN'
fxaddpar, primary_header, 'OBS_ID', planning_data.obs_id, 'Unique identifier for the observation that is associated with the data acquisition', before = 'DATAMIN'
fxaddpar, primary_header, 'TARGET', 'TBD', 'Taget as defined in the SOOP description', before = 'DATAMIN'
fxaddpar, primary_header, 'BSCALE', 1.0, before = 'DATAMIN'
fxaddpar, primary_header, 'BZERO', 0, before = 'DATAMIN'
fxaddpar, primary_header, 'BTYPE', strupcase(metis_datatype[datatype]), 'Science data object type', before = 'DATAMIN'
fxaddpar, primary_header, 'BUNIT', 'DN', 'Units of physical value, after application of BSCALE and BZERO', before = 'DATAMIN'
if datatype le 6 then begin
fxaddpar, primary_header, 'PXBEG1', 1, 'First pixel that has been read out in dimension 1', before = 'COMPRESS'
fxaddpar, primary_header, 'PXBEG2', 1, 'First pixel that has been read out in dimension 2', before = 'COMPRESS'
fxaddpar, primary_header, 'PXEND1', naxis1, 'Last pixel that has been read out in dimension 1', before = 'COMPRESS'
fxaddpar, primary_header, 'PXEND2', naxis2, 'Last pixel that has been read out in dimension 2', before = 'COMPRESS'
endif
fxaddpar, primary_header, 'NBIN1', 1, 'Binning factor in the dimension 1', before = 'COMPRESS'
fxaddpar, primary_header, 'NBIN2', 1, 'Binning factor in the dimension 2', before = 'COMPRESS'
fxaddpar, primary_header, 'NBIN', 1, 'Product of all NBIN values above', before = 'COMPRESS'
fxaddpar, primary_header, 'IDB_VERS', input.idb_version, '', before = 'HDR_VERS'
fxaddpar, primary_header, 'INFO_URL', 'http://metis.oato.inaf.it', 'Link to more information on the instrument data', before = 'HISTORY'
; add checksum and datasum to the fits header
; WARN - should this be done at the end? i don't know
fits_add_checksum, primary_header, image
; add keywords for file history
fxaddpar, primary_header, 'HISTORY', 'Third history line'
; remove unused keywords
sxdelpar, primary_header, 'CAD_BEG'
sxdelpar, primary_header, 'CAD_END'
sxdelpar, primary_header, 'WIDTH'
sxdelpar, primary_header, 'HEIGHT'
sxdelpar, primary_header, 'NB_IMG'
sxdelpar, primary_header, 'N'
sxdelpar, primary_header, 'X_SIZE'
sxdelpar, primary_header, 'Y_SIZE'
sxdelpar, primary_header, 'Z_SIZE'
sxdelpar, primary_header, 'P_BANDS'
sxdelpar, primary_header, 'N_BANDS'
sxdelpar, primary_header, 'ORIG_X'
sxdelpar, primary_header, 'ORIG_Y'
sxdelpar, primary_header, 'FIRSTROW'
mwrfits, image, out_file_name, primary_header, /no_comment, /create, /silent
; if applicable, save the data binary-table extension as it is
if isa(data_bin_table) then mwrfits, data_bin_table, 'output/' + file_name, data_extension_header, /no_comment, /silent
; read the house-keeping telemetry
hk_bin_table = json_parse(input.hk_file_name, /toarray, /tostruct)
; replace raw values with calibrated values in the primary header
; WARN - must be done here?
if datatype eq 0 or datatype eq 3 or datatype eq 5 then begin
fxaddpar, primary_header, 'DAC1POL1', interpol_param(hk_bin_table, 'NIT0E061', obt_avg)
fxaddpar, primary_header, 'DAC2POL1', interpol_param(hk_bin_table, 'NIT0E062', obt_avg)
fxaddpar, primary_header, 'DAC1POL2', interpol_param(hk_bin_table, 'NIT0E064', obt_avg)
fxaddpar, primary_header, 'DAC2POL2', interpol_param(hk_bin_table, 'NIT0E065', obt_avg)
fxaddpar, primary_header, 'DAC1POL3', interpol_param(hk_bin_table, 'NIT0E067', obt_avg)
fxaddpar, primary_header, 'DAC2POL3', interpol_param(hk_bin_table, 'NIT0E068', obt_avg)
fxaddpar, primary_header, 'DAC1POL4', interpol_param(hk_bin_table, 'NIT0E06A', obt_avg)
fxaddpar, primary_header, 'DAC2POL4', interpol_param(hk_bin_table, 'NIT0E06B', obt_avg)
fxaddpar, primary_header, 'TSENSOR', interpol_param(hk_bin_table, 'NIT0E0E0', obt_avg)
fxaddpar, primary_header, 'PMPTEMP', interpol_param(hk_bin_table, 'NIT0L00D', obt_avg)
endif
if datatype eq 1 or datatype eq 4 or datatype eq 6 then begin
fxaddpar, primary_header, 'HVU_SCR', interpol_param(hk_bin_table, 'NIT0E070', obt_avg)
fxaddpar, primary_header, 'HVU_MCP', interpol_param(hk_bin_table, 'NIT0E071', obt_avg)
fxaddpar, primary_header, 'TSENSOR', interpol_param(hk_bin_table, 'NIT0E050', obt_avg)
endif
; build the telemetry extension
hk_extension_header = !null
fxaddpar, hk_extension_header, 'XTENSION', 'BINTABLE', 'Extension type'
fxaddpar, hk_extension_header, 'BITPIX', 8, 'Number of bits per data pixel'
fxaddpar, hk_extension_header, 'NAXIS', 2, 'Number of data axes'
fxaddpar, hk_extension_header, 'NAXIS1', 0, 'Number of pixels along axis 1'
fxaddpar, hk_extension_header, 'NAXIS2', 0, 'Number of pixels along axis 2'
fxaddpar, hk_extension_header, 'PCOUNT', 0, 'Parameter count'
fxaddpar, hk_extension_header, 'GCOUNT', 1, 'Group count'
fxaddpar, hk_extension_header, 'EXTNAME', 'House-keeping', 'Extension name'
mwrfits, hk_bin_table, out_file_name, hk_extension_header, /no_comment, /silent
; write the auxiliary information file
output = { $
file_name: out_file_name, $
l0_file_name: input.l0_file_name, $
log_file_name: '' $
}
openw, unit, 'output/contents.json'
printf, unit, output, /implied_print
free_lun, unit
; json_write, output, 'output/contents.json'
; unload the spice kernels
cspice_unload, kernels
; close the log
journal
end