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

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.1.1release-4.1.2release-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since be88076 was be88076, checked in by lewis, 8 years ago

Correctly read and write run_name

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