fieldpy.stream_gauging.dilution

1 #!/usr/bin/env python2 2 """ 3 This file contains the necessary stuff to calculate the stage 4 discharge relationship for a particular day when dilution experiments 5 were done. 6 7 Workflow: 8 - make input files: 9 - dilution_exps*.maw: the file specifying all the used input files 10 - calibration_*.maw: file with the calibration information 11 - injections_*.maw: file with the injection information 12 - run make_stage_discharge_relation: 13 - check you're happy with the calibrations which will pop-up as 14 plots 15 - check you're happy with the stage-discharge relation which will 16 pop-up as plots 17 - use the returned stage-discharge function to turn the stage 18 timeseries into discharge timeseries 19 """ 20 import copy 21 22 import matplotlib.mlab as mlab 23 import pylab as plt 24 25 from fieldpy.core.experiment_classes import * 26

27 -class SaltInjection(LoggedData):

28 """ 29 Class to hold a single salt injection. 30 """

31 - def __init__(self, inj_file, inj_number, time_interval, 32 conductivity, conductivity_field, 33 all_calib, stage, stage_field, stage_averaging_time=3600):

34 """ 35 After initialisation execute self.post_init to do: 36 - self.make_concentration 37 - self.calc_discharge 38 - self.get_stage 39 40 @type inj_file: sting 41 @param inj_file: A file containing information about salt 42 injections, formatted like 43 test_data/injections_13JUL10.maw 44 45 @type inj_number: int 46 @param inj_number: The index of the injection wrt above file 47 48 @type time_interval: tuple of floats or None 49 50 @param time_interval: The time interval (measured in secs 51 after injection) used for the calcuation 52 of discharge and such. If set to None 53 then (0, 240) is used. 54 55 @type conductivity: L{conductivity.Conductivity} 56 @param conductivity: An instance holding the conductivity data 57 spanning the experiment time. 58 59 @type conductivity_field: string 60 @param conductivity_field: The field of conductivity.data to use 61 62 63 @type all_calib: L{calibration.AllCalibrations} 64 @param all_calib: Instance holding the calibration data and 65 associated conductivity to 66 67 @type stage: L{fieldpy.stream_gauging.stream_gauging.Stage} 68 @param stage: Instance holding the stage 69 70 @type stage_field: string 71 @param stage_field: The field of stage.data to use 72 73 @type stage_averaging_time: float 74 @param stage_averaging_time: the time (in secs) over which to 75 average the stage. (default 3600) 76 """ 77 super(SaltInjection, self).__init__() 78 79 self.filename = inj_file 80 self.inj_data, self.raw_data, self.md = raw_file_readers.read_maw_file(self.filename) 81 self.md.inj_number = inj_number 82 self.md.inj_time = self.inj_data['time'][inj_number] 83 self.md.inj_amount = self.inj_data['amount'][inj_number] 84 85 86 if time_interval is None: 87 self.md.time_interval = (0., 240.) 88 else: 89 self.md.time_interval = time_interval 90 91 #: the calcualted discharge 92 self.md.discharge = None 93 #: the averaged stage 94 self.md.stage = None 95 self.md.stage_field = stage_field 96 self.md.stage_averaging_time = stage_averaging_time 97 98 self._datetime = plt.num2date(self.md.inj_time) 99 100 self.md.conductivity_field = conductivity_field 101 self._conductivity_instance = conductivity 102 self._all_calib_instance = all_calib 103 self._stage = stage 104 105 # fill the self.data with a cut out piece of the conductivity data 106 self.data = copy.deepcopy(self._conductivity_instance.data) 107 self.data = helper_fns.mrec_keep_fields(self.data, ['time', self.md.conductivity_field]) 108 self.cut_out_ec_data() 109 110 self.check()

111

112 - def post_init(self):

113 """ 114 Run after initialisation. 115 116 @todo: figure out a way to handle the make_concentration arguments... 117 @todo: figure out a way to handle the get_stage arguments... 118 """ 119 self.make_concentration() 120 self.calc_discharge() 121 self.get_stage()

122

123 - def check(self):

