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

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.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since beba407 was 7f75a3f, checked in by krzywon, 8 years ago

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