source: sasview/src/sas/sascalc/dataloader/readers/red2d_reader.py @ 3d164b9

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.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 3d164b9 was b699768, checked in by Piotr Rozyczko <piotr.rozyczko@…>, 9 years ago

Initial commit of the refactored SasCalc? module.

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