Select Git revision
metis_l2_prep_vl_generic.pro
metis_l2_prep_vl_generic.pro 12.07 KiB
pro metis_l2_prep_vl_generic
; keyword defining if the detector reference frame must be used for the output
ref_detector = 0
; start the log
journal, 'output/metis_l2_prep_log.txt'
; read the auxiliary file - here we have all the inputs we need
input = json_parse('input/contents.json', /toarray, /tostruct)
; load the spice kernels
load_spice_kernels, input.spice_kernels, kernel_list = kernel_list, kernel_version = kernel_version
journal, 'SPICE kernel files correctly loaded:'
journal, ' SDK version= ' + kernel_version
; read the calibration package
cal_pack = json_parse(input.cal_pack_path + '/index.json', /toarray, /tostruct)
; include the calibration package path into the cal_pack structure
cal_pack = create_struct('path', input.cal_pack_path + path_sep(), cal_pack)
journal, 'Calibration package correctly imported:'
journal, ' version = ' + cal_pack.version
journal, ' validity range = ' + string(cal_pack.validity_range.start, cal_pack.validity_range._end, format = '(A, "-", A)')
; read the primary hdu
data = mrdfits(input.file_name, 0, primary_header, /silent)
; patch to fix the lack of the keyword SEQ_NUM
pol_id = fxpar(primary_header, 'POL_ID', missing = 0)
if pol_id ge 1 and pol_id le 4 then begin
seq_num = fxpar(primary_header, 'SEQ_NUM', missing = 0)
if seq_num eq 0 then begin
obj_cnt = fxpar(primary_header, 'OBJ_CNT')
n_pol = fxpar(primary_header, 'N_POL', missing = 4)
seq_num = ((obj_cnt - 1) / n_pol + 1)
fxaddpar, primary_header, 'SEQ_NUM', seq_num, before = 'POL_ID'
endif
endif
; read the quality matrix
quality_matrix = mrdfits(input.file_name, 'quality matrix', /silent)
; convert the string header into an idl structure
header = fits_hdr2struct(primary_header)
journal, 'L1 FITS file correctly read:'
journal, ' filename = ' + file_basename(input.file_name)
journal, ' datatype = ' + string(header.datatype, format = '(I0)')
journal, ' sess_num = ' + header.sess_num
journal, ' obj_cnt = ' + string(header.obj_cnt, format = '(I0)')
journal, ' pol_id = ' + string(header.pol_id, format = '(I0)')
journal, ' nbin = ' + string(sqrt(header.nbin), format = '(I0)')
; calibration block
file_name_fields = strsplit(header.filename, '_', /extract)
; error is the quadratic relative error
error = 0.
case header.datatype of
0: begin
data = double(data)
; ====================================
; for old l1 data
; data = data * header.nbin * header.ndit
; header.xposure = header.xposure * header.ndit
; header.nsumexp = header.ndit
; ====================================
data = metis_dark_vlda(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
; ====================================
; for data already subtracted of dark
; file = file_basename(header.parent, '.fits') + '_subdark.fits'
; data = mrdfits('VL subdark/' + file, 0, /silent)
; data = rebin(data, header.naxis1, header.naxis2)
; data = data * header.nbin * header.ndit
; history = ['Dark correction: ', ' ' + file]
; ====================================
if header.pol_id eq 0 then begin
btype = 'VL fixed-polarization intensity'
descriptor = 'metis-vl-image'
polarimetric = 1
endif else if header.pol_id ge 1 and header.pol_id le 4 then begin
btype = 'VL fixed-polarization intensity'
descriptor = 'metis-vl-image'
polarimetric = 1
endif else begin
btype = 'VL total brightness'
descriptor = 'metis-vl-tb'
polarimetric = 0
endelse
data = metis_flat_field(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
data = metis_vignetting(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
data = metis_rad_cal(data, header, cal_pack, polarimetric = polarimetric, error = error, quality_matrix = quality_matrix, history = history)
bunit = cal_pack.vl_channel.cal_units
end
3: begin
; calibration of temporal noise images
btype = 'VL temporal standard deviation'
descriptor = file_name_fields[2]
bunit = 'DN'
history = !null
end
5: begin
; calibration of cr/sep log matrices
btype = 'VL cosmic-ray matrix'
descriptor = file_name_fields[2]
bunit = 'DN'
history = !null
end
else: begin
journal, 'Error 01: wrong input data product (expected data types 0, 3, or 5).'
journal
exit, status = 1
end
endcase
; definitions for the primary header ; version of the fits file
version = string(input.l2_version, format = '(I02)')
; fits creation date
date = date_conv(systime(/julian, /utc), 'FITS')
; name of the fits file
file_name = 'solo_L2_' + descriptor + '_' + file_name_fields[3] + '_V' + version + '.fits'
out_file_name = 'output/' + file_name
; adjust the primary header
fxaddpar, primary_header, 'FILENAME', file_name
fxaddpar, primary_header, 'PARENT', file_basename(input.file_name)
fxaddpar, primary_header, 'LEVEL', 'L2'
fxaddpar, primary_header, 'CREATOR', 'metis_l2_prep_vl_generic.pro'
fxaddpar, primary_header, 'VERS_SW', input.sw_version
fxaddpar, primary_header, 'VERS_CAL', cal_pack.version
fxaddpar, primary_header, 'VERSION', version
fxaddpar, primary_header, 'BTYPE', btype
fxaddpar, primary_header, 'BUNIT', bunit
fxaddpar, primary_header, 'DATAMIN', min(data, /nan)
fxaddpar, primary_header, 'DATAMAX', max(data, /nan)
fxaddpar, primary_header, 'WAVEBAND', cal_pack.vl_channel.name
fxaddpar, primary_header, 'XPOSURE', header.xposure
fxaddpar, primary_header, 'NSUMEXP', header.nsumexp
sxdelpar, primary_header, 'BLANK'
; fix planning info keywords
if header.soopname.startswith('unknown') then soopname = 'none' else soopname = header.soopname
if header.obs_mode.startswith('unknown') then obs_mode = 'none' else obs_mode = header.obs_mode
if soopname eq 'none' then sooptype = 'none' else sooptype = header.sooptype
if obs_mode eq 'none' then obs_type = 'none' else obs_type = header.obs_type
if soopname eq 'none' and obs_mode eq 'none' then obs_id = 'none' else obs_id = header.obs_id
if obs_mode eq 'none' then obs_mode = 'METIS_GENERIC_OBS'
fxaddpar, primary_header, 'SOOPNAME', soopname
fxaddpar, primary_header, 'SOOPTYPE', sooptype
fxaddpar, primary_header, 'OBS_MODE', obs_mode
fxaddpar, primary_header, 'OBS_TYPE', obs_type
fxaddpar, primary_header, 'OBS_ID', obs_id
; read the calibration curve to convert pmp raw voltages (dacpol) into effective polarization angles
dacpol_cal = cal_pack.vl_channel.dacpol_cal
if header.pol_id ne 5 then begin
if fix(header.hdr_vers) le 4 then begin
case header.pol_id of
0: dacpol = header.dac1Pol1
1: dacpol = header.dac1Pol1
2: dacpol = header.dac1Pol2
3: dacpol = header.dac1Pol3
4: dacpol = header.dac1Pol4
endcase
endif
if fix(header.hdr_vers) ge 5 then begin
dacpol = header.dac1Pol1
endif
k = where(dacpol_cal.dacpol eq dacpol)
angle = dacpol_cal.angle[k]
angle = float(angle[0])
fxaddpar, primary_header, 'POLANGLE', angle, '[deg] polarization angle', after = 'POL_ID'
endif
; append wcs keywords
wcs = metis_wcs(header, cal_pack, ref_detector = ref_detector)
foreach element, wcs do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
; append solar ephemeris keywords
ephemeris = solo_get_ephemeris(header, cal_pack)
foreach element, ephemeris do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
history = [history, 'Update WCS and solar ephemeris:', ' SKD version = ' + kernel_version]
; update the comment and history keywords
fxaddpar, primary_header, 'COMMENT', 'Uncertainty matrix in the FITS extension is preliminary.'
fxaddpar, primary_header, 'COMMENT', 'Rotate CROTA degrees counter-clockwise to have Solar North up.'
for k = 0, n_elements(history) - 1 do $
fxaddpar, primary_header, 'HISTORY', history[k]
fxaddpar, primary_header, 'HISTORY', 'L2 FITS file created on ' + date
; add checksum and datasum to the fits header
if not ref_detector then data = metis_rectify(data, 'VL')
fits_add_checksum, primary_header, float(data)
mwrfits, float(data), out_file_name, primary_header, /no_comment, /create, /silent
journal, 'Fits file created:'
journal, ' file name = ' + file_basename(out_file_name)
; add the extension with the quality matrix
base_header = primary_header
sxdelpar, base_header, 'SIMPLE'
sxdelpar, base_header, 'EXTEND'
sxdelpar, base_header, 'DATASUM'
sxdelpar, base_header, 'CHECKSUM'
sxdelpar, base_header, 'COMMENT'
sxdelpar, base_header, 'HISTORY'
extension_header = base_header
fxaddpar, extension_header, 'XTENSION', 'IMAGE', 'image extension', before = 'BITPIX'
fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'extension name', before = 'LONGSTRN'
fxaddpar, extension_header, 'BTYPE', 'Pixel quality'
fxaddpar, extension_header, 'BUNIT', 'None'
fxaddpar, extension_header, 'DATAMIN', min(quality_matrix, /nan)
fxaddpar, extension_header, 'DATAMAX', max(quality_matrix, /nan)
fxaddpar, extension_header, 'COMMENT', 'Quality matrix values:'
fxaddpar, extension_header, 'COMMENT', ' NaN = saturated or null L0 pixel counts'
fxaddpar, extension_header, 'COMMENT', ' 0 = unreliable pixel value'
fxaddpar, extension_header, 'COMMENT', ' 1 = good pixel'
if not ref_detector then quality_matrix = metis_rectify(quality_matrix, 'VL')
fits_add_checksum, extension_header, float(quality_matrix)
mwrfits, float(quality_matrix), out_file_name, extension_header, /no_comment, /silent
journal, 'Quality-matrix extension correctly added.'
; add the extension with the error matrix
if not ref_detector then error = metis_rectify(error, 'VL')
error_matrix = data * sqrt(error)
extension_header = base_header
fxaddpar, extension_header, 'XTENSION', 'IMAGE', 'image extension', before = 'BITPIX'
fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN' fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN'
fxaddpar, extension_header, 'BTYPE', 'Absolute error'
fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan)
fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan)
fits_add_checksum, extension_header, float(error_matrix)
mwrfits, float(error_matrix), out_file_name, extension_header, /no_comment, /silent
journal, 'Error-matrix extension correctly added.'
; fix some header keywords
fix_fits_header, out_file_name
; write the auxiliary information file
output = { $
file_name: out_file_name, $
l1_file_name: input.file_name, $
log_file_name: 'output/metis_l2_prep_log.txt' $
}
json_write, output, 'output/contents.json'
; unload the spice kernels
load_spice_kernels, kernel_list = kernel_list, /unload
journal, 'SPICE kernel files unloaded.'
; close the log
journal, 'Exiting without errors.'
journal
exit, status = 0
end