124 """ 125 Does integrity checks: 126 - check the value for the reference_resistor are the same in 127 the conductivity instance and all_calib instance 128 """ 129 super(SaltInjection, self).check() 130 ref_r1 = self._conductivity_instance.md.reference_resistor 131 ref_r2 = self._all_calib_instance.calibrations[0].md.reference_resistor 132 if ref_r1!=ref_r2: 133 raise ValueError('The reference_resistor used in the Conductivity instance is not equal with the one in the AllCalibrations instance.')

134 - def cut_out_ec_data(self):

135 """ 136 Cuts out a piece of the filtered conductivity data lp_ec in 137 self._conductivity_instance. Fills self.data with 138 - time 139 - lp_ec: low pass filtered electict conductivity 140 141 Uses self.md.time_interval for the cut interval 142 143 @todo: adjust here once we got units implemented 144 """ 145 conversion_fac = 1./24./3600. 146 start_time = self.md.inj_time + self.md.time_interval[0]*conversion_fac 147 end_time = self.md.inj_time + sum(self.md.time_interval) *conversion_fac 148 self.cut_time_series([start_time, end_time])

149

150 - def make_concentration(self, bckg_aver_window=4, verbose_=True, 151 stage2dis_method='closest'):

152 """ 153 Appends a field 'conc' to self.data containing the salt 154 concentration derived from the L{calibration.AllCalibrations} 155 instance. 156 157 @type bckg_aver_window: int 158 @param bckg_aver_window: The size of the window (in samples) 159 over which to calulate an average of 160 the conductivity background. 161 162 @note: 163 - there is no provisions if the background at the 164 beginning of the trace differs much from at the end. Though, 165 this should not happen as the salt cloud passes quickly. 166 - there should be a way to choose which conductance2concentration 167 function to use 168 """ 169 self.md.background1 = sum(self.data[self.md.conductivity_field][0:bckg_aver_window])/bckg_aver_window 170 self.md.background2 = sum(self.data[self.md.conductivity_field][-bckg_aver_window:])/bckg_aver_window 171 self.md.background = 0.5*(self.md.background1+self.md.background2) 172 if verbose_: 173 print 'Backgrounds: ', self.md.background1, self.md.background2, self.md.background 174 if stage2dis_method=='closest': 175 conc = self._all_calib_instance.conductance2concentration_closest(self.data['time'], 176 self.data[self.md.conductivity_field], self.md.background, verbose_=True) 177 elif stage2dis_method=='average': 178 conc = self._all_calib_instance.conductance2concentration_average(self.data['time'], 179 self.data[self.md.conductivity_field], self.md.background, verbose_=True) 180 else: 181 raise TypeError('Unknown stage2dis_method') 182 183 184 try: 185 self.data = helper_fns.append_a2mrec(self.data, conc, 'concentration') 186 except: 187 self.data['concentration'] = conc

188

189 - def calc_discharge(self):

190 """ 191 Calculates the discharge by dividing the injected mass by the 192 integrated concentration. Stores it in self.md.discharge 193 """ 194 conversion_factor = 24*3600 195 self.md.discharge = self.md.inj_amount / (self.integrate('concentration') *conversion_factor)

196

197 - def get_stage(self):

198 """ 199 Averages the stage over averaging_time window (in hours) and 200 puts it into self.md.average_stage and self.md.stage_averaging_time. 201 """ 202 averaging_time = self.md.stage_averaging_time 203 dt2 = averaging_time/2./24/3600. 204 205 sli = self._stage.get_ind_as_slice([self.md.inj_time-dt2, self.md.inj_time+dt2]) 206 self.md.average_stage = self._stage.data[self.md.stage_field][sli].mean()

207

208 - def plot_breakthrough():

209 pass

210 211

212 -class DilutionGauging(ExperimentTimeSeries):

213 """ 214 Class for the dilution gauging experiments of one day, with the 215 aim of making a stage discharge relation for that day. 216 """

217 - def __init__(self, inj_file, time_intervals, conductivity, conductivity_field, 218 all_calib, stage, stage_field, stage_averaging_time=3600):

