fieldpy.stream_gauging.calibration

1 """ 2 This file contains the functions to read calibration files. 3 4 Calibration file format is almost identical to TOA5: 5 "TIMESTAMP","CONCENTRATION","SENSOR_EC_READING" 6 "2010-07-13 08:50:00",0.00,0.234 7 """ 8 import numpy as np 9 import datetime 10 import pylab as plt 11 import matplotlib.mlab as mlab 12 13 14 from fieldpy.core.experiment_classes import * 15 from fieldpy.core import raw_file_readers 16

17 -class AllCalibrations(ExperimentData):

18 """ 19 This class can hold several calibrations and provides the 20 following functions to convert sensor readout to concentration: 21 22 - L{readout2concentration_closest} 23 - L{readout2concentration_weighted} 24 """

25 - def __init__(self, list_of_filenames, plot_yes=False):

26 """Initialises the class with 27 28 @type list_of_filenames: list of strings 29 @param list_of_filenames: The file names of the calibration files in .maw file format 30 31 @type plot_yes: boolean 32 @param plot_yes: Plot the linear regression if True (default = False) 33 34 >>> from fieldpy.stream_gauging.calibration import * 35 >>> dir_ = '../stream_gauging/test_data/' 36 >>> fielns = ['calibration_13JUL10.maw', 'calibration_18JUL10_1.maw', 'calibration_18JUL10_2.maw'] 37 >>> ac = AllCalibrations([dir_ + fl for fl in fielns]) 38 Linear regression of conductance vs concentration 39 a=345.72 b=-0.08, root mean squared error= 0.000, R^2 = 0.997 40 Linear regression of conductance vs concentration 41 a=363.19 b=-0.10, root mean squared error= 0.001, R^2 = 0.996 42 Linear regression of conductance vs concentration 43 a=254.78 b=-0.03, root mean squared error= 0.000, R^2 = 0.999 44 """ 45 super(AllCalibrations, self).__init__() 46 self.calibrations = [] 47 self._plot_yes = plot_yes 48 for filename in list_of_filenames: 49 self.calibrations.append(SingleCalibration(filename, self._plot_yes)) 50 # sort them according to when the calibartion was done 51 cmpfn = lambda x, y : cmp(x._datetime, y._datetime) 52 self.calibrations.sort(cmp=cmpfn) 53 self.check()

54

55 - def __repr__(self):

56 st = 'AllCalibrations instance containing:\n' 57 for ca in self.calibrations: 58 st = st + ca.__str__() + '\n' 59 return st[:-1]

60 - def check(self):

61 """ 62 Does integrity checks: 63 - check that all SingleCalibration.md.reference_resistor are the same 64 """ 65 super(AllCalibrations, self).check() 66 ref_r = self.calibrations[0].md.reference_resistor 67 for cali in self.calibrations[1:]: 68 if cali.md.reference_resistor!=ref_r: 69 raise ValueError('The reference resistors of the individual calibrations are not equal!')

70 - def plot_all(self):

71 """ 72 Plots all calibrations. 73 """ 74 axs = self.calibrations[0].plot_it() 75 axs[0].hold(True) 76 axs[1].hold(True) 77 for ca in self.calibrations[1:]: 78 ca.plot_it(axs) 79 axs[0].set_title('Linear Regression between conductance and concentration')

80

81 - def conductance2concentration_average(self, times, conductance, 82 background_cond, verbose_=False):

83 """ 84 Converts the conductance into a concentration using an average of 85 the calibartions. 86 87 @type times: numpy array matplotlib time stamps 88 @param times: times of sensor readout 89 90 @type conductance: numpy array 91 @param conductance: sensor readouts converted to conductance 92 93 @type background_cond: float 94 @param background_cond: the conductance at zero salt concentration 95 96 @type verbose_: boolean 97 @param verbose_: print which calibration was used 98 99 @rtype: numpy array or float 100 @return: concentration 101 """ 102 # average the self.calibrations.md.ar coefficient 103 ar_aver = 0 104 for caa in self.calibrations: 105 ar_aver += caa.md.ar 106 ar_aver = ar_aver/len(self.calibrations) 107 cali2use = SingleCalibration(None) 108 cali2use.md.ar = ar_aver 109 cali2use.md.reference_resistor = self.calibrations[0].md.reference_resistor 110 if verbose_: 111 print 'Used slope is %f' % ar_aver 112 return cali2use.conductance2concentration(conductance, background_cond)

