-
Notifications
You must be signed in to change notification settings - Fork 1
/
load_EISCAT_hdf5.m
185 lines (167 loc) · 6.71 KB
/
load_EISCAT_hdf5.m
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
function [data, metadata] = load_EISCAT_hdf5(varargin)
%% load EISCAT-type hdf5 file.
% Author: L. Cai.
% Create date: 2021-07-20
% the hdf5 file can be downloaded from eiscat.se
% Input:
% file_path: the file full path include the directory and file name
% Output:
% data: structure, selected variables from the hdf5 file. The full
% list of variables can be checked using the function "list_EISCAT_hdf5(file_path)."
% .n_e: electron density
% .n_e_err: errors of the electron density
% .T_i: ion temperature
% .T_i_err: errors of the ion temperature
% ... similar for variables in var_name_dict below
% metadata: structure, selected metadata of the hdf5
% .name_site
% .site
% .pulse_code
% .antenna
% ... the fields can be extended by adding more attributes on
% lines below.
file_path = './example/EISCAT_2021-03-10_beata_ant@uhfa.hdf5';
if length(varargin) == 1
file_path = varargin{1};
end
% Set the variable list, which will be stored in "data"
var_name_dict = struct( ...
'r_XMITloc1', 'XMITloc1', ...
'r_XMITloc2', 'XMITloc2', ...
'r_XMITloc3', 'XMITloc3', ...
'r_RECloc1', 'RECloc1', ...
'r_RECloc2', 'RECloc2', ...
'r_RECloc3', 'RECloc3', ...
'magic_constant', 'Magic_const', ...
'r_SCangle', 'SCangle', ...
'r_m0_2', 'm02', ...
'r_m0_1', 'm01', ...
'az', 'az', ...
'el', 'el', ...
'P_Tx', 'Pt', ...
'height', 'h', ...
'range', 'range', ...
'n_e', 'Ne', ...
'T_i', 'Ti', ...
'T_r', 'Tr', ...
'nu_i', 'Collf', ...
'v_i_los', 'Vi', ...
'comp_mix', 'pm', ...
'comp_O_p', 'po+', ...
'n_e_var', 'var_Ne', ...
'T_i_var', 'var_Ti', ...
'T_r_var', 'var_Tr', ...
'nu_i_var', 'var_Collf', ...
'v_i_los_var', 'var_Vi', ...
'comp_mix_var', 'var_pm', ...
'comp_O_p_var', 'var_po+', ...
'status', 'status', ...
'residual', 'res1' ...
);
% Set the site information for convert the sigle letter to 3-letter
% site code
site_info = struct( ...
'T', 'UHF', ...
'V', 'VHF', ...
'K', 'KIR', ...
'S', 'SOD', ...
'L', 'ESR' ...
);
% load h5 data
h5data.par0d = h5read(file_path, '/data/par0d');
h5data.par1d = h5read(file_path, '/data/par1d');
try
h5data.par2d = h5read(file_path, '/data/par2d');
catch
;
end
% get variable information
try
var_info_list = list_EISCAT_hdf5_variables(file_path);
catch
var_info_list = list_EISCAT_hdf5_variables(file_path,[],{'par0d', 'par1d'});
end
% get the gate numbers (nrec) for each beam. Note: nrec can be a single
% integar or a 1-D array
ind_nrec = strcmp(var_info_list.name, 'nrec');
nrec_group = var_info_list.group{ind_nrec};
nrec_index = var_info_list.index(ind_nrec);
nrec = h5data.(nrec_group)(:, nrec_index);
% load each variables
field_names = fieldnames(var_name_dict);
[ncol, ans] = size(h5data.par1d);
for i=1 : length(field_names)
var_name = field_names{i};
var_name_hdf5 = var_name_dict.(var_name);
% get variable information
ix = strcmp(var_info_list.name, var_name_hdf5);
var_group = var_info_list.group{ix};
var_index = var_info_list.index(ix);
var_value = h5data.(var_group)(:, var_index);
if strcmp(var_group, 'par0d') % for a single number
data.(var_name) = var_value;
elseif strcmp(var_group, 'par1d') % for a 1-D array
data.(var_name) = reshape(var_value, 1, ncol);
elseif strcmp(var_group, 'par2d') % for a 2-D array
if strcmp(nrec_group, 'par0d') % for nrec is a single number
nrow = length(var_value) / ncol;
if nrow ~= nrec
disp("Note: the number of range gates doesn't match nrec!")
end
data.(var_name) = reshape(var_value, nrow, ncol);
elseif strcmp(nrec_group, 'par1d') % for nrec is a 1-D array
num_gates = max(nrec);
var_array = nan(num_gates, ncol);
rec_ind_1 = 1;
for j = 1 : ncol % reshape the 1-D array to a 2-D array
rec_ind_2 = rec_ind_1 + nrec(j) - 1;
var_array(1:nrec(j), j) = var_value(rec_ind_1:rec_ind_2);
rec_ind_1 = rec_ind_2 + 1;
end
data.(var_name) = var_array;
end
end
end
field_names = fieldnames(data);
for i = 1: length(field_names)
var_name = field_names{i};
if contains(var_name, '_var')
var_name_new = replace(var_name, '_var', '_err');
data.(var_name_new) = sqrt(data.(var_name));
end
end
% convert the unix times to datetime
utimes = h5read(file_path, '/data/utime');
utime_1 = utimes(:, 1);
datetime_1 = datetime(utime_1,'ConvertFrom','epochtime','TicksPerSecond',1,'Format','dd-MMM-yyyy HH:mm:ss');
data.datevec_1 = datevec(datetime_1)';
data.datenum_1 = datenum(datetime_1)';
utime_2 = utimes(:, 2);
datetime_2 = datetime(utime_2,'ConvertFrom','epochtime','TicksPerSecond',1,'Format','dd-MMM-yyyy HH:mm:ss');
data.datevec_2 = datevec(datetime_2)';
data.datenum_2 = datenum(datetime_2)';
% Calculate T_e and error
data.T_e = data.T_i .* data.T_r;
data.T_e_err = data.T_e .* sqrt((data.T_i_err./data.T_i).^2 + ...
(data.T_r_err./data.T_r).^2);
data.height = data.height / 1000.;
data.range = data.range / 1000.;
% resize az and el
if length(data.az) == 1
data.az = repmat(data.az, 1, ncol);
end
if length(data.el) == 1
data.el = repmat(data.el, 1, ncol);
end
% load the metadata
names = convertCharsToStrings(strtrim(h5read(file_path, '/metadata/names')));
metadata.name_site = names{2, 2};
metadata.name_expr = names{2, 1};
metadata.site = site_info.(metadata.name_site);
metadata.pulse_code = names{2, 1};
metadata.antenna = names{2, 3};
metadata.r_XMITloc = [data.r_XMITloc1 data.r_XMITloc2 data.r_XMITloc3/1000.];
metadata.r_RECloc = [data.r_RECloc1 data.r_RECloc2 data.r_RECloc3/1000.];
metadata.title = h5read(file_path, '/metadata/schemes/DataCite/Title');
metadata.rawdata_path = h5read(file_path, '/metadata/software/gfd/data_path');
end