219 """ 220 After initialistaion exacute: 221 - self.make_concentration 222 - self.calc_discharge 223 - self.get_stage 224 225 @type inj_file: sting or None 226 @param inj_file: A file containing information about salt 227 injections, formatted like 228 test_data/injections_13JUL10.maw 229 If None an essentially empty instance will be created. 230 231 @type time_intervals: (list of) tuple of floats or None 232 @param time_intervals: The list of time intervals (measured in secs 233 after injection) used for the calcuation 234 of discharge and such. If set to None 235 then (0, 240) is used. 236 237 @type conductivity: L{fieldpy.stream_gauging.conductivity.Conductivity} 238 @param conductivity: An instance holding the conductivtiy data 239 spanning the experiment time. 240 241 @type all_calib: L{fieldpy.stream_gauging.calibration.AllCalibrations} 242 @param all_calib: Instance holding the calibration data and associated 243 conductivity to 244 245 @type stage: L{fieldpy.stream_gauging.stream_gauging.Stage} 246 @param stage: Instance holding the stage 247 248 @type stage_field: string 249 @param stage_field: The field of stage.data to use 250 251 @type stage_averaging_time: float 252 @param stage_averaging_time: the time (in secs) over which to 253 average the stage. (default 3600) 254 """ 255 super(DilutionGauging, self).__init__() 256 if inj_file is None: 257 # do nothing 258 return 259 self.filename = inj_file 260 self.data, self.raw_data, self.md = raw_file_readers.read_maw_file(self.filename) 261 262 self._date = plt.num2date(np.round(self.data['time'][0])) 263 #: number of injections 264 self.md.n_inj = len(self.data) 265 if time_intervals is None: 266 time_intervals = [(0.,240.) for ii in range(self.md.n_inj)] 267 elif type(time_intervals)==type(()): 268 time_intervals = [time_intervals for ii in range(self.md.n_inj)] 269 270 # make all salt injections 271 self.injections = [] 272 dis = [] 273 stages = [] 274 for ii in range(self.md.n_inj): 275 self.injections.append( 276 SaltInjection(self.filename, ii, time_intervals[ii], 277 conductivity, conductivity_field, all_calib, 278 stage, stage_field, stage_averaging_time)) 279 self.injections[-1].post_init() 280 dis.append(self.injections[-1].md.discharge) 281 stages.append(self.injections[-1].md.average_stage) 282 283 # append the stage and discharge to self.data 284 self.data = mlab.rec_append_fields(self.data, 'discharge', np.array(dis)) 285 self.data = mlab.rec_append_fields(self.data, 'stage', np.array(stages)) 286 287 self.make_stage_discharge_relation() 288 self.check()

289

290 - def check(self):

291 """ 292 Does integrity checks: 293 - nothing so far 294 """ 295 super(DilutionGauging, self).check()

296

297 - def make_stage_discharge_relation(self, order_of_poly=2, verbose_=False):

298 """ 299 Makes the stage discharge relation by fitting a polynomial to 300 the stage and discharge record. 301 302 @type order_of_poly: int 303 @param order_of_poly: the order of the polynomial 304 305 @param verbose_: if True 306 """ 307 self.md.stage_dis_coef,residuals, rank, singular_values, rcond \ 308 = np.polyfit(self.data['stage'], self.data['discharge'], 309 order_of_poly, full=True) 310 if verbose_: 311 print 'Fitting a stage-discharge relation\n Coef: %s\n residuals: %s\n rand: %s\n singular values: %s\n condition number: %s' \ 312 % (str(self.md.stage_dis_coef), str(residuals), str(rank), str(singular_values), str(rcond)) 313 ax = self.plot_stage_discharge_relation() 314 ax.hold(True) 315 self.plot_stage_discharge_measurements(ax)

316

317 - def plot_stage_discharge_relation(self, ax=None, color='b'):

318 """ 319 Does as advertised. 320 """ 321 if ax is None: 322 fig = plt.figure() 323 ax = fig.add_subplot(111) 324 ax.plot(self.data['stage'], self.data['discharge'], color+'x', label='Measured') 325 ax.set_xlabel('Stage (m)') 326 ax.set_ylabel('Discharge (m^3/s)') 327 ax.legend() 328 return ax

329

330 - def plot_stage_discharge_measurements(self, ax=None, color='b'):

331 if ax is None: 332 fig = plt.figure() 333 ax = fig.add_subplot(111) 334 stage = np.linspace(np.min(self.data['stage']), np.max(self.data['stage']), 30) 335 dis = self.stage2dis(stage) 336 ax.plot(stage, dis, '-'+color, label='Fit') 337 ax.set_xlabel('Stage (m)') 338 ax.set_ylabel('Discharge (m^3/s)') 339 ax.legend() 340 return ax

