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

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 5a405bd was dcb91cf, checked in by lewis, 7 years ago

Make suggested changes

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