113 114

115 - def readout2concentration_average(self, times, readouts, background_readout, verbose_=False):

116 """ 117 Converts the readout into a concentration using an average of 118 the calibartions. 119 120 @type times: numpy array matplotlib time stamps 121 @param times: times of sensor readout 122 123 @type readouts: numpy array 124 @param readouts: sensor readouts 125 126 @type background_readout: float 127 @param background_readout: the sensor readout at zero salt concentration (default=1) 128 129 @type verbose_: boolean 130 @param verbose_: print which calibration was used 131 132 @rtype: numpy array or float 133 @return: concentration 134 """ 135 # average the self.calibrations.md.ar coefficient 136 ar_aver = 0 137 for caa in self.calibrations: 138 ar_aver += caa.md.ar 139 ar_aver = ar_aver/len(self.calibrations) 140 cali2use = SingleCalibration(None) 141 cali2use.md.ar = ar_aver 142 cali2use.md.reference_resistor = self.calibrations[0].md.reference_resistor 143 if verbose_: 144 print 'Used slope is %f' % ar_aver 145 return cali2use.readout2concentration(readouts, background_readout)

146

147 - def readout2concentration_closest(self, times, readouts, background_readout, verbose_=False):

148 """ 149 Converts the readout into a concentration using the 150 calibration closest in time to the readouts (note only the 151 time of first readout is checked). 152 153 @type times: numpy array matplotlib time stamps 154 @param times: times of sensor readout 155 156 @type readouts: numpy array 157 @param readouts: sensor readouts 158 159 @type background_readout: float 160 @param background_readout: the sensor readout at zero salt concentration (default=1) 161 162 @type verbose_: boolean 163 @param verbose_: print which calibration was used 164 165 @rtype: numpy array or float 166 @return: concentration 167 """ 168 # figure out which calibration to use 169 dt = [np.abs(times[0]-ca.data['time'][0]) for ca in self.calibrations] 170 cali2use = self.calibrations[np.argmin(dt)] 171 if verbose_: 172 print 'Used calibration: ' + cali2use.__repr__() + '\nwith slope ' + cali2use.md.ar.__str__() 173 return cali2use.readout2concentration(readouts, background_readout)

174 - def conductance2concentration_closest(self, times, readouts, background_readout, verbose_=False):

175 """ 176 Converts the readout into a concentration using the 177 calibration closest in time to the readouts (note only the 178 time of first readout is checked). 179 180 @type times: numpy array matplotlib time stamps 181 @param times: times of sensor readout 182 183 @type readouts: numpy array 184 @param readouts: sensor readouts 185 186 @type background_readout: float 187 @param background_readout: the sensor readout at zero salt concentration (default=1) 188 189 @type verbose_: boolean 190 @param verbose_: print which calibration was used 191 192 @rtype: numpy array or float 193 @return: concentration 194 """ 195 # figure out which calibration to use 196 dt = [np.abs(times[0]-ca.data['time'][0]) for ca in self.calibrations] 197 cali2use = self.calibrations[np.argmin(dt)] 198 if verbose_: 199 print 'Used calibration: ' + cali2use.__repr__() + '\nwith slope ' + cali2use.md.ar.__str__() 200 return cali2use.conductance2concentration(readouts, background_readout)

201

202 - def readout2concentration_weighted(self, times, readouts, background_readout):

203 """ 204 Converts the readout into a concentration using the weighted 205 calibration of the two calibrations closest in time to the readout. 206 207 @type times: numpy array 208 @param times: times of sensor readout 209 210 @type readouts: numpy array 211 @param readouts: sensor readouts 212 213 @type background_readout: float 214 @param background_readout: the sensor readout at zero salt concentration (default=1) 215 216 @rtype: numpy array or float 217 @return: concentration 218 """ 219 pass

220 221 222

223 -class SingleCalibration(LoggedData):

224 """ 225 This class holds all the methods to read and process the salt 226 calibation of the electical conductivity (EC) sensor. 227 """

228 - def __init__(self, filename=None, plot_yes=True):

