| 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 | ###################################################################### |
|---|
| 11 | import os |
|---|
| 12 | import numpy |
|---|
| 13 | import math |
|---|
| 14 | from sas.dataloader.data_info import Data2D, Detector |
|---|
| 15 | |
|---|
| 16 | # Look for unit converter |
|---|
| 17 | has_converter = True |
|---|
| 18 | try: |
|---|
| 19 | from sas.data_util.nxsunit import Converter |
|---|
| 20 | except: |
|---|
| 21 | has_converter = False |
|---|
| 22 | |
|---|
| 23 | |
|---|
| 24 | def 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 | |
|---|
| 35 | class 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 |
|---|
