source: sasview/src/sas/sascalc/dataloader/readers/cansas_reader_HDF5.py @ 67b0a99

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 67b0a99 was bbd0f37, checked in by krzywon, 8 years ago

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