| 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 | |
|---|
| 10 | """ |
|---|
| 11 | TXT/IGOR 2D Q Map file reader |
|---|
| 12 | """ |
|---|
| 13 | |
|---|
| 14 | |
|---|
| 15 | import os |
|---|
| 16 | #import sys |
|---|
| 17 | import numpy |
|---|
| 18 | import math |
|---|
| 19 | #import logging |
|---|
| 20 | from DataLoader.data_info import Data2D, Detector |
|---|
| 21 | |
|---|
| 22 | # Look for unit converter |
|---|
| 23 | has_converter = True |
|---|
| 24 | try: |
|---|
| 25 | from data_util.nxsunit import Converter |
|---|
| 26 | except: |
|---|
| 27 | has_converter = False |
|---|
| 28 | |
|---|
| 29 | def 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 | |
|---|
| 39 | class 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'] |
|---|
| 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() |
|---|
| 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 |
|---|
| 79 | |
|---|
| 80 | # Read file |
|---|
| 81 | f = open(filename,'r') |
|---|
| 82 | buf = f.read() |
|---|
| 83 | f.close() |
|---|
| 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 | |
|---|
| 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 |
|---|
| 140 | for line in lines: |
|---|
| 141 | line_num += 1 |
|---|
| 142 | ## Reading the header applies only to IGOR/NIST 2D q_map data files |
|---|
| 143 | # Find setup info line |
|---|
| 144 | if isInfo: |
|---|
| 145 | isInfo = False |
|---|
| 146 | line_toks = line.split() |
|---|
| 147 | # Wavelength in Angstrom |
|---|
| 148 | try: |
|---|
| 149 | wavelength = float(line_toks[1]) |
|---|
| 150 | # Units |
|---|
| 151 | if has_converter == True and \ |
|---|
| 152 | output.source.wavelength_unit != 'A': |
|---|
| 153 | conv = Converter('A') |
|---|
| 154 | wavelength = conv(wavelength, |
|---|
| 155 | units=output.source.wavelength_unit) |
|---|
| 156 | except: |
|---|
| 157 | #Not required |
|---|
| 158 | pass |
|---|
| 159 | # Distance in mm |
|---|
| 160 | try: |
|---|
| 161 | distance = float(line_toks[3]) |
|---|
| 162 | # Units |
|---|
| 163 | if has_converter == True and detector.distance_unit != 'm': |
|---|
| 164 | conv = Converter('m') |
|---|
| 165 | distance = conv(distance, units=detector.distance_unit) |
|---|
| 166 | except: |
|---|
| 167 | #Not required |
|---|
| 168 | pass |
|---|
| 169 | |
|---|
| 170 | # Distance in meters |
|---|
| 171 | try: |
|---|
| 172 | transmission = float(line_toks[4]) |
|---|
| 173 | except: |
|---|
| 174 | #Not required |
|---|
| 175 | pass |
|---|
| 176 | |
|---|
| 177 | if line.count("LAMBDA") > 0: |
|---|
| 178 | isInfo = True |
|---|
| 179 | |
|---|
| 180 | # Find center info line |
|---|
| 181 | if isCenter: |
|---|
| 182 | isCenter = False |
|---|
| 183 | line_toks = line.split() |
|---|
| 184 | # Center in bin number |
|---|
| 185 | center_x = float(line_toks[0]) |
|---|
| 186 | center_y = float(line_toks[1]) |
|---|
| 187 | |
|---|
| 188 | if line.count("BCENT") > 0: |
|---|
| 189 | isCenter = True |
|---|
| 190 | |
|---|
| 191 | # Find data start |
|---|
| 192 | if line.count("Data columns") or line.count("ASCII data") > 0: |
|---|
| 193 | dataStarted = True |
|---|
| 194 | continue |
|---|
| 195 | |
|---|
| 196 | ## Read and get data. |
|---|
| 197 | if dataStarted == True: |
|---|
| 198 | line_toks = line.split() |
|---|
| 199 | if len(line_toks) == 0: |
|---|
| 200 | #empty line |
|---|
| 201 | continue |
|---|
| 202 | # the number of columns must be stayed same |
|---|
| 203 | col_num = len(line_toks) |
|---|
| 204 | break |
|---|
| 205 | # Make numpy array to remove header lines using index |
|---|
| 206 | lines_array = numpy.array(lines) |
|---|
| 207 | |
|---|
| 208 | # index for lines_array |
|---|
| 209 | lines_index = numpy.arange(len(lines)) |
|---|
| 210 | |
|---|
| 211 | # get the data lines |
|---|
| 212 | data_lines = lines_array[lines_index >= (line_num - 1)] |
|---|
| 213 | # Now we get the total number of rows (i.e., # of data points) |
|---|
| 214 | row_num = len(data_lines) |
|---|
| 215 | # make it as list again to control the separators |
|---|
| 216 | data_list = " ".join(data_lines.tolist()) |
|---|
| 217 | # split all data to one big list w/" "separator |
|---|
| 218 | data_list = data_list.split() |
|---|
| 219 | |
|---|
| 220 | # Check if the size is consistent with data, otherwise |
|---|
| 221 | #try the tab(\t) separator |
|---|
| 222 | # (this may be removed once get the confidence |
|---|
| 223 | #the former working all cases). |
|---|
| 224 | if len(data_list) != (len(data_lines)) * col_num: |
|---|
| 225 | data_list = "\t".join(data_lines.tolist()) |
|---|
| 226 | data_list = data_list.split() |
|---|
| 227 | |
|---|
| 228 | # Change it(string) into float |
|---|
| 229 | #data_list = map(float,data_list) |
|---|
| 230 | data_list1 = map(check_point,data_list) |
|---|
| 231 | |
|---|
| 232 | # numpy array form |
|---|
| 233 | data_array = numpy.array(data_list1) |
|---|
| 234 | # Redimesion based on the row_num and col_num, |
|---|
| 235 | #otherwise raise an error. |
|---|
| 236 | try: |
|---|
| 237 | data_point = data_array.reshape(row_num, col_num).transpose() |
|---|
| 238 | except: |
|---|
| 239 | msg = "red2d_reader: Can't read this file: Not a proper file format" |
|---|
| 240 | raise ValueError, msg |
|---|
| 241 | |
|---|
| 242 | ## Get the all data: Let's HARDcoding; Todo find better way |
|---|
| 243 | # Defaults |
|---|
| 244 | dqx_data = numpy.zeros(0) |
|---|
| 245 | dqy_data = numpy.zeros(0) |
|---|
| 246 | qz_data = numpy.zeros(row_num) |
|---|
| 247 | mask = numpy.ones(row_num,dtype=bool) |
|---|
| 248 | # Get from the array |
|---|
| 249 | qx_data = data_point[0] |
|---|
| 250 | qy_data = data_point[1] |
|---|
| 251 | data = data_point[2] |
|---|
| 252 | if col_num > 3: qz_data = data_point[3] |
|---|
| 253 | if col_num > 4: dqx_data = data_point[4] |
|---|
| 254 | if col_num > 5: dqy_data = data_point[5] |
|---|
| 255 | if col_num > 6: mask[data_point[6] < 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(): mask[mask==False] = True |
|---|
| 261 | |
|---|
| 262 | # Store limits of the image in q space |
|---|
| 263 | xmin = numpy.min(qx_data) |
|---|
| 264 | xmax = numpy.max(qx_data) |
|---|
| 265 | ymin = numpy.min(qy_data) |
|---|
| 266 | ymax = numpy.max(qy_data) |
|---|
| 267 | |
|---|
| 268 | # units |
|---|
| 269 | if has_converter == True and output.Q_unit != '1/A': |
|---|
| 270 | xmin = data_conv_q(xmin, units=output.Q_unit) |
|---|
| 271 | xmax = data_conv_q(xmax, units=output.Q_unit) |
|---|
| 272 | ymin = data_conv_q(ymin, units=output.Q_unit) |
|---|
| 273 | ymax = data_conv_q(ymax, units=output.Q_unit) |
|---|
| 274 | |
|---|
| 275 | ## calculate the range of the qx and qy_data |
|---|
| 276 | x_size = math.fabs(xmax - xmin) |
|---|
| 277 | y_size = math.fabs(ymax - ymin) |
|---|
| 278 | |
|---|
| 279 | # calculate the number of pixels in the each axes |
|---|
| 280 | npix_y = math.floor(math.sqrt(len(data))) |
|---|
| 281 | npix_x = math.floor(len(data)/npix_y) |
|---|
| 282 | |
|---|
| 283 | # calculate the size of bins |
|---|
| 284 | xstep = x_size/(npix_x-1) |
|---|
| 285 | ystep = y_size/(npix_y-1) |
|---|
| 286 | |
|---|
| 287 | # store x and y axis bin centers in q space |
|---|
| 288 | x_bins = numpy.arange(xmin,xmax+xstep,xstep) |
|---|
| 289 | y_bins = numpy.arange(ymin,ymax+ystep,ystep) |
|---|
| 290 | |
|---|
| 291 | # get the limits of q values |
|---|
| 292 | xmin = xmin - xstep/2 |
|---|
| 293 | xmax = xmax + xstep/2 |
|---|
| 294 | ymin = ymin - ystep/2 |
|---|
| 295 | ymax = ymax + ystep/2 |
|---|
| 296 | |
|---|
| 297 | #Store data in outputs |
|---|
| 298 | #TODO: Check the lengths |
|---|
| 299 | output.data = data |
|---|
| 300 | output.err_data = numpy.sqrt(numpy.abs(data)) |
|---|
| 301 | output.qx_data = qx_data |
|---|
| 302 | output.qy_data = qy_data |
|---|
| 303 | output.q_data = q_data |
|---|
| 304 | output.mask = mask |
|---|
| 305 | |
|---|
| 306 | output.x_bins = x_bins |
|---|
| 307 | output.y_bins = y_bins |
|---|
| 308 | |
|---|
| 309 | output.xmin = xmin |
|---|
| 310 | output.xmax = xmax |
|---|
| 311 | output.ymin = ymin |
|---|
| 312 | output.ymax = ymax |
|---|
| 313 | |
|---|
| 314 | output.source.wavelength = wavelength |
|---|
| 315 | |
|---|
| 316 | # Store pixel size in mm |
|---|
| 317 | detector.pixel_size.x = pixel_x |
|---|
| 318 | detector.pixel_size.y = pixel_y |
|---|
| 319 | |
|---|
| 320 | # Store the sample to detector distance |
|---|
| 321 | detector.distance = distance |
|---|
| 322 | |
|---|
| 323 | # optional data: if all of dq data == 0, do not pass to output |
|---|
| 324 | if len(dqx_data) == len(qx_data) and dqx_data.any()!=0: |
|---|
| 325 | # if no dqx_data, do not pass dqy_data. |
|---|
| 326 | #(1 axis dq is not supported yet). |
|---|
| 327 | if len(dqy_data) == len(qy_data) and dqy_data.any()!=0: |
|---|
| 328 | output.dqx_data = dqx_data |
|---|
| 329 | output.dqy_data = dqy_data |
|---|
| 330 | |
|---|
| 331 | # Units of axes |
|---|
| 332 | if data_conv_q is not None: |
|---|
| 333 | output.xaxis("\\rm{Q_{x}}", output.Q_unit) |
|---|
| 334 | output.yaxis("\\rm{Q_{y}}", output.Q_unit) |
|---|
| 335 | else: |
|---|
| 336 | output.xaxis("\\rm{Q_{x}}", 'A^{-1}') |
|---|
| 337 | output.yaxis("\\rm{Q_{y}}", 'A^{-1}') |
|---|
| 338 | if data_conv_i is not None: |
|---|
| 339 | output.zaxis("\\rm{Intensity}", output.I_unit) |
|---|
| 340 | else: |
|---|
| 341 | output.zaxis("\\rm{Intensity}","cm^{-1}") |
|---|
| 342 | |
|---|
| 343 | # Store loading process information |
|---|
| 344 | output.meta_data['loader'] = self.type_name |
|---|
| 345 | |
|---|
| 346 | return output |
|---|
| 347 | |
|---|
| 348 | if __name__ == "__main__": |
|---|
| 349 | reader = Reader() |
|---|
| 350 | print reader.read("../test/exp18_14_igor_2dqxqy.dat") |
|---|
| 351 | |
|---|
| 352 | |
|---|
