source: sasview/sansdataloader/src/sans/dataloader/readers/red2d_reader.py @ aac129aa

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 aac129aa was ad8034f, checked in by Jae Cho <jhjcho@…>, 13 years ago

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