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

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalcmagnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 5582b078 was 2b538cd, checked in by krzywon, 7 years ago

Modify documentation and add instance check for numpy arrays in HDF5 reader.

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