229 """Initialises the class. If filename is set to None, the 230 class needs to be initialised by hand. 231 232 @type filename: string 233 @param filename: The file name of the calibration file in .maw file format 234 235 @type plot_yes: boolean 236 @param plot_yes: Plot the linear regression if True (default) 237 238 >>> from fieldpy.stream_gauging.calibration import * 239 >>> sc = SingleCalibration('../stream_gauging/test_data/calibration_13JUL10.maw', plot_yes=False) 240 Linear regression of conductance vs concentration 241 a=345.72 b=-0.08, root mean squared error= 0.000, R^2 = 0.997 242 >>> sc 243 Conductivity sensor calibration at 2010-07-13 10:30 244 """ 245 super(SingleCalibration, self).__init__() 246 if filename is not None: 247 self.filename = filename 248 self.data, self.raw_data, self.md = raw_file_readers.read_maw_file(filename) 249 self._datetime = plt.num2date(self.data['time'][0]) 250 self._plot_yes = plot_yes 251 self.make_calibration() 252 self.check() 253 else: 254 return

255

256 - def __repr__(self):

257 try: 258 return 'Conductivity sensor calibration at %s' % datetime.datetime.strftime(self._datetime, '%Y-%m-%d %H:%M') 259 except: 260 return 'Conductivity sensor calibration at no time'

261 - def __str__(self):

262 try: 263 return 'Conductivity sensor calibration at %s' % datetime.datetime.strftime(self._datetime, '%Y-%m-%d %H:%M') 264 except: 265 return 'Conductivity sensor calibration at no time'

266

267 - def readout2conductance(self, readout):

268 """ 269 Calcualtes the conductance across the sensor electrodes for a 270 given sensor readout. 271 272 @type readout: float or numpy array of floats 273 @param readout: the sensor readout 274 275 @rtype: numpy array or float 276 @return: conductance 277 """ 278 return (1-readout)/(self.md.reference_resistor*readout)

279

280 - def readout2concentration(self, readout, background_readout=None):

281 """ 282 Function which calculates the concentration for a given sensor 283 readout using the least square line fitted to the 284 conductance-concentration relation. If a background_readout 285 is given then this is take into consideration. 286 287 @type readout: float or numpy array of floats 288 @param readout: the sensor readout 289 290 @type background_readout: float 291 @param background_readout: the sensor readout at zero salt concentration (default=1) 292 293 @rtype: numpy array or float 294 @return: concentration 295 """ 296 cond = self.readout2conductance(readout) 297 if background_readout is None: 298 background = None 299 else: 300 background = self.readout2conductance(background_readout) 301 return self.conductance2concentration(cond, background)

302

303 - def conductance2concentration(self, conductance, background=None):

304 """ 305 """ 306 if background is None: 307 return np.polyval([self.md.ar, self.md.br], conductance) 308 else: 309 return np.polyval([self.md.ar, 0], conductance - background)

310

311 - def make_calibration(self):

