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

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.1.2release-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 3563e06 was 7c24685, checked in by krzywon, 7 years ago

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