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

ESS_GUIESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalc
Last change on this file since b8080e1 was b8080e1, checked in by Piotr Rozyczko <rozyczko@…>, 6 years ago

cherry picking sascalc changes from master SASVIEW-996
minor unit test fixes

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