341

342 - def stage2dis(self, stage):

343 """ 344 Converts an array of stage data into discharge. 345 346 @param stage: array of stages 347 348 @return: discharge 349 """ 350 return helper_fns.polyval(self.md.stage_dis_coef, stage)

351

352 - def return_stage2dis_fn(self):

353 """ 354 Returns a stage to discharge function which can be used 355 without the DilutionGauging instance. 356 """ 357 coefs = self.md.stage_dis_coef 358 return lambda stage: helper_fns.polyval(coefs, stage)

359

360 - def plot_all(self):

361 """ 362 Plot the breakthrough curves of all injections. 363 """ 364 plt.figure() 365 for ii,inj in enumerate(self.injections): 366 ax = plt.subplot(self.md.n_inj, 1, ii+1) 367 title = self.raw_data['time'][ii-1] 368 inj.plot_ts('concentration', ax,label=title) 369 plt.legend()

370 #to plot the time from injection start 371 #plt.legend((title)) 372 373

374 -class Stage2Discharge(ExperimentData):

375 """ 376 Class for holding all the dilution gauging experiments of one 377 season, with the aim of providing stage discharge relations for a 378 season. 379 380 Similar in intent to L{calibration.AllCalibrations} in the sense 381 that it unifies several L{dilution_gauging.DilutionGauging} 382 instances. 383 384 It will provide several stage to discharge functions: 385 - L{stage2dis_all_fn} 386 """

387 - def __init__(self, list_of_dilutiongauging_inst, plot_yes=False):

388 """Initialises the class with 389 390 @type list_of_dilutiongauging_inst: list of L{dilution_gauging.DilutionGauging} 391 instances 392 393 @param list_of_dilutiongauging_inst: A list of L{dilution_gauging.DilutionGauging} 394 instances which will be used to make an overall 395 stage2discharge relation 396 397 @type plot_yes: boolean 398 @param plot_yes: Plot the stage-discharge relationships if True (default = False) 399 400 @todo: specify order of polyon for stage2dis 401 """ 402 super(Stage2Discharge, self).__init__() 403 self.dilution_gaugings = list_of_dilutiongauging_inst 404 self._plot_yes = plot_yes 405 # sort them according to when they were done 406 cmpfn = lambda x, y : cmp(x._date, y._date) 407 self.dilution_gaugings.sort(cmp=cmpfn) 408 self.check()

409

410 - def check(self):

411 """ 412 No checks so far 413 """ 414 super(Stage2Discharge, self).check()

415

416 - def plot_comparsion(self):

417 """ 418 plot a comparsion of all stage discharge relation 419 """ 420 dg = self.dilution_gaugings[0] 421 ax = dg.plot_stage_discharge_relation() 422 dg.plot_stage_discharge_measurements(ax) 423 colors = ['r', 'g', 'm', 'k'] 424 for ii,dg in enumerate(self.dilution_gaugings[1:]): 425 dg.plot_stage_discharge_relation(ax, colors[ii]) 426 dg.plot_stage_discharge_measurements(ax, colors[ii])

427

428 - def stage2dis_all_fn(self, stage, dilution_gaugings_to_average):

429 """ 430 Returns a function which converts the stage into a discharge 431 using a fit to all stage discharge mesurements. 432 433 @type stage: numpy array or float 434 @param stage: the stage to convert into discharge 435 436 @type dilution_gaugings_to_average: list of int 437 @param dilution_gaugings_to_average: List of which gaugings of 438 self.dilution_gaugings to use 439 440 @rtype: numpy array or float 441 @return: discharge 442 """ 443 # average the self.calibrations.md.ar coefficient 444 c2a = calibrations_to_average 445 coef_aver = c2a[0] 446 for caa in c2a[1:]: 447 coef_aver += self.calibrations[caa].md.stage_dis_coef 448 coef_aver = coef_aver/len(c2a) 449 # make an empty DilutionGauging instance 450 dil2use = DilutionGauging(None, None, None, None, None) 451 # and fill in the averaged coefficients 452 dil2use.md.stage_dis_coef = coef_aver 453 if verbose_: 454 print 'Used coefficient is' % str(coef_aver) 455 return dil2use.stage2dis(stage 456 )

457

Source Code for Module fieldpy.stream_gauging.dilution_gauging