source: sasview/src/sas/sascalc/dataloader/readers/cansas_reader_HDF5.py @ 6964d44

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 6964d44 was 9687d58, checked in by Piotr Rozyczko <rozyczko@…>, 7 years ago

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