source: sasview/src/sas/sascalc/dataloader/readers/cansas_reader_HDF5.py @ 082239e

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 082239e was 082239e, checked in by krzywon, 7 years ago

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