312 """ 313 Makes the calibaration, afterwards the the method 314 L{readout2concentration} is available to be used. 315 316 @rtype: dict: {function, float, float, float, float} 317 @return: A dictionary containing: 318 - 'a': a in ax+b 319 - 'b': b in ax+b 320 - 'R_squared': R^2 321 - 'rmse': root mean squared error' 322 """ 323 324 # check that the right units are used 325 if ( (not self.md.calibaration_solution_concentration_units=='g/l') 326 and (not self.md.calibaration_solution_concentration_units=='gl^-1')): 327 raise TypeError("self.md.calibaration_solution_concentration_units is not 'g/l' or 'gl^-1' but '%s'" 328 % self.md.calibaration_solution_concentration_units) 329 if not self.md.bucket_size_units=='l': 330 raise TypeError("self.md.bucket_size_units is not 'l' but '%s'" 331 % self.md.bucket_size_units) 332 if not self.md.units[1]=='ml': 333 raise TypeError("self.data['calibaration solution'] (i.e. self.md.units[1]) is not 'ml' but '%s'" 334 % self.md.units[1]) 335 336 # calcualte concentration in bucket 337 amount_of_salt = self.data['calibration_solution']/1e3*self.md.calibaration_solution_concentration/1e3 # amount of salt in bucket in [kg] 338 bucket_size_m3 = self.md.bucket_size/1e3 # bucket size in m^3 339 concentration = amount_of_salt/bucket_size_m3 # concentration in kg/m^3 340 self.data = mlab.rec_append_fields(self.data, 'concentration', concentration) 341 342 # convert sensor reading (which is a resistance ration) into 343 # conductance with the formula cond = 1/(Rf * (reading - 1)) 344 sensor_readout = self.data['sensor_readout'] 345 conductance_in_sensor = self.readout2conductance(sensor_readout) 346 self.data = mlab.rec_append_fields(self.data, 'conductance_in_sensor', conductance_in_sensor) 347 348 # fit a straight line by least squares 349 n_obs = len(concentration) 350 (ar, br) = np.polyfit(conductance_in_sensor, concentration, 1) 351 # ar,residuals, rank, singular_values, rcond = np.polyfit(conductance_in_sensor, concentration, 1, full=True) 352 # print ar, residuals, rank, singular_values, rcond 353 # br = ar[1] 354 # ar = ar[0] 355 predicted_concentration = np.polyval([ar,br], conductance_in_sensor) 356 #compute the squared error 357 mean_squared_error = sum((predicted_concentration-concentration)**2) 358 # root mean squared error 359 rmse = np.sqrt(mean_squared_error)/n_obs 360 # coefficient of determination R^2 361 mean_concentration = sum(concentration)/n_obs 362 Stot = sum((concentration - mean_concentration)**2) 363 Serr = sum((predicted_concentration - concentration)**2) 364 R_squared = 1 - Serr/Stot 365 print('Linear regression of conductance vs concentration') 366 print('a=%.2f b=%.2f, root mean squared error= %.3f, R^2 = %.3f' % (ar, br, rmse, R_squared)) 367 368 self.md.ar = ar 369 self.md.br = br 370 self.md.R_squared = R_squared 371 self.md.rmse = rmse 372 373 374 # plot it 375 if self._plot_yes: 376 self.plot_it()

377

378 - def plot_it(self, axs=None):

379 """ 380 Plot the calibration curve. 381 382 @type axs: list of matplotlib axes object 383 @param axs: two matplotlib axes in which to plot (default None: make a new ones). 384 385 @rtype: list of matplotlib axes object 386 @returns: list of the used matplotlib axes 387 """ 388 # calc all the needed quantities 389 conductance_in_sensor = self.data['conductance_in_sensor'] 390 sensor_readout = self.data['sensor_readout'] 391 concentration = self.data['concentration'] 392 predicted_concentration = np.polyval([self.md.ar,self.md.br], conductance_in_sensor) 393 sr = np.linspace(sensor_readout[-1]-0.05,sensor_readout[0]+0.05, 50) 394 fnval = self.readout2concentration(sr) 395 396 # plot linear regression 397 if axs is None: 398 axs = [] 399 axs.append(plt.subplot(2,1,1)) 400 axs.append(plt.subplot(2,1,2)) 401 axs[0].yaxis.grid(True) 402 axs[0].set_title('Linear Regression between conductance and concentration: a=%.2f b=%.2f, root mean squared error= %.3f, R2 = %.3f' 403 % (self.md.ar, self.md.br, self.md.rmse, self.md.R_squared)) 404 axs[0].plot(conductance_in_sensor, predicted_concentration,'r') 405 axs[0].plot(conductance_in_sensor, concentration,'kx') 406 if np.isnan(self.md.reference_resistor): 407 axs[0].set_xlabel('Conductance (arbitrary units)') 408 else: 409 axs[0].set_xlabel('Conductance (1/%s)' % self.md.reference_resistor_units) 410 axs[0].set_ylabel('Concentration (kg/m3)') 411 axs[0].legend(['regression','data'], loc='best') 412 413 # plot sensor vs concentration 414 axs[1].yaxis.grid(True) 415 axs[1].set_title('Sensor readout vs concentration') 416 axs[1].plot(sr, fnval,'r') 417 axs[1].plot(sensor_readout, predicted_concentration,'.r') 418 axs[1].plot(sensor_readout, concentration,'kx') 419 axs[1].legend(['regression','data'], loc='best') 420 axs[1].set_xlabel('Sensor readout ()') 421 axs[1].set_ylabel('Concentration (kg/m3)') 422 423 return axs

424

Source Code for Module fieldpy.stream_gauging.calibration