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

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 e5cbbce was c94280c, checked in by krzywon, 8 years ago

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