| 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 sans.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 | # Old version NIST files: 0 |
|---|
| 141 | ver = 0 |
|---|
| 142 | for line in lines: |
|---|
| 143 | line_num += 1 |
|---|
| 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 |
|---|
| 153 | if has_converter == True and \ |
|---|
| 154 | output.source.wavelength_unit != 'A': |
|---|
| 155 | conv = Converter('A') |
|---|
| 156 | wavelength = conv(wavelength, |
|---|
| 157 | units=output.source.wavelength_unit) |
|---|
| 158 | except: |
|---|
| 159 | #Not required |
|---|
| 160 | pass |
|---|
| 161 | # Distance in mm |
|---|
| 162 | try: |
|---|
| 163 | distance = float(line_toks[3]) |
|---|
| 164 | # Units |
|---|
| 165 | if has_converter == True and detector.distance_unit != 'm': |
|---|
| 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 | |
|---|
| 179 | if line.count("LAMBDA") > 0: |
|---|
| 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 | |
|---|
| 190 | if line.count("BCENT") > 0: |
|---|
| 191 | isCenter = True |
|---|
| 192 | # Check version |
|---|
| 193 | if line.count("Data columns") > 0: |
|---|
| 194 | if line.count("err(I)") > 0: |
|---|
| 195 | ver = 1 |
|---|
| 196 | # Find data start |
|---|
| 197 | if line.count("ASCII data") > 0: |
|---|
| 198 | dataStarted = True |
|---|
| 199 | continue |
|---|
| 200 | |
|---|
| 201 | ## Read and get data. |
|---|
| 202 | if dataStarted == True: |
|---|
| 203 | line_toks = line.split() |
|---|
| 204 | if len(line_toks) == 0: |
|---|
| 205 | #empty line |
|---|
| 206 | continue |
|---|
| 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 |
|---|
| 217 | data_lines = lines_array[lines_index >= (line_num - 1)] |
|---|
| 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 | |
|---|
| 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). |
|---|
| 229 | if len(data_list) != (len(data_lines)) * col_num: |
|---|
| 230 | data_list = "\t".join(data_lines.tolist()) |
|---|
| 231 | data_list = data_list.split() |
|---|
| 232 | |
|---|
| 233 | # Change it(string) into float |
|---|
| 234 | #data_list = map(float,data_list) |
|---|
| 235 | data_list1 = map(check_point,data_list) |
|---|
| 236 | |
|---|
| 237 | # numpy array form |
|---|
| 238 | data_array = numpy.array(data_list1) |
|---|
| 239 | # Redimesion based on the row_num and col_num, |
|---|
| 240 | #otherwise raise an error. |
|---|
| 241 | try: |
|---|
| 242 | data_point = data_array.reshape(row_num, col_num).transpose() |
|---|
| 243 | except: |
|---|
| 244 | msg = "red2d_reader: Can't read this file: Not a proper file format" |
|---|
| 245 | raise ValueError, msg |
|---|
| 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) |
|---|
| 250 | err_data = numpy.ones(row_num) |
|---|
| 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] |
|---|
| 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 |
|---|
| 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): |
|---|
| 266 | # If all mask elements are False, put all True |
|---|
| 267 | if not mask.any(): mask[mask==False] = True |
|---|
| 268 | |
|---|
| 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 |
|---|
| 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 | |
|---|
| 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 | |
|---|
| 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: |
|---|
| 337 | # if no dqx_data, do not pass dqy_data. |
|---|
| 338 | #(1 axis dq is not supported yet). |
|---|
| 339 | if len(dqy_data) == len(qy_data) and dqy_data.any()!=0: |
|---|
| 340 | # Currently we do not support dq parr, perp. |
|---|
| 341 | # tranfer the comp. to cartesian coord. for newer version. |
|---|
| 342 | if ver != 1: |
|---|
| 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 | |
|---|
| 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 | |
|---|
| 375 | if __name__ == "__main__": |
|---|
| 376 | reader = Reader() |
|---|
| 377 | print reader.read("../test/exp18_14_igor_2dqxqy.dat") |
|---|
| 378 | |
|---|
| 379 | |
|---|
