- Timestamp:
- Aug 1, 2017 5:19:46 AM (7 years ago)
- Branches:
- master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.2.2, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
- Children:
- b2c28a5
- Parents:
- a78433dd
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/sas/sascalc/dataloader/readers/danse_reader.py
r235f514 r713a047 5 5 #This software was developed by the University of Tennessee as part of the 6 6 #Distributed Data Analysis of Neutron Scattering Experiments (DANSE) 7 #project funded by the US National Science Foundation. 7 #project funded by the US National Science Foundation. 8 8 #If you use DANSE applications to do scientific research that leads to 9 9 #publication, we ask that you acknowledge the use of the software with the … … 14 14 import math 15 15 import os 16 import sys17 16 import numpy as np 18 17 import logging 19 from sas.sascalc.dataloader.data_info import Data2D, Detector18 from sas.sascalc.dataloader.data_info import plottable_2D, DataInfo, Detector 20 19 from sas.sascalc.dataloader.manipulations import reader2D_converter 20 from sas.sascalc.dataloader.file_reader_base_class import FileReader 21 from sas.sascalc.dataloader.loader_exceptions import FileContentsException, DataReaderException 21 22 22 23 logger = logging.getLogger(__name__) … … 30 31 31 32 32 class Reader :33 class Reader(FileReader): 33 34 """ 34 35 Example data manipulation … … 40 41 ## Extension 41 42 ext = ['.sans', '.SANS'] 42 43 def read(self, filename=None): 44 """ 45 Open and read the data in a file 46 @param file: path of the file 47 """ 48 49 read_it = False 50 for item in self.ext: 51 if filename.lower().find(item) >= 0: 52 read_it = True 53 54 if read_it: 43 44 def get_file_contents(self): 45 self.current_datainfo = DataInfo() 46 self.current_dataset = plottable_2D() 47 self.output = [] 48 49 loaded_correctly = True 50 error_message = "" 51 52 # defaults 53 # wavelength in Angstrom 54 wavelength = 10.0 55 # Distance in meter 56 distance = 11.0 57 # Pixel number of center in x 58 center_x = 65 59 # Pixel number of center in y 60 center_y = 65 61 # Pixel size [mm] 62 pixel = 5.0 63 # Size in x, in pixels 64 size_x = 128 65 # Size in y, in pixels 66 size_y = 128 67 # Format version 68 fversion = 1.0 69 70 self.current_datainfo.filename = os.path.basename(self.f_open.name) 71 detector = Detector() 72 self.current_datainfo.detector.append(detector) 73 74 self.current_dataset.data = np.zeros([size_x, size_y]) 75 self.current_dataset.err_data = np.zeros([size_x, size_y]) 76 77 read_on = True 78 data_start_line = 1 79 while read_on: 80 line = self.f_open.readline() 81 data_start_line += 1 82 if line.find("DATA:") >= 0: 83 read_on = False 84 break 85 toks = line.split(':') 55 86 try: 56 datafile = open(filename, 'r')57 except:58 raise RuntimeError,"danse_reader cannot open %s" % (filename)59 60 # defaults61 # wavelength in Angstrom62 wavelength = 10.063 # Distance in meter64 distance = 11.065 # Pixel number of center in x66 center_x = 6567 # Pixel number of center in y68 center_y = 6569 # Pixel size [mm]70 pixel = 5.071 # Size in x, in pixels72 size_x = 12873 # Size in y, in pixels74 size_y = 12875 # Format version76 fversion = 1.077 78 output = Data2D()79 output.filename = os.path.basename(filename)80 detector = Detector()81 output.detector.append(detector)82 83 output.data = np.zeros([size_x,size_y])84 output.err_data = np.zeros([size_x, size_y])85 86 data_conv_q = None87 data_conv_i = None88 89 if has_converter == True and output.Q_unit != '1/A':90 data_conv_q = Converter('1/A')91 # Test it92 data_conv_q(1.0, output.Q_unit)93 94 if has_converter == True and output.I_unit != '1/cm':95 data_conv_i = Converter('1/cm')96 # Test it97 data_conv_i(1.0, output.I_unit)98 99 read_on = True100 while read_on:101 line = datafile.readline()102 if line.find("DATA:") >= 0:103 read_on = False104 break105 toks = line.split(':')106 87 if toks[0] == "FORMATVERSION": 107 88 fversion = float(toks[1]) 108 if toks[0] == "WAVELENGTH":89 elif toks[0] == "WAVELENGTH": 109 90 wavelength = float(toks[1]) 110 91 elif toks[0] == "DISTANCE": … … 120 101 elif toks[0] == "SIZE_Y": 121 102 size_y = int(toks[1]) 122 123 # Read the data 124 data = [] 125 error = [] 126 if fversion == 1.0: 127 data_str = datafile.readline() 128 data = data_str.split(' ') 129 else: 130 read_on = True 131 while read_on: 132 data_str = datafile.readline() 133 if len(data_str) == 0: 134 read_on = False 135 else: 136 toks = data_str.split() 137 try: 138 val = float(toks[0]) 139 err = float(toks[1]) 140 if data_conv_i is not None: 141 val = data_conv_i(val, units=output._yunit) 142 err = data_conv_i(err, units=output._yunit) 143 data.append(val) 144 error.append(err) 145 except: 146 logger.info("Skipping line:%s,%s" %(data_str, 147 sys.exc_value)) 148 149 # Initialize 150 x_vals = [] 151 y_vals = [] 152 ymin = None 153 ymax = None 154 xmin = None 155 xmax = None 156 157 # Qx and Qy vectors 158 theta = pixel / distance / 100.0 159 stepq = 4.0 * math.pi / wavelength * math.sin(theta / 2.0) 160 for i_x in range(size_x): 161 theta = (i_x - center_x + 1) * pixel / distance / 100.0 162 qx = 4.0 * math.pi / wavelength * math.sin(theta / 2.0) 163 164 if has_converter == True and output.Q_unit != '1/A': 165 qx = data_conv_q(qx, units=output.Q_unit) 166 167 x_vals.append(qx) 168 if xmin is None or qx < xmin: 169 xmin = qx 170 if xmax is None or qx > xmax: 171 xmax = qx 172 173 ymin = None 174 ymax = None 175 for i_y in range(size_y): 176 theta = (i_y - center_y + 1) * pixel / distance / 100.0 177 qy = 4.0 * math.pi / wavelength * math.sin(theta/2.0) 178 179 if has_converter == True and output.Q_unit != '1/A': 180 qy = data_conv_q(qy, units=output.Q_unit) 181 182 y_vals.append(qy) 183 if ymin is None or qy < ymin: 184 ymin = qy 185 if ymax is None or qy > ymax: 186 ymax = qy 187 188 # Store the data in the 2D array 189 i_x = 0 190 i_y = -1 191 192 for i_pt in range(len(data)): 193 try: 194 value = float(data[i_pt]) 195 except: 196 # For version 1.0, the data were still 197 # stored as strings at this point. 198 msg = "Skipping entry (v1.0):%s,%s" % (str(data[i_pt]), 199 sys.exc_value) 200 logger.info(msg) 201 202 # Get bin number 203 if math.fmod(i_pt, size_x) == 0: 204 i_x = 0 205 i_y += 1 206 else: 207 i_x += 1 208 209 output.data[i_y][i_x] = value 210 if fversion>1.0: 211 output.err_data[i_y][i_x] = error[i_pt] 212 213 # Store all data 214 # Store wavelength 215 if has_converter == True and output.source.wavelength_unit != 'A': 216 conv = Converter('A') 217 wavelength = conv(wavelength, 218 units=output.source.wavelength_unit) 219 output.source.wavelength = wavelength 220 221 # Store distance 222 if has_converter == True and detector.distance_unit != 'm': 223 conv = Converter('m') 224 distance = conv(distance, units=detector.distance_unit) 225 detector.distance = distance 226 227 # Store pixel size 228 if has_converter == True and detector.pixel_size_unit != 'mm': 229 conv = Converter('mm') 230 pixel = conv(pixel, units=detector.pixel_size_unit) 231 detector.pixel_size.x = pixel 232 detector.pixel_size.y = pixel 233 234 # Store beam center in distance units 235 detector.beam_center.x = center_x * pixel 236 detector.beam_center.y = center_y * pixel 237 238 # Store limits of the image (2D array) 239 xmin = xmin - stepq / 2.0 240 xmax = xmax + stepq / 2.0 241 ymin = ymin - stepq /2.0 242 ymax = ymax + stepq / 2.0 243 244 if has_converter == True and output.Q_unit != '1/A': 245 xmin = data_conv_q(xmin, units=output.Q_unit) 246 xmax = data_conv_q(xmax, units=output.Q_unit) 247 ymin = data_conv_q(ymin, units=output.Q_unit) 248 ymax = data_conv_q(ymax, units=output.Q_unit) 249 output.xmin = xmin 250 output.xmax = xmax 251 output.ymin = ymin 252 output.ymax = ymax 253 254 # Store x and y axis bin centers 255 output.x_bins = x_vals 256 output.y_bins = y_vals 257 258 # Units 259 if data_conv_q is not None: 260 output.xaxis("\\rm{Q_{x}}", output.Q_unit) 261 output.yaxis("\\rm{Q_{y}}", output.Q_unit) 262 else: 263 output.xaxis("\\rm{Q_{x}}", 'A^{-1}') 264 output.yaxis("\\rm{Q_{y}}", 'A^{-1}') 265 266 if data_conv_i is not None: 267 output.zaxis("\\rm{Intensity}", output.I_unit) 268 else: 269 output.zaxis("\\rm{Intensity}", "cm^{-1}") 270 271 if not fversion >= 1.0: 272 msg = "Danse_reader can't read this file %s" % filename 273 raise ValueError, msg 274 else: 275 logger.info("Danse_reader Reading %s \n" % filename) 276 277 # Store loading process information 278 output.meta_data['loader'] = self.type_name 279 output = reader2D_converter(output) 280 return output 281 282 return None 103 except ValueError as e: 104 error_message += "Unable to parse {}. Default value used.\n".format(toks[0]) 105 loaded_correctly = False 106 107 # Read the data 108 data = [] 109 error = [] 110 if not fversion >= 1.0: 111 msg = "danse_reader can't read this file {}".format(self.f_open.name) 112 raise FileContentsException(msg) 113 114 for line_num, data_str in enumerate(self.f_open.readlines()): 115 toks = data_str.split() 116 try: 117 val = float(toks[0]) 118 err = float(toks[1]) 119 data.append(val) 120 error.append(err) 121 except ValueError as exc: 122 msg = "Unable to parse line {}: {}".format(line_num + data_start_line, data_str.strip()) 123 raise FileContentsException(msg) 124 125 num_pts = size_x * size_y 126 if len(data) < num_pts: 127 msg = "Not enough data points provided. Expected {} but got {}".format( 128 size_x * size_y, len(data)) 129 raise FileContentsException(msg) 130 elif len(data) > num_pts: 131 error_message += ("Too many data points provided. Expected {0} but" 132 " got {1}. Only the first {0} will be used.\n").format(num_pts, len(data)) 133 loaded_correctly = False 134 data = data[:num_pts] 135 error = error[:num_pts] 136 137 # Qx and Qy vectors 138 theta = pixel / distance / 100.0 139 i_x = np.arange(size_x) 140 theta = (i_x - center_x + 1) * pixel / distance / 100.0 141 x_vals = 4.0 * np.pi / wavelength * np.sin(theta / 2.0) 142 xmin = x_vals.min() 143 xmax = x_vals.max() 144 145 i_y = np.arange(size_y) 146 theta = (i_y - center_y + 1) * pixel / distance / 100.0 147 y_vals = 4.0 * np.pi / wavelength * np.sin(theta / 2.0) 148 ymin = y_vals.min() 149 ymax = y_vals.max() 150 151 self.current_dataset.data = np.array(data, dtype=np.float64).reshape((size_y, size_x)) 152 if fversion > 1.0: 153 self.current_dataset.err_data = np.array(error, dtype=np.float64).reshape((size_y, size_x)) 154 155 # Store all data 156 # Store wavelength 157 if has_converter == True and self.current_datainfo.source.wavelength_unit != 'A': 158 conv = Converter('A') 159 wavelength = conv(wavelength, 160 units=self.current_datainfo.source.wavelength_unit) 161 self.current_datainfo.source.wavelength = wavelength 162 163 # Store distance 164 if has_converter == True and detector.distance_unit != 'm': 165 conv = Converter('m') 166 distance = conv(distance, units=detector.distance_unit) 167 detector.distance = distance 168 169 # Store pixel size 170 if has_converter == True and detector.pixel_size_unit != 'mm': 171 conv = Converter('mm') 172 pixel = conv(pixel, units=detector.pixel_size_unit) 173 detector.pixel_size.x = pixel 174 detector.pixel_size.y = pixel 175 176 # Store beam center in distance units 177 detector.beam_center.x = center_x * pixel 178 detector.beam_center.y = center_y * pixel 179 180 181 self.current_dataset.xaxis("\\rm{Q_{x}}", 'A^{-1}') 182 self.current_dataset.yaxis("\\rm{Q_{y}}", 'A^{-1}') 183 self.current_dataset.zaxis("\\rm{Intensity}", "cm^{-1}") 184 185 self.current_dataset.x_bins = x_vals 186 self.current_dataset.y_bins = y_vals 187 188 # Reshape data 189 x_vals = np.tile(x_vals, (size_y, 1)).flatten() 190 y_vals = np.tile(y_vals, (size_x, 1)).T.flatten() 191 if self.current_dataset.err_data == np.all(np.array(None)) or np.any(self.current_dataset.err_data <= 0): 192 new_err_data = np.sqrt(np.abs(self.current_dataset.data)) 193 else: 194 new_err_data = self.current_dataset.err_data.flatten() 195 196 self.current_dataset.err_data = new_err_data 197 self.current_dataset.qx_data = x_vals 198 self.current_dataset.qy_data = y_vals 199 self.current_dataset.q_data = np.sqrt(x_vals**2 + y_vals**2) 200 self.current_dataset.mask = np.ones(len(x_vals), dtype=bool) 201 202 # Store loading process information 203 self.current_datainfo.meta_data['loader'] = self.type_name 204 205 self.send_to_output() 206 207 if not loaded_correctly: 208 raise DataReaderException(error_message)
Note: See TracChangeset
for help on using the changeset viewer.
