source: sasview/src/sas/sascalc/dataloader/readers/cansas_reader_HDF5.py @ 90d389d

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalc
Last change on this file since 90d389d was 132db16, checked in by Piotr Rozyczko <rozyczko@…>, 8 years ago

Refactored Cansas HDF5 reader to allow multiple instances of SASdata of any dimensionality in a single SASentry and multiple instances of SASentry in a single file. Added unit tests to test all types of data (1 entry | 1 data, 1 entry | many data, many entry | 1 data each, many entr |, many data).
  • Property mode set to 100644
File size: 20.5 KB
Line 
1"""
2    CanSAS 2D data reader for reading HDF5 formatted CanSAS files.
3"""
4
5import h5py
6import numpy as np
7import re
8import os
9import sys
10
11from sas.sascalc.dataloader.data_info import plottable_1D, plottable_2D, Data1D, Data2D, DataInfo, Process, Aperture
12from sas.sascalc.dataloader.data_info import Collimation, TransmissionSpectrum, Detector
13from sas.sascalc.dataloader.data_info import combine_data_info_with_plottable
14
15
16
17class Reader():
18    """
19    A class for reading in CanSAS v2.0 data files. The existing iteration opens Mantid generated HDF5 formatted files
20    with file extension .h5/.H5. Any number of data sets may be present within the file and any dimensionality of data
21    may be used. Currently 1D and 2D SAS data sets are supported, but future implementations will include 1D and 2D
22    SESANS data.
23
24    Any number of SASdata sets may be present in a SASentry and the data within can be either 1D I(Q) or 2D I(Qx, Qy).
25
26    :Dependencies:
27        The CanSAS HDF5 reader requires h5py => v2.5.0 or later.
28    """
29
30    ## CanSAS version
31    cansas_version = 2.0
32    ## Logged warnings or messages
33    logging = None
34    ## List of errors for the current data set
35    errors = None
36    ## Raw file contents to be processed
37    raw_data = None
38    ## Data info currently being read in
39    current_datainfo = None
40    ## SASdata set currently being read in
41    current_dataset = None
42    ## List of plottable1D objects that should be linked to the current_datainfo
43    data1d = None
44    ## List of plottable2D objects that should be linked to the current_datainfo
45    data2d = None
46    ## Data type name
47    type_name = "CanSAS 2.0"
48    ## Wildcards
49    type = ["CanSAS 2.0 HDF5 Files (*.h5)|*.h5"]
50    ## List of allowed extensions
51    ext = ['.h5', '.H5']
52    ## Flag to bypass extension check
53    allow_all = False
54    ## List of files to return
55    output = None
56
57    def read(self, filename):
58        """
59        This is the general read method that all SasView data_loaders must have.
60
61        :param filename: A path for an HDF5 formatted CanSAS 2D data file.
62        :return: List of Data1D/2D objects and/or a list of errors.
63        """
64
65        ## Reinitialize the class when loading a new data file to reset all class variables
66        self.reset_class_variables()
67        ## Check that the file exists
68        if os.path.isfile(filename):
69            basename = os.path.basename(filename)
70            _, extension = os.path.splitext(basename)
71            # If the file type is not allowed, return empty list
72            if extension in self.ext or self.allow_all:
73                ## Load the data file
74                self.raw_data = h5py.File(filename, 'r')
75                ## Read in all child elements of top level SASroot
76                self.read_children(self.raw_data, [])
77                ## Add the last data set to the list of outputs
78                self.add_data_set()
79        ## Return data set(s)
80        return self.output
81
82    def reset_class_variables(self):
83        """
84        Create the reader object and define initial states for class variables
85        """
86        self.current_datainfo = None
87        self.current_dataset = None
88        self.data1d = []
89        self.data2d = []
90        self.raw_data = None
91        self.errors = set()
92        self.logging = []
93        self.output = []
94        self.parent_class = u''
95        self.detector = Detector()
96        self.collimation = Collimation()
97        self.aperture = Aperture()
98        self.process = Process()
99        self.trans_spectrum = TransmissionSpectrum()
100
101    def read_children(self, data, parent_list):
102        """
103        A recursive method for stepping through the hierarchical data file.
104
105        :param data: h5py Group object of any kind
106        :param parent: h5py Group parent name
107        """
108
109        ## Loop through each element of the parent and process accordingly
110        for key in data.keys():
111            ## Get all information for the current key
112            value = data.get(key)
113            if value.attrs.get(u'canSAS_class') is not None:
114                class_name = value.attrs.get(u'canSAS_class')
115            else:
116                class_name = value.attrs.get(u'NX_class')
117            if class_name is not None:
118                class_prog = re.compile(class_name)
119            else:
120                class_prog = re.compile(value.name)
121
122            if isinstance(value, h5py.Group):
123                self.parent_class = class_name
124                parent_list.append(key)
125                ## If this is a new sasentry, store the current data sets and create a fresh Data1D/2D object
126                if class_prog.match(u'SASentry'):
127                    self.add_data_set(key)
128                elif class_prog.match(u'SASdata'):
129                    self._initialize_new_data_set(parent_list)
130                ## Recursion step to access data within the group
131                self.read_children(value, parent_list)
132                self.add_intermediate()
133                parent_list.remove(key)
134
135            elif isinstance(value, h5py.Dataset):
136                ## If this is a dataset, store the data appropriately
137                data_set = data[key][:]
138
139                for data_point in data_set:
140                    ## Top Level Meta Data
141                    unit = self._get_unit(value)
142                    if key == u'definition':
143                        self.current_datainfo.meta_data['reader'] = data_point
144                    elif key == u'run':
145                        self.current_datainfo.run.append(data_point)
146                    elif key == u'title':
147                        self.current_datainfo.title = data_point
148                    elif key == u'SASnote':
149                        self.current_datainfo.notes.append(data_point)
150
151                    ## I and Q Data
152                    elif key == u'I':
153                        if type(self.current_dataset) is plottable_2D:
154                            self.current_dataset.data = np.append(self.current_dataset.data, data_point)
155                            self.current_dataset.zaxis("Intensity", unit)
156                        else:
157                            self.current_dataset.y = np.append(self.current_dataset.y, data_point)
158                            self.current_dataset.yaxis("Intensity", unit)
159                    elif key == u'Idev':
160                        if type(self.current_dataset) is plottable_2D:
161                            self.current_dataset.err_data = np.append(self.current_dataset.err_data, data_point)
162                        else:
163                            self.current_dataset.dy = np.append(self.current_dataset.dy, data_point)
164                    elif key == u'Q':
165                        self.current_dataset.xaxis("Q", unit)
166                        if type(self.current_dataset) is plottable_2D:
167                            self.current_dataset.q = np.append(self.current_dataset.q, data_point)
168                        else:
169                            self.current_dataset.x = np.append(self.current_dataset.x, data_point)
170                    elif key == u'Qy':
171                        self.current_dataset.yaxis("Q_y", unit)
172                        self.current_dataset.qy_data = np.append(self.current_dataset.qy_data, data_point)
173                    elif key == u'Qydev':
174                        self.current_dataset.dqy_data = np.append(self.current_dataset.dqy_data, data_point)
175                    elif key == u'Qx':
176                        self.current_dataset.xaxis("Q_x", unit)
177                        self.current_dataset.qx_data = np.append(self.current_dataset.qx_data, data_point)
178                    elif key == u'Qxdev':
179                        self.current_dataset.dqx_data = np.append(self.current_dataset.dqx_data, data_point)
180                    elif key == u'Mask':
181                        self.current_dataset.mask = np.append(self.current_dataset.mask, data_point)
182
183                    ## Sample Information
184                    elif key == u'Title' and self.parent_class == u'SASsample':
185                        self.current_datainfo.sample.name = data_point
186                    elif key == u'thickness' and self.parent_class == u'SASsample':
187                        self.current_datainfo.sample.thickness = data_point
188                    elif key == u'temperature' and self.parent_class == u'SASsample':
189                        self.current_datainfo.sample.temperature = data_point
190
191                    ## Instrumental Information
192                    elif key == u'name' and self.parent_class == u'SASinstrument':
193                        self.current_datainfo.instrument = data_point
194                    elif key == u'name' and self.parent_class == u'SASdetector':
195                        self.detector.name = data_point
196                    elif key == u'SDD' and self.parent_class == u'SASdetector':
197                        self.detector.distance = float(data_point)
198                        self.detector.distance_unit = unit
199                    elif key == u'SSD' and self.parent_class == u'SAScollimation':
200                        self.collimation.length = data_point
201                        self.collimation.length_unit = unit
202                    elif key == u'name' and self.parent_class == u'SAScollimation':
203                        self.collimation.name = data_point
204
205                    ## Process Information
206                    elif key == u'name' and self.parent_class == u'SASprocess':
207                        self.process.name = data_point
208                    elif key == u'Title' and self.parent_class == u'SASprocess':
209                        self.process.name = data_point
210                    elif key == u'description' and self.parent_class == u'SASprocess':
211                        self.process.description = data_point
212                    elif key == u'date' and self.parent_class == u'SASprocess':
213                        self.process.date = data_point
214                    elif self.parent_class == u'SASprocess':
215                        self.process.notes.append(data_point)
216
217                    ## Transmission Spectrum
218                    elif key == u'T' and self.parent_class == u'SAStransmission_spectrum':
219                        self.trans_spectrum.transmission.append(data_point)
220                    elif key == u'Tdev' and self.parent_class == u'SAStransmission_spectrum':
221                        self.trans_spectrum.transmission_deviation.append(data_point)
222                    elif key == u'lambda' and self.parent_class == u'SAStransmission_spectrum':
223                        self.trans_spectrum.wavelength.append(data_point)
224
225                    ## Other Information
226                    elif key == u'wavelength' and self.parent_class == u'SASdata':
227                        self.current_datainfo.source.wavelength = data_point
228                        self.current_datainfo.source.wavelength.unit = unit
229                    elif key == u'radiation' and self.parent_class == u'SASsource':
230                        self.current_datainfo.source.radiation = data_point
231                    elif key == u'transmission' and self.parent_class == u'SASdata':
232                        self.current_datainfo.sample.transmission = data_point
233
234                    ## Everything else goes in meta_data
235                    else:
236                        new_key = self._create_unique_key(self.current_datainfo.meta_data, key)
237                        self.current_datainfo.meta_data[new_key] = data_point
238
239            else:
240                ## I don't know if this reachable code
241                self.errors.add("ShouldNeverHappenException")
242
243    def add_intermediate(self):
244        """
245        This method stores any intermediate objects within the final data set after fully reading the set.
246
247        :param parent: The NXclass name for the h5py Group object that just finished being processed
248        """
249
250        if self.parent_class == u'SASprocess':
251            self.current_datainfo.process.append(self.process)
252            self.process = Process()
253        elif self.parent_class == u'SASdetector':
254            self.current_datainfo.detector.append(self.detector)
255            self.detector = Detector()
256        elif self.parent_class == u'SAStransmission_spectrum':
257            self.current_datainfo.trans_spectrum.append(self.trans_spectrum)
258            self.trans_spectrum = TransmissionSpectrum()
259        elif self.parent_class == u'SAScollimation':
260            self.current_datainfo.collimation.append(self.collimation)
261            self.collimation = Collimation()
262        elif self.parent_class == u'SASaperture':
263            self.collimation.aperture.append(self.aperture)
264            self.aperture = Aperture()
265        elif self.parent_class == u'SASdata':
266            if type(self.current_dataset) is plottable_2D:
267                self.data2d.append(self.current_dataset)
268            elif type(self.current_dataset) is plottable_1D:
269                self.data1d.append(self.current_dataset)
270
271    def final_data_cleanup(self):
272        """
273        Does some final cleanup and formatting on self.current_datainfo and all data1D and data2D objects and then
274        combines the data and info into Data1D and Data2D objects
275        """
276
277        ## Type cast data arrays to float64
278        if len(self.current_datainfo.trans_spectrum) > 0:
279            spectrum_list = []
280            for spectrum in self.current_datainfo.trans_spectrum:
281                spectrum.transmission = np.delete(spectrum.transmission, [0])
282                spectrum.transmission = spectrum.transmission.astype(np.float64)
283                spectrum.transmission_deviation = np.delete(spectrum.transmission_deviation, [0])
284                spectrum.transmission_deviation = spectrum.transmission_deviation.astype(np.float64)
285                spectrum.wavelength = np.delete(spectrum.wavelength, [0])
286                spectrum.wavelength = spectrum.wavelength.astype(np.float64)
287                if len(spectrum.transmission) > 0:
288                    spectrum_list.append(spectrum)
289            self.current_datainfo.trans_spectrum = spectrum_list
290
291        ## Append errors to dataset and reset class errors
292        self.current_datainfo.errors = self.errors
293        self.errors.clear()
294
295        ## Combine all plottables with datainfo and append each to output
296        ## Type cast data arrays to float64 and find min/max as appropriate
297        for dataset in self.data2d:
298            dataset.data = np.delete(dataset.data, [0])
299            dataset.data = dataset.data.astype(np.float64)
300            dataset.err_data = np.delete(dataset.err_data, [0])
301            dataset.err_data = dataset.err_data.astype(np.float64)
302            dataset.mask = np.delete(dataset.mask, [0])
303            if dataset.qx_data is not None:
304                dataset.qx_data = np.delete(dataset.qx_data, [0])
305                dataset.xmin = np.min(dataset.qx_data)
306                dataset.xmax = np.max(dataset.qx_data)
307                dataset.qx_data = dataset.qx_data.astype(np.float64)
308            if dataset.dqx_data is not None:
309                dataset.dqx_data = np.delete(dataset.dqx_data, [0])
310                dataset.dqx_data = dataset.dqx_data.astype(np.float64)
311            if dataset.qy_data is not None:
312                dataset.qy_data = np.delete(dataset.qy_data, [0])
313                dataset.ymin = np.min(dataset.qy_data)
314                dataset.ymax = np.max(dataset.qy_data)
315                dataset.qy_data = dataset.qy_data.astype(np.float64)
316            if dataset.dqy_data is not None:
317                dataset.dqy_data = np.delete(dataset.dqy_data, [0])
318                dataset.dqy_data = dataset.dqy_data.astype(np.float64)
319            if dataset.q_data is not None:
320                dataset.q_data = np.delete(dataset.q_data, [0])
321                dataset.q_data = dataset.q_data.astype(np.float64)
322            zeros = np.ones(dataset.data.size, dtype=bool)
323            try:
324                for i in range (0, dataset.mask.size - 1):
325                    zeros[i] = dataset.mask[i]
326            except:
327                self.errors.add(sys.exc_value)
328            dataset.mask = zeros
329            ## Calculate the actual Q matrix
330            try:
331                if dataset.q_data.size <= 1:
332                    dataset.q_data = np.sqrt(dataset.qx_data * dataset.qx_data + dataset.qy_data * dataset.qy_data)
333            except:
334                dataset.q_data = None
335            final_dataset = combine_data_info_with_plottable(dataset, self.current_datainfo)
336            self.output.append(final_dataset)
337
338        for dataset in self.data1d:
339            if dataset.x is not None:
340                dataset.x = np.delete(dataset.x, [0])
341                dataset.x = dataset.x.astype(np.float64)
342                dataset.xmin = np.min(dataset.x)
343                dataset.xmax = np.max(dataset.x)
344            if dataset.y is not None:
345                dataset.y = np.delete(dataset.y, [0])
346                dataset.y = dataset.y.astype(np.float64)
347                dataset.ymin = np.min(dataset.y)
348                dataset.ymax = np.max(dataset.y)
349            if dataset.dx is not None:
350                dataset.dx = np.delete(dataset.dx, [0])
351                dataset.dx = dataset.dx.astype(np.float64)
352            if dataset.dxl is not None:
353                dataset.dxl = np.delete(dataset.dxl, [0])
354                dataset.dxl = dataset.dxl.astype(np.float64)
355            if dataset.dxw is not None:
356                dataset.dxw = np.delete(dataset.dxw, [0])
357                dataset.dxw = dataset.dxw.astype(np.float64)
358            if dataset.dy is not None:
359                dataset.dy = np.delete(dataset.dy, [0])
360                dataset.dy = dataset.dy.astype(np.float64)
361            final_dataset = combine_data_info_with_plottable(dataset, self.current_datainfo)
362            self.output.append(final_dataset)
363
364    def add_data_set(self, key=""):
365        """
366        Adds the current_dataset to the list of outputs after preforming final processing on the data and then calls a
367        private method to generate a new data set.
368
369        :param key: NeXus group name for current tree level
370        """
371
372        if self.current_datainfo and self.current_dataset:
373            self.final_data_cleanup()
374        self.data1d = []
375        self.data2d = []
376        self.current_datainfo = DataInfo()
377
378    def _initialize_new_data_set(self, parent_list = None):
379        """
380        A private class method to generate a new 1D or 2D data object based on the type of data within the set.
381        Outside methods should call add_data_set() to be sure any existing data is stored properly.
382
383        :param parent_list: List of names of parent elements
384        """
385
386        if parent_list is None:
387            parent_list = []
388        if self._find_intermediate(parent_list, "Qx"):
389            self.current_dataset = plottable_2D()
390        else:
391            x = np.array(0)
392            y = np.array(0)
393            self.current_dataset = plottable_1D(x, y)
394        self.current_datainfo.filename = self.raw_data.filename
395
396    def _find_intermediate(self, parent_list, basename=""):
397        """
398        A private class used to find an entry by either using a direct key or knowing the approximate basename.
399
400        :param parent_list: List of parents to the current level in the HDF5 file
401        :param basename: Approximate name of an entry to search for
402        :return:
403        """
404
405        entry = False
406        key_prog = re.compile(basename)
407        top = self.raw_data
408        for parent in parent_list:
409            top = top.get(parent)
410        for key in top.keys():
411            if (key_prog.match(key)):
412                entry = True
413                break
414        return entry
415
416    def _create_unique_key(self, dictionary, name, numb=0):
417        """
418        Create a unique key value for any dictionary to prevent overwriting
419        Recurses until a unique key value is found.
420
421        :param dictionary: A dictionary with any number of entries
422        :param name: The index of the item to be added to dictionary
423        :param numb: The number to be appended to the name, starts at 0
424        :return: The new name for the dictionary entry
425        """
426        if dictionary.get(name) is not None:
427            numb += 1
428            name = name.split("_")[0]
429            name += "_{0}".format(numb)
430            name = self._create_unique_key(dictionary, name, numb)
431        return name
432
433    def _get_unit(self, value):
434        """
435        Find the unit for a particular value within the h5py dictionary
436
437        :param value: attribute dictionary for a particular value set
438        :return: unit for the value passed to the method
439        """
440        unit = value.attrs.get(u'units')
441        if unit == None:
442            unit = value.attrs.get(u'unit')
443        ## Convert the unit formats
444        if unit == "1/A":
445            unit = "A^{-1}"
446        elif unit == "1/cm":
447            unit = "cm^{-1}"
448        return unit
Note: See TracBrowser for help on using the repository browser.