#!/usr/bin/env python # # This is AOOS_MBS11_convert.py # # This script is written for Python 2.5.2. # It converts the (semi) raw data delivered by # the CR10X data logger at the Mass Balance Site # to meaningful numbers. # It is a rewrite of Daniel Pringle's Matlab script # AOOS_MBS08_convert.m which served the same purpose. # # This script was written by Chris Petrich in January 2009 # and modified and adjusted several times since. # Last modified on May 29, 2011. # from numpy import log, zeros, where, reshape import numpy as np from scipy import select import sys, copy # # TODO: currently (2006--2010), the data logger corrects the # distances measured by the pingers with the ambient temperature. # This is fine if the temperature measurement is accurate # (the water temperature reading seems to be 0.2 C low) and not # outside the range of the sensor (i.e. above -40 C). I think # it would be better to have the data logger do all the logging # and the processing scripts do all the processing. # # Probe modifications in the 2010 season: # * new thermistor strings from CRREL, 2x15 strings # * now dedicated control resistor for each string # # * The Y-connector between upper and lower pinger mast broke # again during installation in 2010 (this already happened # in 2008, with a new Y-connector installed in 2009). ## # Reference point at the pinger mast is the upper rim of # the Y-connector between upper and lower section. In # previous years, the bolt hole was used. # # Distance between sea ice--snow interface and face of # upward looking sounder (UWS1): UWS1_ref = 2.09+0.85 # Distance between sea ice--snow interface and face of # downward looking sounder (UWS2): UWS2_ref = 6.30-2.90 # # SR50, ID=2 # height above snow--ice interface SR50_mast = 1.358+0.04 # # SR50A with ID=0 (located at albedo radiometers) SR50A_ID0 = 1.366+0.04 # after adjustedment around day 140 -- this value # is an estimate (May 29, 2011) SR50A_ID0_adjusted = .65 # # SR50A with ID=1 (located at under/in-ice radiometers) SR50A_ID1 = 1.298+0.04 # # the temperature offset of the underwater thermometer # T107_offset = 0.5 # offset of the air temperature probe: T500_offset = 0 # invalid numers are indicated as NaN_equiv = -9999 # input matrix: MBS # output matrix: OUT fn_in = 'BRW_MBS_2011_raw.txt' fn_out= 'BRW_MBS_2011.txt' try: # input file name: fn_in = sys.argv[1] # output file name: fn_out = sys.argv[2] except: print "expected two parameters: (raw) input file name and output file name" print "using", fn_in, "and", fn_out, ", respectively" # read data logger data # tolerant to incomplete data: want_cols = 46 MBS = None f = open(fn_in,'rb') try: for line in f: line = line.strip() data = line.split(',') if len(data) != want_cols: # incomplete line continue if int(data[0]) != 110: # wrong array ID continue entry = np.nan * np.ones((1,len(data))) for idx in range(len(data)): entry[0,idx] = float(data[idx]) if MBS == None: MBS = copy.copy(entry) else: MBS = np.append(MBS, entry, axis=0) finally: f.close() # dictionary of output headers header={} # if we have only one entry then MBS will be a vector # which is a column array in numpy/python. However, we want a row! # hence, we'll have to reshape the vector into an array # of our liking if len(MBS.shape) == 1: # this is a vector (1D) MBS = reshape(MBS, (1,MBS.shape[0])) # remove entries with wrong year: want_year = 2011 MBS = MBS[MBS[:,1]==want_year,:] # create output array with all-zero entries OUT = zeros(MBS.shape, float) # The mass balance site runs on Alaska standard time (non-daylight savings) # output data is in UTC # column 0: array ID in MBS (e.g. 109), fractional day of year in OUT (UTC) day = MBS[:,2] hour = (MBS[:,3]+.5)//100 minute = MBS[:,3] - hour*100 # do time zone adjustment to UTC: hour = hour + 9 day = where([hour >=24], day+1, day) hour = where([hour >=24], hour-24, hour) #if hour >= 24: # day = day +1 # hour = hour - 24 OUT[:,0] = day + (hour+minute/60.)/24. header[0] = 'day (UTC)' # column 1: year OUT[:,1] = MBS[:,1] header[1] = 'year' # column 2: day (UTC) OUT[:,2] = day header[2] = 'doy' # column 3: hour and minute: hhmm (UTC) OUT[:,3] = hour * 100 + minute header[3] = 'HHMM' # column 4: data logger temperature (deg C) OUT[:,4] = MBS[:,4] header[4] = 'T_CR10X (C)' # column 5: Battery voltage (V) OUT[:,5] = MBS[:,5] header[5] = 'V (V)' # column 6: Water temperature (T107 probe) (deg C) # this temperature is used within the CR10X # to correct the speed of sound. # # DP: "The data logger program modifies the Benthos sounder # output to account for temperature dependence of speed of # sound in water. Underwater sounder range is calculated # using v(sound, water) = 1500 m/s. To correct for salinity, # density and temperature, one can use # (MacKenzie, K.V. JASA vol. 70, No. 3 Sept 1981, 807-812): # c[m/s] = 1449 + 4.6 T[C] + 1.34(S[psu] - 35) + 0.00163 z[m]. # For us, we simplified to c[m/s] = 1449 + 4.6 T[C] # using Twater[C] from T107 probe." OUT[:,6] = MBS[:,6] + T107_offset OUT[:,6] = where(OUT[:,6] > -100., OUT[:,6], NaN_equiv) header[6] = 'T_water (C)' # NOTE: since late Jan 2011, the T107 readings are obviously bogus, # varying by >+-0.2 C instead of <+-0.1C # --> correct Benthos distances by assuming a constant speed of # sound. Benthos_correct = (1449. + 4.6 * (-1.8)) / (1449. + 4.6 * MBS[:,6]) # column 7: relative humidity (CS500 probe) OUT[:,7] = MBS[:,7] header[7] = 'RH (%)' # column 8: air temperature (CS500 probe); underranges for T<=-39.66 C OUT[:,8] = MBS[:,8] + T500_offset header[8] = 'T_air (C)' # columns 9--11: snow depth from SR50(a) # This includes a correction based on air temperature from CS500 # applied within the data logger. # the set-up of pinger 0 (albedo stand) got adjusted during the season OUT[:,9] = where(OUT[:,0] < 140., SR50A_ID0 - MBS[:,9], SR50A_ID0_adjusted - MBS[:,9]) #OUT[:, 9] = SR50A_ID0 - MBS[:, 9] OUT[:,10] = SR50A_ID1 - MBS[:,10] OUT[:,11] = SR50_mast - MBS[:,11] # SID=2 header[ 9] = 'h_snow_0 (m)' header[10] = 'h_snow_1 (m)' header[11] = 'h_snow (m)' # SR50A ID0 and ID1 fell into the snow --> mask that period of time OUT[:, 9] = select([ (OUT[:,0] >= 124.9)*(OUT[:,0] <= 141.73)==0 ], [ OUT[:, 9] ], NaN_equiv) OUT[:,10] = select([ (OUT[:,0] >= 124.9)*(OUT[:,0] <= 139.97)==0 ], [ OUT[:,10] ], NaN_equiv) # column 12: upward looking benthos sounder: ice thickness (m) OUT[:,12] = UWS1_ref - MBS[:,12] * Benthos_correct header[12] = 'h_ice (m)' # column 13: downward looking benthos sounder: water level (m) # measured with reference to (growth season) snow--ice # interface OUT[:,13] = UWS2_ref + MBS[:,13] * Benthos_correct header[13] = 'h_water (m)' # now convert thermistor readings # Voltage measured over thermistor in bridge circuit with 10 kOhm resistor. # this is new as of 2010: # * 2x15 strings instead of 3x10 # * dedicated control resistor at virtual position 16 # * equations are documented in the metafile # * use nominal resistance and voltage # * YSI 44033 (prob the same as before), fitted range: # -35 to +5 in 1deg steps; and values at +10, -40, -45, -50. V0 = 2500. # mV excitation voltage of bridge circuit R_ref = 10000. (A,B,C) =(1.46843e-03, 2.38214e-04, 1.01608e-07) # copy voltages from control bridge verbatim OUT[:,14] = MBS[:,29] # upper string OUT[:,15] = MBS[:,45] # lower string header[14] = 'V_ctrl1 (mV)' header[15] = 'V_ctrl2 (mV)' # upper thermistor string: for i in range(14, 29): R = R_ref / ( (V0/MBS[:,i]) -1. ) OUT[:,i+2] = -273.15+ 1./(A + B * log(R) + C * (log(R))**3) # bottom thermistor string: for i in range(30, 45): R = R_ref / ( (V0/MBS[:,i]) -1. ) OUT[:,i+1] = -273.15+ 1./(A + B * log(R) + C * (log(R))**3) # there's a chance the voltage is close to V0 in which case # R is close to infty and OUT is -273.15 # --> remove obviously unrealistic temperatures for i in range(16, 46): OUT[:,i] = where(OUT[:,i] > -65., OUT[:,i], NaN_equiv) header[i] = 'T_%+.2i (C)' % ((7 - (i-16) )*10) # check validity of output data # specifically, invalidate unreasonable values for pingers # better: base this #def validify_rangers(distance): # if distance < 0: # distance = NaN_equiv # if distance > 10: # distance = NaN_equiv # return distance # vectorize function on a scalar: #vec_validify_rangers = vectorize(validify_rangers) # check all rangers: #OUT[:,9:14] = vec_validify_rangers(OUT[:,9:14]) # distance from pinger is between 0.5 m (minimum required for operation) # and 10 m (changed to 5 m on May 30, 2009 after the sealevel sounder returned # two readings of 8.4 m while the reading really increased from 3.2 m to # 3.3 m) # "out of range" reading is 99.99, "error" reading is 0 OUT[:,9:14] = select([ (MBS[:,9:14] > 0.5) * (MBS[:,9:14] < 5.) ],[ OUT[:,9:14] ], NaN_equiv) # remove air temperature and humidity data before installation # and equilibration of CS500 # humidity OUT[:,7] = select([ OUT[:,0] >= 13+23/60.+15/3600. ], [ OUT[:,7] ], NaN_equiv) # air temperature OUT[:,8] = select([ OUT[:,0] >= 13+23/60.+15/3600. ], [ OUT[:,8] ], NaN_equiv) # CS500 not connected if temperature is bogus # NB does not read lower then -39.66C ! OUT[ OUT[:,8]<-40 ,7] = NaN_equiv OUT[ OUT[:,8]<-40 ,8] = NaN_equiv # remove all data before field installation of the probe: for i in range(4,46): OUT[:,i] = select([ OUT[:,0] > 12. ], [ OUT[:,i] ], NaN_equiv) # make sure we don't have "nan" in the array OUT[:,:] = where(OUT[:,:] == OUT[:,:], OUT[:,:] , NaN_equiv ) # sort output by time (first column) # NB mergesort does preserves the order of entries of identical time OUT=OUT[ OUT[:,0].argsort(kind='mergesort'), ] # insert rows of NaN for each missing entry, # assuming a time step of log_dt = 15 # minutes idx = 0 while idx < len(OUT[:,0])-1: if OUT[idx+1,0] - OUT[idx,0] > (log_dt+.1) / (60.*24.): #if OUT[idx,3] == 0: # there is at least one entry missing after # idx newrow = list() newrow.append(OUT[idx,0]+log_dt/(60.*24)) (year, day, HHMM) = OUT[idx,1:4] HH = HHMM // 100 MM = HHMM - 100* HH MM += int(log_dt) if MM > 59: MM -= 60 HH += 1 if HH > 23: HH -= 24 day += 1 if day > 365: year += 1 day = 1 HHMM = HH*100 + MM newrow.append(year) newrow.append(day) newrow.append(HHMM) for i in range(len(OUT[0,4:])): newrow.append(-9999) newrow = np.array(newrow) newrow = reshape(newrow, (1,newrow.shape[0])) OUT = np.append( np.append(OUT[:(idx+1),], newrow, axis=0), OUT[(idx+1):,], axis=0 ) idx += 1 # sort output by time (first column) # NB mergesort does preserves the order of entries of identical time OUT=OUT[ OUT[:,0].argsort(kind='mergesort'), ] # and make sure entries are unique while True: # keep only one entry for each time try: doubles = np.nonzero(abs(OUT[1:,0] - OUT[:-1,0]) < 1e-8)[0] try: if len(doubles)>1: doubles = doubles[0] except: pass except: # ok, none there. break try: doubles = int(doubles) except: # empty list or None # finished break OUT = np.append( OUT[:doubles,], OUT[(doubles+1):,], axis=0) # output result # with up to 7 significant digits, tab delimited (%.7g) # NB we need 7 digits only to make the date significant to the minute! # everything else is sufficient with 3 decimal places at the most. # ==> output as %.3f --> date precise to 1.5 minutes (15 minutes = # 0.01 day). Temperatures to 1 mK, distances to 1 mm. #savetxt(fn_out, OUT, fmt='%.3f', delimiter='\t') # format each column individually f = open(fn_out,'wb') try: # write header line line = "#" for i in range(len(OUT[0,:])): if i != 0: line +='\t' line += header[i] f.write(line+'\n') # write data for i in range(len(OUT[:,0])): if OUT[i,0] < 12.: # this was testing continue s="" # doy, year, day, time s+="%8.4f\t%4.0f\t%3.0f\t%4.4i\t" % \ (OUT[i,0], OUT[i,1], OUT[i,2], int(OUT[i,3]+.1)) # T_CR10X if OUT[i, 4] > -9000.: s+="%.3f\t" % OUT[i, 4] else: s+="%.0f\t" % OUT[i, 4] # Vbatt if OUT[i, 5] > -9000.: s+="%.2f\t" % OUT[i, 5] else: s+="%.0f\t" % OUT[i, 5] # T107 if OUT[i, 6] > -9000.: s+="%.3f\t" % OUT[i, 6] else: s+="%.0f\t" % OUT[i, 6] # RH if OUT[i, 7] >= 0.: s+="%5.1f\t" % OUT[i, 7] else: s+="%.0f\t" % OUT[i, 7] # T_air if OUT[i, 8] > -9000.: s+="%.3f\t" % OUT[i, 8] else: s+="%.0f\t" % OUT[i, 8] # snow pingers 0, 1, and 2: s+="%g\t%g\t%g\t" % \ (OUT[i,9], OUT[i,10], OUT[i,11]) # underwater sounders s+="%g\t%g\t" % \ (OUT[i,12], OUT[i,13]) # reference voltages s+="%g\t%g" % \ (OUT[i,14], OUT[i,15]) # ref. and thermistor string thermistors for j in range(16, 46): if OUT[i, j] < -9000.: s+="\t%.0f" % OUT[i, j] else: s+="\t%.3f" % OUT[i, j] f.write(s+'\n') finally: f.close()