source: sasview/src/sas/sascalc/dataloader/readers/red2d_reader.py @ 83db1cc

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 83db1cc was 2f85af7, checked in by lewis, 7 years ago

Refactor red_2d_reader to use FileReader? class
  • Property mode set to 100644
File size: 11.8 KB
Line 
1"""
2    TXT/IGOR 2D Q Map file reader
3"""
4#####################################################################
5#This software was developed by the University of Tennessee as part of the
6#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
7#project funded by the US National Science Foundation.
8#See the license text in license.txt
9#copyright 2008, University of Tennessee
10######################################################################
11import os
12import numpy as np
13import math
14from sas.sascalc.dataloader.data_info import plottable_2D, DataInfo, Detector
15from sas.sascalc.dataloader.file_reader_base_class import FileReader
16from sas.sascalc.dataloader.loader_exceptions import FileContentsException
17
18# Look for unit converter
19has_converter = True
20try:
21    from sas.sascalc.data_util.nxsunit import Converter
22except:
23    has_converter = False
24
25
26def check_point(x_point):
27    """
28    check point validity
29    """
30    # set zero for non_floats
31    try:
32        return float(x_point)
33    except:
34        return 0
35
36
37class Reader(FileReader):
38    """ Simple data reader for Igor data files """
39    ## File type
40    type_name = "IGOR/DAT 2D Q_map"
41    ## Wildcards
42    type = ["IGOR/DAT 2D file in Q_map (*.dat)|*.DAT"]
43    ## Extension
44    ext = ['.DAT', '.dat']
45
46    def write(self, filename, data):
47        """
48        Write to .dat
49
50        :param filename: file name to write
51        :param data: data2D
52        """
53        import time
54        # Write the file
55        try:
56            fd = open(filename, 'w')
57            t = time.localtime()
58            time_str = time.strftime("%H:%M on %b %d %y", t)
59
60            header_str = "Data columns are Qx - Qy - I(Qx,Qy)\n\nASCII data"
61            header_str += " created at %s \n\n" % time_str
62            # simple 2D header
63            fd.write(header_str)
64            # write qx qy I values
65            for i in range(len(data.data)):
66                fd.write("%g  %g  %g\n" % (data.qx_data[i],
67                                            data.qy_data[i],
68                                           data.data[i]))
69        finally:
70            fd.close()
71
72    def get_file_contents(self):
73        # Read file
74        buf = self.f_open.read()
75        self.f_open.close()
76        # Instantiate data object
77        self.current_dataset = plottable_2D()
78        self.current_datainfo = DataInfo()
79        self.current_datainfo.filename = os.path.basename(self.f_open.name)
80        self.current_datainfo.detector.append(Detector())
81
82        # Get content
83        data_started = False
84
85        ## Defaults
86        lines = buf.split('\n')
87        x = []
88        y = []
89
90        wavelength = None
91        distance = None
92        transmission = None
93
94        pixel_x = None
95        pixel_y = None
96
97        is_info = False
98        is_center = False
99
100        # Remove the last lines before the for loop if the lines are empty
101        # to calculate the exact number of data points
102        count = 0
103        while (len(lines[len(lines) - (count + 1)].lstrip().rstrip()) < 1):
104            del lines[len(lines) - (count + 1)]
105            count = count + 1
106
107        #Read Header and find the dimensions of 2D data
108        line_num = 0
109        # Old version NIST files: 0
110        ver = 0
111        for line in lines:
112            line_num += 1
113            ## Reading the header applies only to IGOR/NIST 2D q_map data files
114            # Find setup info line
115            if is_info:
116                is_info = False
117                line_toks = line.split()
118                # Wavelength in Angstrom
119                try:
120                    wavelength = float(line_toks[1])
121                    # Units
122                    if has_converter == True and \
123                    self.current_datainfo.source.wavelength_unit != 'A':
124                        conv = Converter('A')
125                        wavelength = conv(wavelength,
126                                          units=self.current_datainfo.source.wavelength_unit)
127                except:
128                    #Not required
129                    pass
130                # Distance in mm
131                try:
132                    distance = float(line_toks[3])
133                    # Units
134                    if has_converter == True and self.current_datainfo.detector[0].distance_unit != 'm':
135                        conv = Converter('m')
136                        distance = conv(distance,
137                            units=self.current_datainfo.detector[0].distance_unit)
138                except:
139                    #Not required
140                    pass
141
142                # Distance in meters
143                try:
144                    transmission = float(line_toks[4])
145                except:
146                    #Not required
147                    pass
148
149            if line.count("LAMBDA") > 0:
150                is_info = True
151
152            # Find center info line
153            if is_center:
154                is_center = False
155                line_toks = line.split()
156                # Center in bin number
157                center_x = float(line_toks[0])
158                center_y = float(line_toks[1])
159
160            if line.count("BCENT") > 0:
161                is_center = True
162            # Check version
163            if line.count("Data columns") > 0:
164                if line.count("err(I)") > 0:
165                    ver = 1
166            # Find data start
167            if line.count("ASCII data") > 0:
168                data_started = True
169                continue
170
171            ## Read and get data.
172            if data_started == True:
173                line_toks = line.split()
174                if len(line_toks) == 0:
175                    #empty line
176                    continue
177                # the number of columns must be stayed same
178                col_num = len(line_toks)
179                break
180        # Make numpy array to remove header lines using index
181        lines_array = np.array(lines)
182
183        # index for lines_array
184        lines_index = np.arange(len(lines))
185
186        # get the data lines
187        data_lines = lines_array[lines_index >= (line_num - 1)]
188        # Now we get the total number of rows (i.e., # of data points)
189        row_num = len(data_lines)
190        # make it as list again to control the separators
191        data_list = " ".join(data_lines.tolist())
192        # split all data to one big list w/" "separator
193        data_list = data_list.split()
194
195        # Check if the size is consistent with data, otherwise
196        #try the tab(\t) separator
197        # (this may be removed once get the confidence
198        #the former working all cases).
199        if len(data_list) != (len(data_lines)) * col_num:
200            data_list = "\t".join(data_lines.tolist())
201            data_list = data_list.split()
202
203        # Change it(string) into float
204        #data_list = map(float,data_list)
205        data_list1 = map(check_point, data_list)
206
207        # numpy array form
208        data_array = np.array(data_list1)
209        # Redimesion based on the row_num and col_num,
210        #otherwise raise an error.
211        try:
212            data_point = data_array.reshape(row_num, col_num).transpose()
213        except:
214            msg = "red2d_reader can't read this file: Incorrect number of data points provided."
215            raise FileContentsException(msg)
216        ## Get the all data: Let's HARDcoding; Todo find better way
217        # Defaults
218        dqx_data = np.zeros(0)
219        dqy_data = np.zeros(0)
220        err_data = np.ones(row_num)
221        qz_data = np.zeros(row_num)
222        mask = np.ones(row_num, dtype=bool)
223        # Get from the array
224        qx_data = data_point[0]
225        qy_data = data_point[1]
226        data = data_point[2]
227        if ver == 1:
228            if col_num > (2 + ver):
229                err_data = data_point[(2 + ver)]
230        if col_num > (3 + ver):
231            qz_data = data_point[(3 + ver)]
232        if col_num > (4 + ver):
233            dqx_data = data_point[(4 + ver)]
234        if col_num > (5 + ver):
235            dqy_data = data_point[(5 + ver)]
236        #if col_num > (6 + ver): mask[data_point[(6 + ver)] < 1] = False
237        q_data = np.sqrt(qx_data*qx_data+qy_data*qy_data+qz_data*qz_data)
238
239        # Extra protection(it is needed for some data files):
240        # If all mask elements are False, put all True
241        if not mask.any():
242            mask[mask == False] = True
243
244        # Store limits of the image in q space
245        xmin = np.min(qx_data)
246        xmax = np.max(qx_data)
247        ymin = np.min(qy_data)
248        ymax = np.max(qy_data)
249
250        ## calculate the range of the qx and qy_data
251        x_size = math.fabs(xmax - xmin)
252        y_size = math.fabs(ymax - ymin)
253
254        # calculate the number of pixels in the each axes
255        npix_y = math.floor(math.sqrt(len(data)))
256        npix_x = math.floor(len(data) / npix_y)
257
258        # calculate the size of bins
259        xstep = x_size / (npix_x - 1)
260        ystep = y_size / (npix_y - 1)
261
262        # store x and y axis bin centers in q space
263        x_bins = np.arange(xmin, xmax + xstep, xstep)
264        y_bins = np.arange(ymin, ymax + ystep, ystep)
265
266        # get the limits of q values
267        xmin = xmin - xstep / 2
268        xmax = xmax + xstep / 2
269        ymin = ymin - ystep / 2
270        ymax = ymax + ystep / 2
271
272        #Store data in outputs
273        #TODO: Check the lengths
274        self.current_dataset.data = data
275        if (err_data == 1).all():
276            self.current_dataset.err_data = np.sqrt(np.abs(data))
277            self.current_dataset.err_data[self.current_dataset.err_data == 0.0] = 1.0
278        else:
279            self.current_dataset.err_data = err_data
280
281        self.current_dataset.qx_data = qx_data
282        self.current_dataset.qy_data = qy_data
283        self.current_dataset.q_data = q_data
284        self.current_dataset.mask = mask
285
286        self.current_dataset.x_bins = x_bins
287        self.current_dataset.y_bins = y_bins
288
289        self.current_dataset.xmin = xmin
290        self.current_dataset.xmax = xmax
291        self.current_dataset.ymin = ymin
292        self.current_dataset.ymax = ymax
293
294        self.current_datainfo.source.wavelength = wavelength
295
296        # Store pixel size in mm
297        self.current_datainfo.detector[0].pixel_size.x = pixel_x
298        self.current_datainfo.detector[0].pixel_size.y = pixel_y
299
300        # Store the sample to detector distance
301        self.current_datainfo.detector[0].distance = distance
302
303        # optional data: if all of dq data == 0, do not pass to output
304        if len(dqx_data) == len(qx_data) and dqx_data.any() != 0:
305            # if no dqx_data, do not pass dqy_data.
306            #(1 axis dq is not supported yet).
307            if len(dqy_data) == len(qy_data) and dqy_data.any() != 0:
308                # Currently we do not support dq parr, perp.
309                # tranfer the comp. to cartesian coord. for newer version.
310                if ver != 1:
311                    diag = np.sqrt(qx_data * qx_data + qy_data * qy_data)
312                    cos_th = qx_data / diag
313                    sin_th = qy_data / diag
314                    self.current_dataset.dqx_data = np.sqrt((dqx_data * cos_th) * \
315                                                 (dqx_data * cos_th) \
316                                                 + (dqy_data * sin_th) * \
317                                                  (dqy_data * sin_th))
318                    self.current_dataset.dqy_data = np.sqrt((dqx_data * sin_th) * \
319                                                 (dqx_data * sin_th) \
320                                                 + (dqy_data * cos_th) * \
321                                                  (dqy_data * cos_th))
322                else:
323                    self.current_dataset.dqx_data = dqx_data
324                    self.current_dataset.dqy_data = dqy_data
325
326        # Units of axes
327        self.current_dataset.xaxis("\\rm{Q_{x}}", 'A^{-1}')
328        self.current_dataset.yaxis("\\rm{Q_{y}}", 'A^{-1}')
329        self.current_dataset.zaxis("\\rm{Intensity}", "cm^{-1}")
330
331        # Store loading process information
332        self.current_datainfo.meta_data['loader'] = self.type_name
333
334        self.send_to_output()
Note: See TracBrowser for help on using the repository browser.