will apply a notch filter to the array bx by finding the nearest peak around the supplied notch locations. The filter is a zero-phase Chebyshev type 1 bandstop filter with minimal ripples.
combineEdifiles(edifile1,edifile2,read1,read2) will combine edifile1 with edifile2 according to read1 and read2. It will combine frequencies, impedance and tipper.
Note nread1 is from the start for edifile1 and nread2 is from end for edifile2.
edifile1 = full path to edifile with the largest frequencies edifile2 = full path to edifile with smallest frequencies, longest
periods
combineFewFiles(dirpath,station,starttime,endtime,cachelength, complst=[‘BX’,’BY’,’BZ’,’EX’,’EY’],d=0) Will combine files in a directory path (dirpath) that have a given start and end time in the form of (HHMMSS). It looks for files at cachelength then decimates the data by d. Input:
dirpath = directory path to folder where files to combine reside station = station name starttime = start time as 6 character string hhmmss endtime = end time as 6 character string hhmmss cachelength = cacherate of each file as 6 character string hhmmss complst = components to combine d = decimation factor, if no decimation enter as 0 or 1
will combine files from different days into one file.
dirpath = directory path where station files are station = name of station daylst = list of information to combine files with a dictionary for
- each day with keys:
- day = utc day (thee character string) start = start time in (‘hhmmss’) stop = end time in(‘hhmmss’) filt = ‘Filtered’ if time series was filtered
cacherate = length of file (‘hhmmss’) dec = decimation factor (0 for none) complst = list of components to combine fdict = dictionary for adaptiveNotchFilter
convertCounts2Units(filenames,eyn=’n’,lpyn=’n’,egain=1.0,dlgain=1.0,exlen=100.,eylen=100., magtype=’lp’,zadj=2) will convert a set of files given by filenames with parameters defined: filenames, eyn => electric channels converted y or n lpyn => long period magnetic channels converted y or n egain => electric gain dlgain => data logger gain exlen => ex length (m) eylen => ey length (m) magtype => magnetic sensor type lp for longperiod or bb for broadband zadj => bz adjusting parameter for bartington sensor
Convert the electric field from counts to units of microV/m. efield is a list, dlgain is gain applied by data logger(verylow=2.5, low=1, high=.1), egain is interface box gain (10,100),dlength is the dipole length of that component.
convertfiles will convert data of counts from data logger to units.
Convert the magnetic field from counts to units of microV/nT. bfield is a list of numbers. dlain is amount of gain applied by data logger(verylow=2.5,low=1, high=.1)
Will decimate a function by the factor m. First an 8th order Cheybechev type I filter with a cuttoff frequency of .8/m is applied in both directions to minimize any phase distortion and remove any aliasing. If m is greater than 10, decimatef will be called multiple times.
Will apply a sinc filter of width w to the function f by multipling in the frequency domain. Returns filtered function
getStationInfo(stationinfofile,station) returns info for the nominated station from the stationinfofile as a dictionary with key words in the hdrinfo.
imp2resphase(z,zvar,freq) will convert impedances z[4,3,len(freq)] to resistivities (ohm-m) and phases (deg) as well as the errors of each. Note the phase is calculated using arctan2 putting the phase in the correct quadrant. The xy component is placed into the positive quadrant by adding 180 deg. INPUTS:
z=[[zxxreal,zxximag,zxxvar],...]], shape=(4,3,len(freq) OR z=[[[zxxr+1j*zxxi,zxyr+1j*zxyi],[zyx,zyy]]], shape=(len(freq),2,2) freq=frequency array zvar=[[[zxxvar,zyyvar],[zyxvar,zyyvar]]], shape=(len(freq),2,2) ffactor=fudge factor to make sure resistivities are the right order
makeDayFoldersFromObservatoryData will take observatory data and put it into day folders for processing.
startday = UTC start day as a string of length 3 d = decimation factor of collected data sampling frequency to get to
observatory sampling frequency. So if you sampled at 1000 Hz and the observatory data is at 1 sec, d=1000
ndays = number of days, only used if observatory data is in one big file df = sampling frequency of observatory data in Hz complst = components to make into files, can add ‘BZ’ if you like station = station name of observatory data rowskip = number of rows to skip from observatory file, usually the
header ends at row 6.
makeKML will make a kml file for Google Earth to plot station locations relatively quickly
normalizeL2(f) will return the function f normalized by the L2 norm -> f/(sqrt(sum(abs(x_i)^2))).
Open an MT file, convert counts to units and return an 1D-array. Make sure filename includes the entire path. gain (verylow=2.5,low=1, high=.1), egain same as gain, dlength is dipole length in m of EX,EY. magtype is the magnetic measurment type ‘lp’ for long period and ‘bb’ for broadband, zadj is an adjustment for lp instruments on the z channel
padzeros(f) will return a function that is padded with zeros to the next power of 2 for faster processing for fft or to length npad if given.
readStationInfo(stationinfofile) will read in a .txt (tab delimited) or .csv(comma delimited)file that has the following information: station name, latitude, longitude, elevation, date collected, components measured (a number: 4 for ex,ey,bx,by, 5 for ex,ey,bx,by,bz, 6 for ex,ey,bx,by,bz,tp), magnetic type (bb or lp), dipole lengths (m), data logger gain (very low=2.5,low=1,high=.1),interface box gain (10 or 100), bz correction for longperiod data and the bx,by,bz components as they were measured in the field. Use headers: station, lat, lon, elev, date, mcomps, magtype, ex, ey, dlgain, igain, lpbzcor, magori. Returns a list of the dictionaries with found information.
readedi(edifile) will read in an edi file written in a format given by format given by http://www.dias.ie/mtnet/docs/ediformat.txt. Returns: lat,lon,frequency,Z[zreal+i*zimag],Zvar,tipper,tippervar (if applicable)
removePeriodicNoise will take average a window of length noise period and average the signal for as many windows that can fit within the data. This averaged window is convolved with a series of delta functions at each window location to create a noise time series. This is then subtracted from the data to get a ‘noise free’ time series.
dt = time sample rate (s) noiseperiods = a list of estimated periods with a range of values to
look around [[noiseperiod1,df1]...] where df1 is a fraction value find the peak about noiseperiod1 must be less than 1.
cacherate = length of file (min)
removeStaticShift(edifile,stol=.2,dm=1000) will remove static shift by calculating the median of respones of near by stations, within dm. If the ratio of the station response to the median on either side of 1+-stol then the impedance tensor for that electric component will be corrected for static shift.
rewriteedi(edifile) will rewrite an edifile say if it needs to be rotated or distortion removed.
edifile = full path to edifile to be rewritten znew = impedance tensor if a new one has been created zvarnew = errors in impedance tensor if a new one has been created freqnew = new frequency list if one has been created newfile = ‘y’ for yes or ‘n’ for no if you want a new file to be
created.
tipnew = new tipper array tipvarnew = new tipper error array thetar = rotation angle counter clockwise (N=0, E=-90) dr = ‘n’ if no distortion removal and ‘y’ for distortion removal
rotateB(bx,by) will Rotate the magnetic field such that Bx is pointing to magnetic north and by to geomagnetic east. Assumes setup was orthogonal. Returns list of strings that are 8 significant digits.
Rotate the electric field to geographic north, ex & ey are lists, declination is in degrees
sigfigs(numstr,digits=8,fmt=’g’) will return a string with the proper amount of significant digits for the input number, can be str or float.