diff --git a/FilterFunctions.py b/FilterFunctions.py new file mode 100644 index 0000000..44e5f04 --- /dev/null +++ b/FilterFunctions.py @@ -0,0 +1,97 @@ +from scipy.signal import butter,filtfilt +from scipy import signal,fftpack +from scipy.signal import freqz +import numpy as np +import matplotlib.pyplot as plt + +def nextpow2(i): + n = 1 + while n < i: + n *= 2 + return n + +def butter_bandpass(lowcut, highcut, fs, order): + nyq = 0.5 * fs + low = lowcut / nyq + high = highcut / nyq + b, a = butter(order, [low, high], btype='bandpass', analog=False) + return b, a + +def butter_bandpass_filter(data, lowcut, highcut, fs, order): + b, a = butter_bandpass(lowcut, highcut, fs, order=order) + #w, h = freqz(b, a, worN=2000) + #plt.plot((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order) + y = filtfilt(b,a,data) + return y + +def butter_lowpass(highcut, fs, order): + nyq = 0.5 * fs + high = highcut / nyq + b, a = butter(order, high, btype='lowpass', analog=False) + return b, a + +def butter_lowpass_filter(data, highcut, fs, order): + b, a = butter_lowpass(highcut, fs, order=order) + y = filtfilt(b,a,data) + return y + +def butter_highpass(lowcut, fs, order): + nyq = 0.5 * fs + low = lowcut / nyq + b, a = butter(order, low, btype='highpass', analog=False) + return b, a + +def butter_highpass_filter(data, lowcut, fs, order): + b, a = butter_highpass(lowcut, fs, order=order) + y = filtfilt(b,a,data) + return y + +def filtre_hautbas_calcul(var_filtre, data, cut, fs, ordre): + if var_filtre == 1: + datain_filtered = butter_lowpass_filter(data, cut, fs, ordre) + elif var_filtre == 2: + datain_filtered = butter_highpass_filter(data, cut, fs, ordre) + num = len(datain_filtered) + N = nextpow2(num) + A_filt = fftpack.fft(datain_filtered, N) + A_shift_filt = fftpack.fftshift(A_filt) + Xpos_filt = A_shift_filt[N//2:] + return datain_filtered, Xpos_filt, N + +def filtre_passband_calcul(twtt,dataset, lowcut, highcut, fs, ordre, numero_trace, trace_pour_afficher): + data = dataset[numero_trace] + datain_filtered = butter_bandpass_filter(data, lowcut, highcut, fs, ordre) + num=len(datain_filtered) + N = nextpow2(num) + A_filt = fftpack.fft(datain_filtered, N) + A_shift_filt = fftpack.fftshift(A_filt) + Xpos_filt = A_shift_filt[N//2:] + + if numero_trace == trace_pour_afficher : + + b, a = butter_bandpass(lowcut, highcut, fs, ordre) + w, h = freqz(b, a, worN=2000) + + plt.subplot(231) + plt.plot(twtt, data.tolist()[0]) + plt.xlabel("Temps (ns)") + plt.ylabel("Amplitude") + plt.title('Signal') + + plt.subplot(232) + plt.plot((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % ordre) + plt.xlabel("Fréquence (GHz)") + plt.ylabel("Gain (dB)") + plt.title('Filtre de Butterworth') + + plt.subplot(233) + plt.plot(twtt, datain_filtered.tolist()[0]) + plt.xlabel("Temps (ns)") + plt.ylabel("Amplitude") + plt.title('Signal filtré') + + plt.waitforbuttonpress(0) # this will wait for indefinite time + plt.close() + + + return datain_filtered, Xpos_filt, N diff --git a/__main__.py b/__main__.py new file mode 100644 index 0000000..545000f --- /dev/null +++ b/__main__.py @@ -0,0 +1,52 @@ +import sys +import tkinter as tk + +from gprpy.gprpyGUI import GPRPyApp +from gprpy.gprpyCWGUI import GPRPyCWApp + +def main(args=None): + + if len(sys.argv) < 2: + mode = input("Profile [p] or Common Midpoint / WARR [c]?") + else: + mode = sys.argv[1][0] + + + + if mode == 'p': + rightcol=9 + figrowsp=19+1 + + root = tk.Tk() + + for col in range(rightcol): + root.columnconfigure(col, weight=1) + for row in range(figrowsp): + root.rowconfigure(row, weight=1) + + app = GPRPyApp(root) + + root.mainloop() + + elif mode == 'c' or mode == 'w': + rightcol=10 + figrowsp=15+1 + + root = tk.Tk() + + for col in range(rightcol): + root.columnconfigure(col, weight=1) + for row in range(figrowsp): + root.rowconfigure(row, weight=1) + + app = GPRPyCWApp(root) + + root.mainloop() + + else: + print("You need to chose either profile [p] or CMP/WARR [c] mode.") + + + +if __name__ == "__main__": + main() diff --git a/gprIO_DZT.py b/gprIO_DZT.py new file mode 100644 index 0000000..1512620 --- /dev/null +++ b/gprIO_DZT.py @@ -0,0 +1,231 @@ +import struct +import numpy as np +#import re # Regular expressions +from datetime import datetime +import math +from itertools import takewhile + +####from contents.py in readgssi softwware #### + +# a dictionary of standard gssi antenna codes and frequencies + +ANT = { + '100MHz': 100, + '200MHz': 200, + '270MHz': 270, + '350MHz': 350, + '400MHz': 400, + '500MHz': 500, + '800MHz': 800, + '900MHz': 900, + '1600MHz': 1600, + '2000MHz': 2000, + '2300MHz': 2300, + '2600MHz': 2600, + '3200': 'adjustable', + '3200MLF': 'adjustable', + 'gprMa': 'adjustable', # gprMax support + 'GSSI': 'adjustable', # support for issue #11 + 'CUSTOM': 'adjustable', + '120D' : 120, + '400D':400, + '3207': 100, + '3207AP': 100, + '5106': 200, + '5106A': 200, + '50300': 300, + '350': 350, + '350HS': 350, + 'D400HS': 350, + '50270': 270, + '50270S': 270, + 'D50300': 300, + '5103': 400, + '5103A': 400, + '50400': 400, + '50400S': 400, + '800': 800, + 'D50800': 800, + '3101': 900, + '3101A': 900, + '51600': 1600, + '51600S': 1600, + 'SS MINI': 1600, + '62000': 2000, + '62000-003': 2000, + '62300': 2300, + '62300XT': 2300, + '52600': 2600, + '52600S': 2600, +} + +############################ + +def readdzt(filename): + ''' + Reads a GSSI .DZT data file. + + INPUT: + filename data file name including .DZT extension + + OUTPUT: + data data matrix whose columns contain the traces + info dict with information from the header + + Thanks to Ian Nesbitt for pointing out extended headers and + providing the documentation file. + ''' + + # Documentation file is DZT.File.Format.6-14-16.pdf + + info = {} + + fid = open(filename,'rb'); + + + # H is unsigned int 16 (ushort = uint16) + # h is short (int16) + # I is unsigned int 32 (uint = uint32) + # i is int32 + # f is float + + # All of the following information is from DZT.File.Format.6-14-16.pdf + # Provided by Ian Nesbitt + + minheadsize = 1024 + infoareasize = 128 + + info['name'] = fid.name + + rh_tag = struct.unpack('h', fid.read(2))[0] # Pos 00 + + # Size of the header + rh_data = struct.unpack('h', fid.read(2))[0] # Pos 02 + + # Samples per trace + rh_nsamp = struct.unpack('h', fid.read(2))[0] # Pos 04 + info["Samples per traces :"] = rh_nsamp + + # Bits per word + rh_bits = struct.unpack('h', fid.read(2))[0] # Pos 06 + + # Binary offset + rh_zero = struct.unpack('h', fid.read(2))[0] # Pos 08 + + # Scans per second + rhf_sps = struct.unpack('f', fid.read(4))[0] # Pos 10 + info["Scans per second : "] = rhf_sps + + # Scans per meter + rhf_spm = struct.unpack('f', fid.read(4))[0] # Pos 14 + info["Scans per meter : "] = rhf_spm + + # Meters per mark + rhf_mpm = struct.unpack('f', fid.read(4))[0] # Pos 18 + info['Meters per mark : '] = rhf_mpm + + # Startposition [ns] + rhf_position = struct.unpack('f', fid.read(4))[0] # Pos 22 + info["Startposition (ns) : "] = rhf_position + + # length of trace [ns] + rhf_range = struct.unpack('f', fid.read(4))[0] # Pos 26 + info["Length of trace (ns) : "] = rhf_range + + # frequence d echantillonage [GHs] + info["Sample frequency (GHz) : "]= rh_nsamp/rhf_range + + # Number of passes + rh_npass = struct.unpack('h', fid.read(2))[0] # Pos 30 + info["Number of passes : "] = rh_npass + + # Creation date and time + rhb_cdt = struct.unpack('f', fid.read(4))[0] # Pos 32 + #info["Creation date and time : "] = rhb_cdt + + # Last modified date & time + rhb_mdt = struct.unpack('f', fid.read(4))[0] # Pos 36 + #info["Last modified date and time : "] = rhb_mdt + + # no idea + rh_mapOffset = struct.unpack('h', fid.read(2))[0] # Pos 40 + #info['rh_mapOffset'] = rh_mapOffset + # No idea + rh_mapSize = struct.unpack('h',fid.read(2))[0] # Pos 42 + #info['rh_mapSize'] = rh_mapSize + + # offset to text + rh_text = struct.unpack('h',fid.read(2))[0] # Pos 44 + + # Size of text + rh_ntext = struct.unpack('h',fid.read(2))[0] # Pos 46 + + # offset to processing history + rh_proc = struct.unpack('h',fid.read(2))[0] # Pos 48 + + # size of processing history + rh_nproc = struct.unpack('h',fid.read(2))[0] # Pos 50 + + # number of channels + rh_nchan = struct.unpack('h',fid.read(2))[0] # Pos 52 + info['Number of channels : ']=rh_nchan + # ... and more stuff we don't really need + + #inspired by readgssi software + fid.seek(98) # start of antenna section + rh_ant = fid.read(14).decode('utf-8').split('\x00')[0] + info['Antenna : ']= rh_ant.rsplit('x')[0] + + # info['Antenna frequency (MHz) : '] = ANT[info['Antenna : ']] + try: + info['Antenna frequency (MHz) : '] = ANT[info['Antenna : ']] + + except KeyError: + print('Antenna frequency couln t be read') + + fid.close() + + # offset will tell us how long the header is in total + # There could be a minimal header of 1024 bits + # (very old files may have had 512 bits) + if rh_data < minheadsize: + offset = minheadsize*rh_data + else: + offset = minheadsize*rh_nchan + + # Define data type based on words per bit + # From documentation: Eight byte and sixteen byte samples are + # unsigned integers. Thirty-two bit samples are signed integers. + if rh_bits == 8: + datatype = 'uint8' # unsigned char + elif rh_bits == 16: + datatype = 'uint16' # unsigned int + elif rh_bits == 32: + datatype = 'int32' + + + # Read the entire file + vec = np.fromfile(filename,dtype=datatype) + + headlength = offset/(rh_bits/8) + + # Only use the data, discard the header + datvec = vec[int(headlength):] + + # Turn unsigned integers into signed integers + # Only necessary where unsigned + if rh_bits == 8 or rh_bits == 16: + datvec = datvec - (2**rh_bits)/2.0 + + # reshape into matrix + data = np.reshape(datvec,[int(len(datvec)/rh_nsamp),rh_nsamp]) + + data = np.asmatrix(data) + data_transposed = data.transpose() + data_transposed[0,:]=0 + data_transposed[1,:]=0 + print('data',data_transposed) + #print('min, max', np.min(data_transposed),np.max(data_transposed)) + + + return data_transposed, info \ No newline at end of file diff --git a/gprpy.py b/gprpy.py new file mode 100644 index 0000000..b9e33a7 --- /dev/null +++ b/gprpy.py @@ -0,0 +1,1619 @@ +import sys +sys.path.insert(0, "/home/aline/GPRPY_test/gprpy") +import os +import numpy as np + +import pickle +import gprpy.toolbox.gprIO_DT1 as gprIO_DT1 +import gprpy.toolbox.gprIO_DZT as gprIO_DZT +import gprpy.toolbox.gprIO_BSQ as gprIO_BSQ +import gprpy.toolbox.gprIO_MALA as gprIO_MALA +import gprpy.toolbox.gprpyTools as tools +import gprpy.toolbox.readmrq as readmrq +#from onclick_event import onclick +#from scipy.fftpack import rfft, irfft, fftfreq +from scipy import signal,fftpack +import matplotlib.pyplot as plt +import gprpy.toolbox.FilterFunctions as filters + +try: + import gprpy.irlib.external.mig_fk as mig_fk +except: + print("No fk migration in public version") +import copy +import scipy.interpolate as interp +from pyevtk.hl import gridToVTK + + +class gprpyProfile: + ''' + Ground penetrating radar data processing and visualization class + for common-offset profiles. + ''' + + def __init__(self,filename=None): + ''' + Initialization for a gprpyProfile object. Initialization can be + empty or with a provided filename for the GPR data. + + INPUT: + filename data file name. Currently supported formats: + .gpr (GPRPy), .DT1 (SnS), .DZT (GSSI), .rd3 (MALA), + and ENVI standard BSQ. + ''' + + self.history = ["mygpr = gp.gprpyProfile()"] + + # Initialize previous for undo + self.previous = {} + + # # Initialisation for onclick + # self.frequencies_ = [] + + if filename is not None: + self.importdata(filename) + + def importdata(self,filename, marker = None): + ''' + Loads .gpr (native GPRPy), .DT1 (Sensors and Software), + .DZT (GSSI), .GPRhdr (ENVI standard BSQ), .rad (MALA) + data files and populates all the gprpyProfile fields. + + INPUT: + filename name of the .gpr, .DT1, dt1, .DZT, .GPRhdr, dat, + rd3, or .rad file you want to import. + The header file name and the data file name + have to be the same! + ''' + # Initialisation for onclick + self.frequencies_ = [] + self.t0_ = 0 + + + file_name, file_ext = os.path.splitext(filename) + + if file_ext==".DT1" or file_ext==".HD" or file_ext==".dt1" or file_ext==".hd": + if file_ext==".DT1" or file_ext==".HD": + self.data=gprIO_DT1.readdt1(file_name + ".DT1") + self.info=gprIO_DT1.readdt1Header(file_name + ".HD") + else: + self.data=gprIO_DT1.readdt1(file_name + ".dt1") + self.info=gprIO_DT1.readdt1Header(file_name + ".hd") + + self.profilePos = np.linspace(self.info["Start_pos"], + self.info["Final_pos"], + self.info["N_traces"]) + + self.twtt = np.linspace(self.info["TZ_at_pt"], + self.info["Total_time_window"], + self.info["N_pts_per_trace"]) + + self.antsep = self.info["Antenna_sep"] # Set to m in the loading routine + self.velocity = None + self.depth = None + self.maxTopo = None + self.minTopo = None + self.threeD = None + self.data_pretopo = None + self.twtt_pretopo = None + #frequence d'echantillonage en Hz + self.freq_ech = self.info["N_pts_per_trace"]/(self.info["Total_time_window"]-self.info["TZ_at_pt"])*10E9 + # Initialize previous + self.initPrevious() + + # Put what you did in history + histstr = "mygpr.importdata('%s')" %(filename) + self.history.append(histstr) + + + elif file_ext==".DZT": + + self.data, self.info = gprIO_DZT.readdzt(filename) + + + if marker: + # self.profilePos = self.info["Startposition (ns) : "]+np.linspace(0.0, + # self.data.shape[1]/self.info["Scans per second : "], + # self.data.shape[1]) + #marker = input('marker path and file name : ') + self.profilePos = readmrq.read_mrk(marker, self.data.shape[1]) + print(self.profilePos) + + elif self.info["Scans per meter : "] != 0: + self.profilePos = np.linspace(0.0, + self.data.shape[1]/self.info["Scans per meter : "], + self.data.shape[1]) + print('datashape[1]', self.data.shape[1]) + + else : + print("PROFIL LENGHT UNKNOWN, PLEASE CREATE A MARKER FILE") + + self.twtt = np.linspace(0,self.info["Length of trace (ns) : "],self.info["Samples per traces :"]) + + self.antsep = 0 + self.velocity = None + self.depth = None + self.maxTopo = None + self.minTopo = None + self.threeD = None + self.data_pretopo = None + self.twtt_pretopo = None + + #sample frequency in GHz + self.freq_ech = self.info["Sample frequency (GHz) : "] + + # Initialize previous + self.initPrevious() + + # Put what you did in history + histstr = "mygpr.importdata('%s')" %(filename) + self.history.append(histstr) + + #### @author: aline + #read info for readdzt + for i in self.info: + print(i, self.info[i]) + print ('Number of traces : ', self.data.shape[1]) + #### + + + elif file_ext==".GPRhdr" or file_ext==".dat": + # ENVI standard BSQ file + self.data, self.info = gprIO_BSQ.readBSQ(file_name) + + self.profilePos = float(self.info["dx"])*np.arange(0,int(self.info["columns"])) + self.twtt = np.linspace(0,float(self.info["time_window"]),int(self.info["lines"])) + + self.antsep = 0 + self.velocity = None + self.depth = None + self.maxTopo = None + self.minTopo = None + self.threeD = None + self.data_pretopo = None + self.twtt_pretopo = None + #frequence d'echantillonage en Hz + self.freq_ech = int(self.info["lines"])/float(self.info["time_window"])*10E9 + # Initialize previous + self.initPrevious() + + # Put what you did in history + histstr = "mygpr.importdata('%s')" %(filename) + self.history.append(histstr) + + + elif file_ext==".rad" or file_ext==".rd3" or file_ext==".rd7": + self.data, self.info = gprIO_MALA.readMALA(file_name) + + self.twtt = np.linspace(0,float(self.info["TIMEWINDOW"]),int(self.info["SAMPLES"])) + self.profilePos = float(self.info["DISTANCE INTERVAL"])*np.arange(0,self.data.shape[1]) + + self.antsep = self.info["ANTENNA SEPARATION"] + self.velocity = None + self.depth = None + self.maxTopo = None + self.minTopo = None + self.threeD = None + self.data_pretopo = None + self.twtt_pretopo = None + #frequence d'echantillonage en Hz + self.freq_ech = int(self.info["SAMPLES"])/float(self.info["TIMEWINDOW"])*10E9 + # Initialize previous + self.initPrevious() + + # Put what you did in history + histstr = "mygpr.importdata('%s')" %(filename) + self.history.append(histstr) + + + + elif file_ext==".gpr": + ## Getting back the objects: + with open(filename, 'rb') as f: + data, info, profilePos, twtt, history, antsep, velocity, depth, maxTopo, minTopo, threeD, data_pretopo, twtt_pretopo = pickle.load(f) + self.data = data + self.info = info + self.profilePos = profilePos + self.twtt = twtt + self.history = history + self.antsep = antsep + self.velocity = velocity + self.depth = depth + self.maxTopo = maxTopo + self.minTopo = minTopo + self.threeD = threeD + self.data_pretopo = data_pretopo + self.twtt_pretopo = twtt_pretopo + + # Initialize previous + self.initPrevious() + + else: + print("Can only read dt1, DT1, hd, HD, DZT, dat, GPRhdr, rad, rd3, rd7, and gpr files") + + #sauvegarde de la liste des tps et des données pour pick frequencies (graphe) + self.twtt_before_set_zero = self.twtt + self.data_before_set_zero = self.data + + + + def give_info(self): + ''' + @author: aline + prints info read in file header, only for DZT format + Remplaced in the 'import data' function + Parameters + ---------- + + Returns : + ------- + None. + ''' + + for i in self.info: + print(i, self.info[i]) + print ('Number of traces : ', self.data.shape[1]) + + + + def showHistory(self): + ''' + Prints out processing and visualization history of a data set. + ''' + for i in range(0,len(self.history)): + print(self.history[i]) + + def writeHistory(self,outfilename="myhistory.py"): + ''' + Turns the processing and visualization history into a Python script. + The full path names are saved in the Python script. You can edit the + Python script after saving to remove the full path names. + + INPUT: + outfilename filename for Python script + ''' + with open(outfilename,"w") as outfile: + outfile.write("# Automatically generated by GPRPy\nimport gprpy.gprpy as gp\n") + for i in range(0,len(self.history)): + outfile.write(self.history[i] + "\n") + + def undo(self): + ''' + Undoes the last processing step and removes that step fromt he history. + ''' + self.frequencies_ = self.previous["frequencies"] + self.t0_ = self.previous["time_zero"] + self.data = self.previous["data"] + self.twtt = self.previous["twtt"] + self.info = self.previous["info"] + self.profilePos = self.previous["profilePos"] + self.velocity = self.previous["velocity"] + self.depth = self.previous["depth"] + self.maxTopo = self.previous["maxTopo"] + self.minTopo = self.previous["minTopo"] + self.threeD = self.previous["threeD"] + self.data_pretopo = self.previous["data_pretopo"] + self.twtt_pretopo = self.previous["twtt_pretopo"] + # Make sure to not keep deleting history + # when applying undo several times. + histsav = copy.copy(self.previous["history"]) + del histsav[-1] + self.history = histsav + print("undo") + + + def initPrevious(self): + ''' + Initialization of data strucure that contains the step + before the most recent action. + ''' + self.previous["frequencies"] = [] + self.previous["time_zero"] = self.t0_ + self.previous["data"] = self.data + self.previous["twtt"] = self.twtt + self.previous["info"] = self.info + self.previous["profilePos"] = self.profilePos + self.previous["velocity"] = self.velocity + self.previous["depth"] = self.depth + self.previous["maxTopo"] = self.maxTopo + self.previous["minTopo"] = self.minTopo + self.previous["threeD"] = self.threeD + self.previous["data_pretopo"] = self.data_pretopo + self.previous["twtt_pretopo"] = self.twtt_pretopo + histsav = copy.copy(self.history) + self.previous["history"] = histsav + + + + def save(self,filename): + ''' + Saves the processed data together with the processing and visualization + history. Warning: The history stored in this file will contain the full + path to the file. + + INPUT: + filename name for .gpr file + ''' + # Saving the objects: + # Want to force the file name .gpr + file_name, file_ext = os.path.splitext(filename) + if not(file_ext=='.gpr'): + filename = filename + '.gpr' + with open(filename, 'wb') as f: + pickle.dump([self.data, self.info, self.profilePos, self.twtt, + self.history, self.antsep, self.velocity, self.depth, + self.maxTopo, self.minTopo, self.threeD, self.data_pretopo, + self.twtt_pretopo], f) + print("Saved " + filename) + # Add to history string + histstr = "mygpr.save('%s')" %(filename) + self.history.append(histstr) + + + # This is a helper function + def prepProfileFig(self, color="gray", contrast=1.0, yrng=None, xrng=None, asp=None): + ''' + This is a helper function. + It prepares the plot showing the processed profile data. + + INPUT: + color "gray", or "bwr" for blue-white-red, + or any other Matplotlib color map [default: "gray"] + contrast Factor to increase contrast by reducing color range. + [default = 1.0] + yrng y-axis range to show [default: None, meaning "everything"] + xrng x-axis range to show [default: None, meaning "everything"] + asp aspect ratio [default: None, meaning automatic] + + OUTPUT: + contrast contrast value used to prepare the figure + color color value used to prepare the figure + yrng yrng value used to prepare the figure + xrng xrng value used to prepare the figure + asp asp value used to prepare the figure + + ''' + dx=self.profilePos[3]-self.profilePos[2] + dt=self.twtt[3]-self.twtt[2] + stdcont = np.nanmax(np.abs(self.data)[:]) + + if self.velocity is None: + plt.imshow(self.data,cmap=color,extent=[min(self.profilePos)-dx/2.0, + max(self.profilePos)+dx/2.0, + max(self.twtt)+dt/2.0, + min(self.twtt)-dt/2.0], + aspect="auto",vmin=-stdcont/contrast, vmax=stdcont/contrast) + plt.gca().set_ylabel("two-way travel time [ns]") + plt.gca().invert_yaxis() + if yrng is not None: + yrng=[np.max(yrng),np.min(yrng)] + else: + yrng=[np.max(self.twtt),np.min(self.twtt)] + + elif self.maxTopo is None: + dy=dt*self.velocity + plt.imshow(self.data,cmap=color,extent=[min(self.profilePos)-dx/2.0, + max(self.profilePos)+dx/2.0, + max(self.depth)+dy/2.0, + min(self.depth)-dy/2.0], + aspect="auto",vmin=-stdcont/contrast, vmax=stdcont/contrast) + plt.gca().set_ylabel("depth [m]") + plt.gca().invert_yaxis() + if yrng is not None: + yrng=[np.max(yrng),np.min(yrng)] + else: + yrng=[np.max(self.depth),np.min(self.depth)] + + else: + dy=dt*self.velocity + plt.imshow(self.data,cmap=color,extent=[min(self.profilePos)-dx/2.0, + max(self.profilePos)+dx/2.0, + self.minTopo-max(self.depth)-dy/2.0, + self.maxTopo-min(self.depth)+dy/2.0], + aspect="auto",vmin=-stdcont/contrast, vmax=stdcont/contrast) + plt.gca().set_ylabel("elevation [m]") + if yrng is None: + yrng=[self.minTopo-np.max(self.depth),self.maxTopo-np.min(self.depth)] + + + if xrng is None: + xrng=[min(self.profilePos),max(self.profilePos)] + + if yrng is not None: + plt.ylim(yrng) + + if xrng is not None: + plt.xlim(xrng) + + if asp is not None: + plt.gca().set_aspect(asp) + + plt.gca().get_xaxis().set_visible(True) + plt.gca().get_yaxis().set_visible(True) + plt.gca().set_xlabel("profile position [m]") + plt.gca().xaxis.tick_top() + plt.gca().xaxis.set_label_position('top') + + return contrast, color, yrng, xrng, asp + + + def showProfile(self, **kwargs): + ''' + Plots the profile using Matplotlib. + You need to run .show() afterward to show it + + INPUT: + color "gray", or "bwr" for blue-white-red, + or any other Matplotlib color map [default: "gray"] + contrast Factor to increase contrast by reducing color range. + [default = 1.0] + yrng y-axis range to show [default: None, meaning "everything"] + xrng x-axis range to show [default: None, meaning "everything"] + asp aspect ratio [default: None, meaning automatic] + + ''' + self.prepProfileFig(**kwargs) + plt.show(block=False) + + + def printProfile(self, figname, dpi=600, **kwargs): + ''' + Creates a pdf of the profile. + + INPUT: + figname file name for the pdf + dpi dots per inch resolution [default: 600 dpi] + color "gray", or "bwr" for blue-white-red, + or any other Matplotlib color map [default: "gray"] + contrast Factor to increase contrast by reducing color range. + [default = 1.0] + yrng y-axis range to show [default: None, meaning "everything"] + xrng x-axis range to show [default: None, meaning "everything"] + asp aspect ratio [default: None, meaning automatic] + + ''' + contrast, color, yrng, xrng, asp = self.prepProfileFig(**kwargs) + plt.savefig(figname, format='pdf', dpi=dpi) + plt.close('all') + # Put what you did in history + if asp is None: + histstr = "mygpr.printProfile('%s', color='%s', contrast=%g, yrng=[%g,%g], xrng=[%g,%g], dpi=%d)" %(figname,color,contrast,yrng[0],yrng[1],xrng[0],xrng[1],dpi) + else: + histstr = "mygpr.printProfile('%s', color='%s', contrast=%g, yrng=[%g,%g], xrng=[%g,%g], asp=%g, dpi=%d)" %(figname,color,contrast,yrng[0],yrng[1],xrng[0],xrng[1],asp,dpi) + self.history.append(histstr) + + + ####### Processing ####### + def pick_frequencies(self, trace_pour_afficher): + """ + @author: aline + Plot FFT of a chosen trace given in parameter + Click twice in the plot to chose cutting frequencies + To use before filtering with bandpass filter + + Parameters + ---------- + trace_pour_afficher : TYPE int + DESCRIPTION. the number of the trace to plot + + Returns + ------- + self.frequencies_ : TYPE list + DESCRIPTION. [fmin, fmax] for bandpass filter + + """ + #x = self.twtt_before_set_zero + #y = self.data_before_set_zero.transpose()[0] + x = self.twtt + y = self.data.transpose()[0] + + num = len(x) + N = filters.nextpow2(num) + dt = x[1]-x[0] + df = 1/(N*dt) + Nyq = 1/(2*dt) + + fpos = np.linspace(0,Nyq-df,N//2) + #fpos = np.linspace(0,Nyq-df,num//2) + #fpos2 = np.linspace(0,2*Nyq-df,N) + #fpos3 = np.linspace(-Nyq,Nyq-df,N) + + A = fftpack.fft(y,N).tolist()[0] + A_shift = fftpack.fftshift(A) + + #Xpos = A_shift[n:] + Xpos = A_shift[N//2:] + + fig, ax = plt.subplots() + plt.plot(fpos,abs(np.real(Xpos))) + plt.xlabel("Frequency (GHz)") + plt.title('FFT trace %i' %trace_pour_afficher) + #plt.plot(fpos2,A) + #plt.plot(fpos3, A_shift) + + + def onclick(event): + #global frequencies + #pick des fréquences en GHz + print('frequence (GHz): ', event.xdata) + self.frequencies_.append(event.xdata) + + if len(self.frequencies_)==2 : + fig.canvas.mpl_disconnect(cid) + plt.close() + return self.frequencies_ + + cid = fig.canvas.mpl_connect('button_press_event', onclick) + + plt.show(block=True) + + + def pick_time_zero(self, trace_pour_afficher): + """ + @author: aline + Plot a trace (number of trace given in parameter) + Click on the signal to chose time zero + To use before doing 'set time zero' + + Parameters + ---------- + trace_pour_afficher : TYPE int + DESCRIPTION. numéro de trace à afficher + + Returns + ------- + self.t0_ : TYPE float + DESCRIPTION. the number of the trace to plot + + """ + #trace_pour_afficher = int(input('Numéro de la trace à afficher : ')) + x = self.twtt + y = self.data.transpose()[int(trace_pour_afficher-1)] + + fig, ax = plt.subplots() + plt.plot(x, y.tolist()[0]) + plt.xlim((0,80)) + + def onclick(event): + t0 = event.xdata + print('t0 : ', t0 , 'ns') + + fig.canvas.mpl_disconnect(cid) + plt.close() + self.t0_ = t0 + return self.t0_ + + cid = fig.canvas.mpl_connect('button_press_event', onclick) + + plt.show(block=True) + + + def bandpass(self, trace_pour_afficher): + """ + @author: aline + Butterworth bandpass filter + To use after picking cutting frequencies (with pick_frequencies) + Plots the trace chosen before and after filtering and Butterworth filter applied + + Returns + ------- + None. + + """ + # Store previous state for undo + self.storePrevious() + + fmin = self.frequencies_[0] #GHz + fmax = self.frequencies_[1] #GHz + + trace_filtered = [] + for n_trace in range(0, self.data.shape[1]): + datain_filtered, Xpos_filt, N = filters.filtre_passband_calcul(self.twtt,self.data.transpose(), fmin, fmax, self.freq_ech, 2, n_trace, trace_pour_afficher) + trace_filtered.append(datain_filtered[0]) + + #conversion de la liste d'array en matrice + transposée + self.data = np.asmatrix(trace_filtered).transpose() + print("data filtered") + + # Put in history + histstr = "mygpr.filter(%g,%g)" %(fmin,fmax) + self.history.append(histstr) + + + def adjProfile(self,minPos,maxPos): + ''' + Adjusts the length of the profile. + + INPUT: + minPos starting position of the profile + maxpos end position of the profile + ''' + # Store previous state for undo + self.storePrevious() + # set new profile positions + self.profilePos = np.linspace(minPos,maxPos,len(self.profilePos)) + # Put what you did in history + histstr = "mygpr.adjProfile(%g,%g)" %(minPos,maxPos) + self.history.append(histstr) + + + def flipProfile(self): + ''' + Flips the profile left-to-right (start to end) + ''' + # Flips the profile left to right (start to end) + self.storePrevious() + self.data=np.flip(self.data,1) + if self.data_pretopo is not None: + self.data_pretopo = np.flip(self.data_pretopo,1) + histstr = "mygpr.flipProfile()" + self.history.append(histstr) + + + def alignTraces(self): + ''' + Aligns the traces in the profile such that their maximum + amplitudes align at the average two-way travel time of the + maximum amplitudes. + ''' + # Store previous state for undo + self.storePrevious() + self.data = tools.alignTraces(self.data) + # Put what you did in history + histstr = "mygpr.alignTraces()" + self.history.append(histstr) + + + + def cut(self,minPos,maxPos): + ''' + Removes all data outside of the profile positions between + minPos and maxPos. + + INPUT: + minPos starting position of data to keep + maxPos end position of data to keep + ''' + # Store previous state for undo + self.storePrevious() + zeroind = np.abs(self.profilePos - minPos).argmin() + maxind = np.abs(self.profilePos - maxPos).argmin() + self.data = self.data[:,zeroind:(maxind+1)] + self.profilePos=self.profilePos[zeroind:(maxind+1)] + if self.data_pretopo is not None: + self.data_pretopo = self.data_pretopo[:,zeroind:(maxind+1)] + # Put into history string + histstr = "mygpr.cut(%g,%g)" %(minPos,maxPos) + self.history.append(histstr) + + + def setZeroTime(self): + #def setZeroTime(self, newZeroTime): + ''' + Deletes all data recorded before newZeroTime and + sets newZeroTime to zero. + + INPUT: + newZeroTime The new zero-time + ''' + + # Store previous state for undo + self.storePrevious() + + newZeroTime = self.t0_ + + # Find index of value that is nearest to newZeroTime + zeroind = np.abs(self.twtt - newZeroTime).argmin() + # Cut out everything before + self.twtt = self.twtt[zeroind:] - newZeroTime + # Set first value to 0 + self.twtt[0] = 0 + self.data = self.data[zeroind:,:] + + # Put what you did in history + histstr = "mygpr.setZeroTime(%g)" %(newZeroTime) + self.history.append(histstr) + + + def dewow(self,window): + ''' + Subtracts from each sample along each trace an + along-time moving average. + + Can be used as a low-cut filter. + + INPUT: + window length of moving average window + [in "number of samples"] + ''' + # Store previous state for undo + self.storePrevious() + self.data = tools.dewow(self.data,window) + # Put in history + histstr = "mygpr.dewow(%d)" %(window) + self.history.append(histstr) + + + def smooth(self,window): + ''' + Replaces each sample along each trace with an + along-time moving average. + + Can be used as high-cut filter. + + INPUT: + window length of moving average window + [in "number of samples"] + ''' + # Store previous state for undo + self.storePrevious() + self.data = tools.smooth(self.data,window) + # Put in history + histstr = "mygpr.smooth(%d)" %(window) + self.history.append(histstr) + + + def remMeanTrace(self,ntraces): + ''' + Subtracts from each trace the average trace over + a moving average window. + + Can be used to remove horizontal arrivals, + such as the airwave. + + INPUT: + ntraces window width; over how many traces + to take the moving average. + ''' + # Store previous state for undo + self.storePrevious() + # apply + self.data = tools.remMeanTrace(self.data,ntraces) + # Put in history + histstr = "mygpr.remMeanTrace(%d)" %(ntraces) + self.history.append(histstr) + + + def profileSmooth(self,ntraces,noversample): + ''' + First creates copies of each trace and appends the copies + next to each trace, then replaces each trace with the + average trace over a moving average window. + + Can be used to smooth-out noisy reflectors appearing + in neighboring traces, or simply to increase the along-profile + resolution by interpolating between the traces. + + For example: To increase the along-profile resolution smoothly + by a factor of 4: use + + mygpr.profileSmooth(4,4) + + INPUT: + ntraces window width [in "number of samples"]; + over how many traces to take the moving average. + noversample how many copies of each trace + ''' + # Store previous state for undo + self.storePrevious() + self.data,self.profilePos = tools.profileSmooth(self.data,self.profilePos, + ntraces,noversample) + # Put in history + histstr = "mygpr.profileSmooth(%d,%d)" %(ntraces,noversample) + self.history.append(histstr) + + + def tpowGain(self,power=0.0): + ''' + Apply a t-power gain to each trace with the given exponent. + + INPUT: + power exponent + ''' + # Store previous state for undo + self.storePrevious() + # apply tpowGain + self.data = tools.tpowGain(self.data,self.twtt,power) + # Put in history + histstr = "mygpr.tpowGain(%g)" %(power) + self.history.append(histstr) + + def agcGain(self,window=10): + ''' + Apply automated gain controll (AGC) by normalizing the energy + of the signal over a given window width in each trace + + INPUT: + window window width [in "number of samples"] + ''' + # Store previous state for undo + self.storePrevious() + # apply agcGain + self.data = tools.agcGain(self.data,window) + # Put in history + histstr = "mygpr.agcGain(%d)" %(float(window)) + self.history.append(histstr) + + + def setVelocity(self,velocity): + ''' + Provide the subsurface RMS radar velocity + + INPUT: + velocity subsurface RMS velocity [in m/ns] + ''' + # Store previous state for undo + self.storePrevious() + + self.velocity = velocity + self.depth = self.twtt * velocity/2.0 + + # Put in history + histstr = "mygpr.setVelocity(%g)" %(velocity) + self.history.append(histstr) + + + def antennaSep(self): + ''' + Corrects for distortions of arrival times caused by the + separation of the antennae. + + For this to work properly, you must have set the velocity + and you must have set the zero time to the beginning of the + arrival of the airwave. + + ''' + + # Store previous state for undo + self.storePrevious() + + # Take into account that the airwave first break + # is after the airwave has already traveled the + # antenna separation with the speed of light 0.3 m/ns. + # And we only look at half the + # two-way travel time. Hence divide by two + t0 = self.twtt/2 + self.antsep/(2*0.3) + + # t0 is when the waves left the transmitter antenna. + # To be able to calculate the depth time from the + # single-way travel time we need to shift the time reference + # frame. Lets set it "arriving at midpoint time", so + ta = t0 + self.antsep/(2*self.velocity) + # Later we will need to undo this reference frame transformation + + # Now use the pythagorean relationship between single-way travel + # time to the depth point and the depth time td + tad = np.sqrt( ta**2 - (self.antsep/(2*self.velocity))**2 ) + + # We are still in the "arriving at midpoint" time frame ta + # To transform ta into depth time td, we need to shift it back + # by the time it took for the ground wave to get to the midpoint. + # This makes sense because the times before the groundwave got to the + # midpoint will not actually be underground in the sense of: + # No travel into depth has been recorded at the receiver. + # These "arrivals" will just be shifted into "negative arrival times" + # and hence "negative depth" + td = tad - self.antsep/(2*self.velocity) + + # Finally, translate time into depth + self.depth = td*self.velocity + + # And update the two-way travel time + self.twtt = td + + # Put in history + histstr = "mygpr.antennaSep()" + self.history.append(histstr) + + + def fkMigration(self): + ''' + Apply Stolt's f-k migration to the profile. Requires the + velocity to be set. + + This is a wrapper function for the migration code + imported from Nat Wilson's irlib software. + ''' + # Store previous state for undo + self.storePrevious() + # apply migration + dt=self.twtt[3]-self.twtt[2] + #dx=self.profilePos[1]-self.profilePos[0] + dx=(self.profilePos[-1]-self.profilePos[0])/(len(self.profilePos)-1) + # fkmig sets x profile to start at zero but resamples + self.data,self.twtt,migProfilePos=mig_fk.fkmig(self.data,dt,dx,self.velocity) + self.profilePos = migProfilePos + self.profilePos[0] + + # Put in history + histstr = "mygpr.fkMigration()" + self.history.append(histstr) + + + def truncateY(self,maxY): + ''' + Delete all data after y-axis position maxY. + + INPUT: + maxY maximum y-axis position for data to be kept + ''' + # Store previous state for undo + self.storePrevious() + if self.velocity is None: + maxtwtt = maxY + maxind = np.argmin( np.abs(self.twtt-maxY) ) + self.twtt = self.twtt[0:maxind] + # Set the last value to maxY + self.twtt[-1] = maxY + self.data = self.data[0:maxind,:] + else: + maxtwtt = maxY*2.0/self.velocity + maxind = np.argmin( np.abs(self.twtt-maxtwtt) ) + self.twtt = self.twtt[0:maxind] + # Set the last value to maxtwtt + self.twtt[-1] = maxtwtt + self.data = self.data[0:maxind,:] + self.depth = self.depth[0:maxind] + self.depth[-1] = maxY + # Put in history + histstr = "mygpr.truncateY(%g)" %(maxY) + self.history.append(histstr) + + + + def topoCorrect(self,topofile,delimiter=','): + ''' + Correct for topography along the profile by shifting each + Trace up or down depending on a provided ASCII text file + containing topography information. + + The topography data file can either have + two columns: profile position and topography, + or three columns: X, Y, and topography, or Easting, Northing, topography + + In the topo text file, units of profile position (or northing and easting) + and of the topography (or elevation) need to be in meters! + + Requires the velocity to be set. + + INPUT: + topofile file name for ASCII text topography information + delimiter how the entries are delimited (by comma, or by tab) + [default: ',']. To set tab: delimiter='\t' + ''' + if self.velocity is None: + print("First need to set velocity!") + return + # Store previous state for undo + self.storePrevious() + self.data_pretopo = self.data + self.twtt_pretopo = self.twtt + topoPos, topoVal, self.threeD = tools.prepTopo(topofile,delimiter,self.profilePos[0]) + self.data, self.twtt, self.maxTopo, self.minTopo = tools.correctTopo(self.data, + velocity=self.velocity, + profilePos=self.profilePos, + topoPos=topoPos, + topoVal=topoVal, + twtt=self.twtt) + # Put in history + if delimiter is ',': + histstr = "mygpr.topoCorrect('%s')" %(topofile) + else: + histstr = "mygpr.topoCorrect('%s',delimiter='\\t')" %(topofile) + self.history.append(histstr) + + + + def exportVTK(self,outfile,gpsinfo,delimiter=',',thickness=0,aspect=1.0,smooth=True, win_length=51, porder=3): + ''' + Turn processed profile into a VTK file that can be imported in + Paraview or MayaVi or other VTK processing and visualization tools. + + If three-dimensional topo information is provided (X,Y,Z or + Easting, Northing, Elevation), then the profile will be exported + in its three-dimensional shape. + + INPUT: + outfile file name for the VTK file + gpsinfo EITHER: n x 3 matrix containing x, y, and z or + Easting, Northing, Elevation information + OR: file name for ASCII text file containing this + information + delimiter if topo file is provided: delimiter (by comma, or by tab) + [default: ',']. To set tab: delimiter='\t' + thickness If you want your profile to be exported as a + three-dimensional band with thickness, enter thickness + in meters [default: 0] + aspect aspect ratio in case you want to exaggerate z-axis. + default = 1. I recommend leaving this at 1 and using + your VTK visualization software to set the aspect for + the representation. + smooth Want to smooth the profile's three-dimensional alignment + instead of piecewise linear? [Default: True] + win_length If smoothing, the window length for + scipy.signal.savgol_filter [default: 51] + porder If smoothing, the polynomial order for + scipy.signal.savgol_filter [default: 3] + ''' + # If gpsmat is a filename, we first need to load the file: + if type(gpsinfo) is str: + gpsmat = np.loadtxt(gpsinfo,delimiter=delimiter) + else: + gpsmat = gpsinfo + + # First get the x,y,z positions of our data points + x,y,z = tools.prepVTK(self.profilePos,gpsmat,smooth,win_length,porder) + z = z*aspect + if self.velocity is None: + downward = self.twtt*aspect + else: + downward = self.depth*aspect + Z = np.reshape(z,(len(z),1)) - np.reshape(downward,(1,len(downward))) + + + if thickness: + ZZ = np.tile(np.reshape(Z, (1,Z.shape[0],Z.shape[1])), (2,1,1)) + else: + ZZ = np.tile(np.reshape(Z, (1,Z.shape[0],Z.shape[1])), (1,1,1)) + + # This is if we want everything on the x axis. + #X = np.tile(np.reshape(self.profilePos,(len(self.profilePos),1)),(1,len(downward))) + #XX = np.tile(np.reshape(X, (X.shape[0],1,X.shape[1])), (1,2,1)) + #YY = np.tile(np.reshape([-thickness/2,thickness/2],(1,2,1)), (len(x),1,len(downward))) + + # To create a 3D grid with a width, calculate the perpendicular direction, + # normalize it, and add it to xvals and yvals as below. + # To figure this out, just drar the profile point-by-point, and at each point, + # draw the perpendicular to the segment and place a grid point in each perpendicular + # direction + # + # x[0]-px[0], x[1]-px[1], x[2]-px[2], ..... + # xvals = x[0] , x[1] , x[2] , ..... + # x[0]+px[0], x[1]+px[1], x[2]+px[2], ..... + # + # y[0]+py[0], y[1]+py[1], y[2]+py[2], ..... + # yvals = y[0] , y[1] , y[2] , ..... + # y[0]-py[0], y[1]-py[1], y[2]-py[2], ..... + # + # Here, the [px[i],py[i]] vector needs to be normalized by the thickness + if thickness: + pvec = np.asarray([(y[0:-1]-y[1:]).squeeze(), (x[1:]-x[0:-1]).squeeze()]) + pvec = np.divide(pvec, np.linalg.norm(pvec,axis=0)) * thickness/2.0 + # We can't calculate the perpendicular direction at the last point + # let's just set it to the same as for the second-to-last point + pvec = np.append(pvec, np.expand_dims(pvec[:,-1],axis=1) ,axis=1) + X = np.asarray([(x.squeeze()-pvec[0,:]).squeeze(), (x.squeeze()+pvec[0,:]).squeeze()]) + Y = np.asarray([(y.squeeze()+pvec[1,:]).squeeze(), (y.squeeze()-pvec[1,:]).squeeze()]) + else: + X = np.asarray([x.squeeze()]) + Y = np.asarray([y.squeeze()]) + + # Copy-paste the same X and Y positions for each depth + XX = np.tile(np.reshape(X, (X.shape[0],X.shape[1],1)), (1,1,ZZ.shape[2])) + YY = np.tile(np.reshape(Y, (Y.shape[0],Y.shape[1],1)), (1,1,ZZ.shape[2])) + + if self.maxTopo is None: + data=self.data.transpose() + else: + data=self.data_pretopo.transpose() + + data = np.asarray(data) + data = np.reshape(data,(1,data.shape[0],data.shape[1])) + data = np.tile(data, (2,1,1)) + + # Remove the last row and column to turn it into a cell + # instead of point values + data = data[0:-1,0:-1,0:-1] + + nx=2-1 + ny=len(x)-1 + nz=len(downward)-1 + datarray = np.zeros(nx*ny*nz).reshape(nx,ny,nz) + datarray[:,:,:] = data + + gridToVTK(outfile,XX,YY,ZZ, cellData ={'gpr': datarray}) + + # Put in history + if gpsinfo is None: + histstr = "mygpr.exportVTK('%s',aspect=%g)" %(outfile,aspect) + else: + if type(gpsinfo) is str: + if delimiter is ',': + histstr = "mygpr.exportVTK('%s',gpsinfo='%s',thickness=%g,delimiter=',',aspect=%g,smooth=%r, win_length=%d, porder=%d)" %(outfile,gpsinfo,thickness,aspect,smooth,win_length,porder) + else: + histstr = "mygpr.exportVTK('%s',gpsinfo='%s',thickness=%g,delimiter='\\t',aspect=%g,smooth=%r, win_length=%d, porder=%d)" %(outfile,gpsinfo,thickness,aspect,smooth,win_length,porder) + else: + if delimiter is ',': + histstr = "mygpr.exportVTK('%s',gpsinfo=mygpr.threeD,thickness=%g,delimiter=',',aspect=%g,smooth=%r, win_length=%d, porder=%d)" %(outfile,thickness,aspect,smooth,win_length,porder) + else: + histstr = "mygpr.exportVTK('%s',gpsinfo=mygpr.threeD,thickness=%g,delimiter='\\t',aspect=%g,smooth=%r, win_length=%d, porder=%d)" %(outfile,thickness,aspect,smooth,win_length,porder) + + self.history.append(histstr) + + def storePrevious(self): + ''' + Stores the current state of the profile and history in the + "previous" variable to be able to apply undo. + ''' + self.previous["data"] = self.data + self.previous["twtt"] = self.twtt + self.previous["info"] = self.info + self.previous["profilePos"] = self.profilePos + self.previous["history"] = self.history + self.previous["velocity"] = self.velocity + self.previous["depth"] = self.depth + self.previous["maxTopo"] = self.maxTopo + self.previous["threeD"] = self.threeD + self.previous["data_pretopo"] = self.data_pretopo + self.previous["twtt_pretopo"] = self.twtt_pretopo + + + + + +class gprpyCW(gprpyProfile): + ''' + Ground penetrating radar data processing and visualization class + for common midpoint or wide angle reflection and refraction data + + Inherits all of the gprpyProfile class functions but not all + of these functions may be useful here. + ''' + def __init__(self,filename=None,dtype=None): + ''' + Initialization for a gprpyCW object. Initialization can be + empty or with a provided filename for the GPR data and + a data type. + + INPUT: + filename data file name. Currently supported formats: + .gpr (GPRPy), .DT1 (SnS), .DZT (GSSI), .rd3 (MALA), + and ENVI standard BSQ. + dtype data type. Either "CMP" or "WARR" + ''' + # Inheriting the initializer from the gprpyProfile class + super().__init__(None) + self.history = ["mygpr = gp.gprpyCW()"] + # Initialize previous for undo + self.previous = {} + self.dtype = dtype + # Stacked amplitude plots + self.linStAmp = None + self.hypStAmp = None + # Picked lines and hyperbolae + self.lins = list() + self.hyps = list() + + if (filename is not None) and (dtype is not None): + self.importdata(filename,dtype) + + + def storePrevious(self): + ''' + Stores the current state of the profile and history in the + "previous" variable to be able to apply undo. + ''' + self.previous["data"] = self.data + self.previous["twtt"] = self.twtt + self.previous["info"] = self.info + self.previous["profilePos"] = self.profilePos + self.previous["history"] = self.history + self.previous["dtype"] = self.dtype + self.previous["hypStAmp"] = self.hypStAmp + self.previous["linStAmp"] = self.linStAmp + self.previous["lins"] = self.lins + self.previous["hyps"] = self.hyps + + + + + def importdata(self,filename,dtype): + ''' + Loads .gpr (native GPRPy), .DT1 (Sensors and Software), + .DZT (GSSI), .GPRhdr (ENVI standard BSQ), .rad (MALA) + data files and populates all the gprpyProfile fields. + + INPUT: + filename name of the .gpr, .DT1, .DZT, .GPRhdr, dat, rd3, + or .rad file you want to import. + The header file name and the data file name + have to be the same! + dtype data type. Either "CMP" or "WARR + + ''' + print(filename) + super().importdata(filename) + self.dtype = dtype + self.vVals = None + # Remove the history string from the super-class importing + del self.history[-1] + # Put what you did in history + histstr = "mygpr.importdata('%s',dtype='%s')" %(filename,dtype) + self.history.append(histstr) + + + + + #def setZeroTimeCW(self,newZeroTime): + def setZeroTimeCW(self,newZeroTime): + ''' + Deletes all data recorded before newZeroTime and + sets newZeroTime to zero. + + INPUT: + newZeroTime The new zero-time + ''' + # Store previous state for undo + self.storePrevious() + # Find index of value that is nearest to newZeroTime + + zeroind = np.abs(self.twtt - newZeroTime).argmin() + # Cut out everything before + self.twtt = self.twtt[zeroind:] - newZeroTime + #self.data = self.data[zeroind:,:] + # Put what you did in history + histstr = "mygpr.setZeroTime(%g)" %(newZeroTime) + self.history.append(histstr) + + + def normalize(self): + ''' + Divides each trace by its total energy to counteract the + loss of energy for wider antennae separations. + ''' + # Store previous state for undo + self.storePrevious() + # Calculate norm of each trace and divide each trace by it + self.data = np.divide(self.data,np.maximum(np.linalg.norm(self.data,axis=0),1e-8)) + print("normalized data set") + # Put what you did in history + histstr = "mygpr.normalize()" + self.history.append(histstr) + + + def linStackedAmplitude(self,vmin=0.01,vmax=0.35,vint=0.01): + ''' + Calculates the linear stacked amplitudes for each two-way + travel time sample and the provided velocity range + by summing the pixels of the data that follow a line given + by the two-way travel time zero offset and the velocity. + + INPUT: + vmin minimal velocity for which to calculate the + stacked amplitude, in m/ns [default = 0.01 m/ns] + vmax maximum velocity for which to calculate the + stacked amplitude, in m/ns [default = 0.35 m/ns] + vint velocity intervall, in m/ns [default = 0.01 m/ns] + ''' + # Store previous state for undo + self.storePrevious() + self.vVals = np.arange(vmin,vmax+vint,vint) + if self.dtype is "WARR": + typefact = 1 + elif self.dtype is "CMP": + typefact = 2 + self.linStAmp = tools.linStackedAmplitude(self.data,self.profilePos,self.twtt,self.vVals,self.twtt,typefact) + print("calculated linear stacked amplitude") + # Put what you did in history + histstr = "mygpr.linStackedAmplitude(vmin=%g,vmax=%g,vint=%g)" %(vmin,vmax,vint) + self.history.append(histstr) + + + def hypStackedAmplitude(self,vmin=0.01,vmax=0.35,vint=0.01): + ''' + Calculates the hyperbolic stacked amplitudes for each two-way + travel time sample and the provided velocity range + by summing the pixels of the data that follow a hyperbola given + by the two-way travel time apex and the velocity. + + INPUT: + vmin minimal velocity for which to calculate the + stacked amplitude, in m/ns [default = 0.01 m/ns] + vmax maximum velocity for which to calculate the + stacked amplitude, in m/ns [default = 0.35 m/ns] + vint velocity intervall, in m/ns [default = 0.01 m/ns] + ''' + # Store previous state for undo + self.storePrevious() + self.vVals = np.arange(vmin,vmax+vint,vint) + if self.dtype is "WARR": + typefact = 1 + elif self.dtype is "CMP": + typefact = 2 + self.hypStAmp = tools.hypStackedAmplitude(self.data,self.profilePos,self.twtt,self.vVals,self.twtt,typefact) + print("calculated hyperbola stacked amplitude") + # Put what you did in history + histstr = "mygpr.hypStackedAmplitude(vmin=%g,vmax=%g,vint=%g)" %(vmin,vmax,vint) + self.history.append(histstr) + + + def addLin(self,zerotwtt,vel): + ''' + Adds an observed line given by its zero-offset two-way travel + time and velocity to the list of lines. + + INPUT: + zerotwtt the zero-offset two-way travel time (intercept) + of the observed line + vel the velocity (inverse slope) of the observed line + ''' + # Store previous state for undo + self.storePrevious() + self.lins.append([zerotwtt,vel]) + # Put what you did in history + histstr = "mygpr.addLin(zerotwtt=%g,vel=%g)" %(zerotwtt,vel) + self.history.append(histstr) + + + def addHyp(self,zerotwtt,vel): + ''' + Adds an observed hyperbola given by its apex two-way + travel time and velocity to the list of lines. + + INPUT: + zerotwtt the apex two-way travel time of the observed line + vel the velocity of the observed line + ''' + # Store previous state for undo + self.storePrevious() + self.hyps.append([zerotwtt,vel]) + # Put what you did in history + histstr = "mygpr.addHyp(zerotwtt=%g,vel=%g)" %(zerotwtt,vel) + self.history.append(histstr) + + def remLin(self): + ''' + Removes the most recently added observed line from the list + of observed lines + ''' + # Store previous state for undo + self.storePrevious() + del self.lins[-1] + # Put what you did in history + histstr = "mygpr.remLin()" + self.history.append(histstr) + + def remHyp(self): + ''' + Removes the most recently added observed hyperbola from the list + of observed lines + ''' + # Store previous state for undo + self.storePrevious() + del self.hyps[-1] + # Put what you did in history + histstr = "mygpr.remHyp()" + self.history.append(histstr) + + + + + # This is a helper function + def prepCWFig(self, contrast=1.0, color="gray", yrng=None, xrng=None, showlnhp=False): + ''' + This is a helper function. + It prepares the plot showing the processed CMP or WARR data. + + INPUT: + color "gray", or "bwr" for blue-white-red, + or any other Matplotlib color map [default: "gray"] + contrast Factor to increase contrast by reducing color range. + [default = 1.0] + yrng y-axis range to show [default: None, meaning "everything"] + xrng x-axis range to show [default: None, meaning "everything"] + showlnhp show the observed lines and hyperbolae from the list + [default: False] + + OUTPUT: + contrast contrast value used to prepare the figure + color color value used to prepare the figure + yrng yrng value used to prepare the figure + xrng xrng value used to prepare the figure + showlnhp showlnhp value used to prepare the figure + ''' + dx=self.profilePos[3]-self.profilePos[2] + dt=self.twtt[3]-self.twtt[2] + stdcont = np.nanmax(np.abs(self.data)[:]) + + plt.imshow(self.data,cmap=color,extent=[min(self.profilePos)-dx/2.0, + max(self.profilePos)+dx/2.0, + max(self.twtt)+dt/2.0, + min(self.twtt)-dt/2.0], + aspect="auto",vmin=-stdcont/contrast, vmax=stdcont/contrast) + plt.gca().set_ylabel("two-way travel time [ns]") + plt.gca().invert_yaxis() + if yrng is not None: + yrng=[np.max(yrng),np.min(yrng)] + else: + yrng=[np.max(self.twtt),np.min(self.twtt)] + plt.ylim(yrng) + + if xrng is None: + xrng=[min(self.profilePos),max(self.profilePos)] + plt.xlim(xrng) + + plt.gca().get_xaxis().set_visible(True) + plt.gca().get_yaxis().set_visible(True) + if self.dtype == "WARR": + plt.gca().set_xlabel("antenna separation [m]") + typefact=1 + elif self.dtype == "CMP": + plt.gca().set_xlabel("distance from midpoint [m]") + typefact=2 + plt.gca().xaxis.tick_top() + plt.gca().xaxis.set_label_position('top') + + # Show hyperbolae / lines if you want + if showlnhp: + if self.lins: + for lin in self.lins: + time = lin[0] + typefact*self.profilePos/lin[1] + plt.plot(self.profilePos,time,linewidth=2,color='yellow') + plt.plot(self.profilePos,time,linewidth=1,color='black') + if self.hyps: + x2 = np.power(typefact*self.profilePos,2.0) + for hyp in self.hyps: + time = np.sqrt(x2 + 4*np.power(hyp[0]/2.0 * hyp[1],2.0))/hyp[1] + plt.plot(self.profilePos,time,linewidth=2,color='yellow') + plt.plot(self.profilePos,time,linewidth=1,color='black') + + return contrast, color, yrng, xrng, showlnhp + + + + def prepStAmpFig(self, whichstamp="lin", saturation=1.0, yrng=None, vrng=None): + ''' + This is a helper function. + It prepares the plot showing the stacked amplitudes results. + + INPUT: + whichstamp is this for the linear ("lin") or hyperbolic ("hyp") + stacked amplitudes + saturation Factor to increase contrast by reducing color range. + [default = 1.0] + yrng y-axis range to show [default: None, meaning "everything"] + vrng velocities (x-axis) range to show + [default: None, meaning "everything"] + + OUTPUT: + whichstamp whichstamp value used to prepare the figure + saturation saturation value used to prepare the figure + yrng yrng value used to prepare the figure + vrng vrng value used to prepare the figure + ''' + dt=self.twtt[3]-self.twtt[2] + dv=self.vVals[3]-self.vVals[2] + if whichstamp == "lin": + stamp = self.linStAmp + title = "linear stacked amplitude" + elif whichstamp == "hyp": + stamp = self.hypStAmp + title = "hyperbolic stacked amplitude" + else: + stamp = None + + if stamp is not None: + stdcont = np.nanmax(np.abs(stamp)[:]) + plt.imshow(np.flipud(np.abs(stamp)), cmap='inferno', + extent=[np.min(self.vVals)-dv/2.0, + np.max(self.vVals)+dv/2.0, + np.min(self.twtt)-dt/2.0, + np.max(self.twtt)+dt/2.0], + aspect='auto', + vmin=0, vmax=stdcont/saturation) + + if yrng is not None: + yrng=[np.max(yrng),np.min(yrng)] + else: + yrng=[np.max(self.twtt),np.min(self.twtt)] + plt.ylim(yrng) + + if vrng is None: + vrng=[np.min(self.vVals),np.max(self.vVals)] + plt.xlim(vrng) + + + plt.gca().set_xlabel("velocity [m/ns]") + plt.gca().set_ylabel("two-way travel time [ns]") + #plt.gca().invert_yaxis() + plt.gca().set_title(title) + plt.gca().get_xaxis().set_visible(True) + plt.gca().get_yaxis().set_visible(True) + plt.gca().get_xaxis().set_ticks_position('both') + plt.gca().get_yaxis().set_ticks_position('both') + + return whichstamp, saturation, yrng, vrng + + + + def showCWFig(self, **kwargs): + ''' + Plots the CMP or WARR data using Matplotlib. + You need to run .show() afterward to show it + + INPUT: + color "gray", or "bwr" for blue-white-red, + or any other Matplotlib color map [default: "gray"] + contrast Factor to increase contrast by reducing color range. + [default = 1.0] + yrng y-axis range to show [default: None, meaning "everything"] + xrng x-axis range to show [default: None, meaning "everything"] + showlnhp show the observed lines and hyperbolae from the list + [default: False] + ''' + self.prepCWFig(**kwargs) + plt.show(block=False) + + + def showStAmpFig(self, **kwargs): + ''' + Plots the stacked amplitude results using Matplotlib. + You need to run .show() afterward to show it + + INPUT: + whichstamp is this for the linear ("lin") or hyperbolic ("hyp") + stacked amplitudes + saturation Factor to increase contrast by reducing color range. + [default = 1.0] + yrng y-axis range to show [default: None, meaning "everything"] + vrng velocities (x-axis) range to show + [default: None, meaning "everything"] + ''' + self.prepStAmpFig(**kwargs) + plt.show(block=False) + + + + def printCWFigure(self, figname, dpi=600, **kwargs): + ''' + Creates a pdf of the common midpoint or wide angle reflection + and refraction data. + + INPUT: + figname file name for the pdf + dpi dots per inch resolution [default: 600 dpi] + color "gray", or "bwr" for blue-white-red, + or any other Matplotlib color map [default: "gray"] + contrast Factor to increase contrast by reducing color range. + [default = 1.0] + yrng y-axis range to show [default: None, meaning "everything"] + xrng x-axis range to show [default: None, meaning "everything"] + showlnhp show the observed lines and hyperbolae from the list + [default: False] + ''' + + + contrast, color, yrng, xrng, showlnhp = self.prepCWFig(**kwargs) + plt.savefig(figname, format='pdf', dpi=dpi) + plt.close('all') + # Put what you did in history + histstr = "mygpr.printCWFigure('%s', color='%s', contrast=%g, yrng=[%g,%g], xrng=[%g,%g], showlnhp=%r, dpi=%d)" %(figname,color,contrast,yrng[0],yrng[1],xrng[0],xrng[1],showlnhp,dpi) + self.history.append(histstr) + + + def printStAmpFigure(self, figname, dpi=600, **kwargs): + ''' + Creates a pdf of the the stacked amplitude results. + + INPUT: + figname file name for the pdf + dpi dots per inch resolution [default: 600 dpi] + whichstamp is this for the linear ("lin") or hyperbolic ("hyp") + stacked amplitudes + saturation Factor to increase contrast by reducing color range. + [default = 1.0] + yrng y-axis range to show [default: None, meaning "everything"] + vrng velocities (x-axis) range to show + [default: None, meaning "everything"] + ''' + + whichstamp, saturation, yrng, vrng = self.prepStAmpFig(**kwargs) + plt.savefig(figname, format='pdf', dpi=dpi) + plt.close('all') + histstr = "mygpr.printStAmpFigure('%s', whichstamp='%s', saturation=%g, yrng=[%g,%g], vrng=[%g,%g])" %(figname, whichstamp, saturation, yrng[0], yrng[1], vrng[0], vrng[1]) + self.history.append(histstr) + + + + + + diff --git a/gprpyGUI.py b/gprpyGUI.py new file mode 100644 index 0000000..cc48266 --- /dev/null +++ b/gprpyGUI.py @@ -0,0 +1,970 @@ +# import sys +# if sys.version_info[0] < 3: +# import Tkinter as tk +# from Tkinter import filedialog as fd +# else: +# import tkinter as tk +# from tkinter import filedialog as fd + +import sys +import tkinter as tk +from tkinter import filedialog as fd +from tkinter import simpledialog as sd +from tkinter import messagebox as mesbox +import matplotlib as mpl +mpl.use('TkAgg') +from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg +from matplotlib.figure import Figure +import gprpy.gprpy as gp +import numpy as np +import gprpy.toolbox.splash as splash +import os +import Pmw +import scipy.interpolate as interp + + + + +colsp=2 +rightcol=9 +halfwid=6 +figrowsp=20+1 +figcolsp=9 + + +class GPRPyApp: + ''' + GPRPy class for graphical user interface for GPR profile data + ''' + + def __init__(self,master): + self.window = master + + # Set up for high-resolution screens + normscrwidt=1024 + normscrhigt=768 + scrwidt=master.winfo_screenwidth() + scrhigt=master.winfo_screenheight() + # These to use if operating system doesn't automatically adjust + #self.widfac=scrwidt/normscrwidt + #self.highfac=scrhigt/normscrhigt + self.widfac=normscrwidt/normscrhigt + self.highfac=1 + fontfac=(normscrwidt/normscrhigt)/(scrwidt/scrhigt) + + master.title("GPRPy") + + # Variables specific to GUI + self.balloon = Pmw.Balloon() + self.picking = False + self.delimiter = None + self.grid = False + + # Initialize the gprpy + proj = gp.gprpyProfile() + + # Show splash screen + fig=Figure(figsize=(8*self.widfac,5*self.highfac)) + a=fig.add_subplot(111) + dir_path = os.path.dirname(os.path.realpath(__file__)) + splash.showSplash(a,dir_path,self.widfac,self.highfac,fontfac) + + # Set font size for screen res + mpl.rcParams.update({'font.size': mpl.rcParams['font.size']*self.widfac}) + a.tick_params(direction='out',length=6*self.widfac,width=self.highfac) + + a.get_xaxis().set_visible(False) + a.get_yaxis().set_visible(False) + canvas = FigureCanvasTkAgg(fig, master=self.window) + canvas.get_tk_widget().grid(row=2,column=0,columnspan=figcolsp,rowspan=figrowsp,sticky='nsew') + + canvas.draw() + + + + ## Visualization Buttons + + # Undo Button + undoButton = tk.Button( + text="undo", + command=lambda : [self.resetYrng(proj), + self.undo(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + undoButton.config(height = 1, width = 2*halfwid) + undoButton.grid(row=0, column=0, sticky='nsew',rowspan=2) + self.balloon.bind(undoButton, + '"Undoes" the most recent processing step and\n' + 'sets the data back to its previous state.\n' + 'This also removes the most recent processing\n' + 'step from the history. Does not revert\n' + 'visualization settings such as "set x-range"\n' + 'etc.') + + # Full view + FullButton = tk.Button( + text="full view", fg="black", + command=lambda : [self.setFullView(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + FullButton.config(height = 1, width = 2*halfwid) + FullButton.grid(row=0, column=1, sticky='nsew',rowspan=2) + self.balloon.bind(FullButton,"Resets x- and y-axis limits to full data.") + + + # Grid button + GridButton = tk.Button( + text="grid", fg="black", + command=lambda : [self.toggleGrid(), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + GridButton.config(height = 1, width = 2*halfwid) + GridButton.grid(row=0, column=2, sticky='nsew',rowspan=2) + self.balloon.bind(GridButton,"Toggles grid on/off.") + + + # X range + XrngButton = tk.Button( + text="set x-range", fg="black", + command=lambda : [self.setXrng(), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + XrngButton.config(height = 1, width = 2*halfwid) + XrngButton.grid(row=0, column=3, sticky='nsew',rowspan=2) + self.balloon.bind(XrngButton,"Set the x-axis display limits.") + + + # Y range + YrngButton = tk.Button( + text="set y-range", fg="black", + command=lambda : [self.setYrng(), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + YrngButton.config(height = 1, width = 2*halfwid) + YrngButton.grid(row=0, column=4, sticky='nsew',rowspan=2) + self.balloon.bind(YrngButton,"Set the y-axis display limits.") + + + # Aspect + AspButton = tk.Button( + text="aspect ratio", fg="black", + command=lambda : [self.setAspect(), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + AspButton.config(height = 1, width = 2*halfwid) + AspButton.grid(row=0, column=5, sticky='nsew',rowspan=2) + self.balloon.bind(AspButton, "Set the aspect ratio between x- and y-axis.") + + + # Contrast + contrtext = tk.StringVar() + contrtext.set("contrast") + contrlabel = tk.Label(master, textvariable=contrtext,height = 1,width = 2*halfwid) + contrlabel.grid(row=0, column=6, sticky='nsew') + self.balloon.bind(contrlabel,"Set color saturation") + self.contrast = tk.DoubleVar() + contrbox = tk.Entry(master, textvariable=self.contrast, width=2*halfwid) + contrbox.grid(row=1, column=6, sticky='nsew') + #contr.set("1.0") + self.contrast.set("1.0") + + + # Mode switch for figure color + self.color=tk.StringVar() + self.color.set("gray") + colswitch = tk.OptionMenu(master,self.color,"gray","bwr") + colswitch.grid(row=0, column=7, sticky='nsew',rowspan=2) + self.balloon.bind(colswitch, + "Choose between gray-scale\n" + "and red-white-blue (rwb)\n" + "data representation.") + + + # Refreshing plot + plotButton = tk.Button( + text="refresh plot", + command=lambda : self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)) + plotButton.config(height = 1, width = 2*halfwid) + plotButton.grid(row=0, column=8, sticky='nsew',rowspan=2) + self.balloon.bind(plotButton, + "Refreshes the figure after changes\n" + "in the visualization settings. Also\n" + "removes any plotted hyperbolae.") + + + + + ## Methods buttons + + # Load data + LoadButton = tk.Button( + text="import data", fg="black", + command=lambda : [self.loadData(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + LoadButton.config(height = 1, width = 2*halfwid) + LoadButton.grid(row=0, column=rightcol, sticky='nsew',columnspan=colsp,rowspan=2) + self.balloon.bind(LoadButton,"Load .gpr, .DT1, or .DZT data.") + + + # Adjust profile length; if trigger wheel is not good + AdjProfileButton = tk.Button( + text="adj profile", fg="black", + command=lambda : [self.adjProfile(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + AdjProfileButton.config(height = 1, width = 2*halfwid) + AdjProfileButton.grid(row=2, column=rightcol, sticky='nsew',columnspan=colsp) + self.balloon.bind(AdjProfileButton, + "Adjust the profile length to \n" + "known start and end positions\n" + "and/or flip the profile horizontally\n" + "(left to right)") + + # pick zero time + PickZeroTimeButton = tk.Button( + text="Pick t0", fg="black", + command=lambda : [self.pickZeroTime(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + PickZeroTimeButton.config(height = 1, width = halfwid) + PickZeroTimeButton.grid(row=3, column=rightcol, sticky='nsew') + self.balloon.bind(PickZeroTimeButton, + "Pick time zero for a chosen trace \n" + "To do before setting zero time \n" + "Created by AE") + + + # Set new zero time + SetZeroTimeButton = tk.Button( + text="set t0", fg="black", + command=lambda : [self.setZeroTime(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + SetZeroTimeButton.config(height = 1, width = halfwid) + SetZeroTimeButton.grid(row=3, column=rightcol+1, sticky='nsew') + self.balloon.bind(SetZeroTimeButton, + "Set the two-way travel time \n" + "that corresponds to the surface.\n" + "/!\ Pick zero time first with button Pick") + + + + # TimeZero Adjust = align traces + TrAlignButton = tk.Button( + text="align traces", fg="black", + command=lambda : [proj.alignTraces(), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + TrAlignButton.config(height = 1, width = 2*halfwid) + TrAlignButton.grid(row=4, column=rightcol, sticky='nsew',columnspan=colsp) + self.balloon.bind(TrAlignButton, + 'Automatically shifts each trace up or down\n' + 'such that the maximum aplitudes of the individual\n' + 'traces align. Can lead to problems when the maxima\n' + 'are not in the air waves. If the results are bad,\n' + 'use the "undo" button.') + + + + # truncate Y + truncYButton = tk.Button( + text="truncate Y", fg="black", + command=lambda : [self.truncateY(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + truncYButton.config(height = 1, width = 2*halfwid) + truncYButton.grid(row=5, column=rightcol, sticky='nsew',columnspan=colsp) + self.balloon.bind(truncYButton, + "Remove data points at arrival times\n" + "later than the chosen value. If velocity\n" + "is given: remove data points at depths greater\n" + "than the chosen value") + + + + # Cut + cutButton = tk.Button( + text="cut profile", fg="black", + command=lambda : [self.cut(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + cutButton.config(height = 1, width = 2*halfwid) + cutButton.grid(row=6, column=rightcol, sticky='nsew',columnspan=colsp) + self.balloon.bind(cutButton, + "trims data to desired along-profile range.") + + + + # Dewow + DewowButton = tk.Button( + text="dewow", fg="black", + command=lambda : [self.dewow(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + DewowButton.config(height = 1, width = 2*halfwid) + DewowButton.grid(row=7, column=rightcol, sticky='nsew',columnspan=colsp) + self.balloon.bind(DewowButton, + "Trace-wise low-cut filter. Removes\n" + "from each trace a running mean of\n" + "chosen window width.") + + + # Rem mean trace + remMeanTraceButton = tk.Button( + text="rem mean tr", fg="black", + command=lambda : [self.remMeanTrace(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + remMeanTraceButton.config(height = 1, width = 2*halfwid) + remMeanTraceButton.grid(row=8, column=rightcol, sticky='nsew',columnspan=colsp) + self.balloon.bind(remMeanTraceButton, + "Removes from each trace the average\n" + "of its surrounding traces. This can be\n" + "useful to remove air waves, ground\n" + "waves, or horizontal features.") + + + # Smooth + SmoothButton = tk.Button( + text="smooth (temp)", fg="black", + command=lambda : [self.smooth(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + SmoothButton.config(height = 1, width = 2*halfwid) + SmoothButton.grid(row=9, column=rightcol, sticky='nsew',columnspan=colsp) + self.balloon.bind(SmoothButton, + "Trace-wise high-cut filter. Replaces\n" + "each sample within a trace by a\n" + "running mean of chosen window width.") + + + + + + # profile Smoothing Button + profSmButton = tk.Button( + text="profile smoothing", fg="black", + command=lambda : [self.profileSmooth(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + profSmButton.config(height = 1, width = 2*halfwid) + profSmButton.grid(row=10, column=rightcol, sticky='nsew',columnspan=colsp) + self.balloon.bind(profSmButton, + "First oversamples the profile (makes 'n' copies\n" + "of each trace) and then replaces each trace by\n" + "the mean of its neighboring 'm' traces." + ) + + + + # Gain + tpowButton = tk.Button( + text="tpow", fg="black", + command=lambda : [self.tpowGain(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + tpowButton.config(height=1, width=halfwid) + tpowButton.grid(row=11, column=rightcol, sticky='nsew') + self.balloon.bind(tpowButton, + "t-power gain. Increases the power of the\n" + "signal by a factor of (two-way travel time)^p,\n" + "where the user provides p. This gain is often\n" + "less aggressive than agc.") + + + agcButton = tk.Button( + text="agc",fg="black", + command=lambda : [self.agcGain(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + agcButton.config(height=1, width=halfwid) + agcButton.grid(row=11, column=rightcol+1, sticky='nsew') + self.balloon.bind(agcButton, + "Automatic gain controll. Normalizes the power\n" + "of the signal per given sample window along\n" + "each trace.") + + # show hyperbola + hypButton = tk.Button( + text="show hyperb", fg="black", + command=lambda : [self.showHyp(proj,a), canvas.draw()]) + hypButton.config(height = 1, width = 2*halfwid) + hypButton.grid(row=12, column=rightcol, sticky='nsew',columnspan=colsp) + self.balloon.bind(hypButton, + "Draws a hyperbola depending on profile position,\n" + "two-way travel time, and estimated velocity. This\n" + "can be used to find the subsurface velocity when\n" + "a hyperbola is visible in the data.\n" + "The plotted hyperbola will disappear when the image\n" + "is refreshed.") + + + + # Set Velocity + setVelButton = tk.Button( + text="set vel", fg="black", + command=lambda : [self.setVelocity(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + setVelButton.config(height = 1, width = halfwid) + setVelButton.grid(row=13, column=rightcol, sticky='nsew',columnspan=1) + self.balloon.bind(setVelButton, + "Set the known subsurface radar velocity. This will\n" + "turn the y-axis from two-way travel time to depth.\n" + "This step is necessary for topographic correction.") + + # Correct for antenna separation + antennaSepButton = tk.Button( + text="ant sep", fg="black", + command=lambda : [self.antennaSep(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + antennaSepButton.config(height = 1, width = halfwid) + antennaSepButton.grid(row=13, column=rightcol+1, sticky='nsew',columnspan=1) + self.balloon.bind(antennaSepButton, + "If the antenna offset is provided, this corrects for\n" + "distortion of arrival times near the surface due to\n" + "the separation of transmitter and receiver antenna." + "You must have picked the first break of the airwave\n" + "for this to function properly and the velocity must be set.") + + # pickFrequenciesButton + pickFrequenciesButton = tk.Button( + text="Pick f", fg="black", + command=lambda : [self.pickFrequencies(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + pickFrequenciesButton.config(height = 1, width = halfwid) + pickFrequenciesButton.grid(row=14, column=rightcol, sticky='nsew') + self.balloon.bind(pickFrequenciesButton, + "pick frequencies from one FFT trace\n" + "Created by A.E \n") + + + # filterButton + filterButton = tk.Button( + text="bandpass", fg="black", + command=lambda : [self.bandpass(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + filterButton.config(height = 1, width = halfwid) + filterButton.grid(row=14, column=rightcol+1, sticky='nsew') + self.balloon.bind(filterButton, + "Bandpass filter\n" + "/!\ Pick frequencies first with button Pick\n" + "/!\ CLICK on figure to close it, DO NOT CLOSE IT MANUALLY \n" + "Created by A.R \n" + "Modified by A.E \n") + # # Migration Button + # migButton = tk.Button( + # text="fk migration", fg="black", + # command=lambda : [self.fkMigration(proj), + # self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + # migButton.config(height = 1, width = 2*halfwid) + # migButton.grid(row=14, column=rightcol, sticky='nsew',columnspan=colsp) + # self.balloon.bind(migButton, + # "Stolt fk migration using a code originally written\n" + # "in Matlab for the CREWES software package.\n" + # "Translated into Python 2 by Nat Wilson.\n" + # "\n" + # "Not included in the public version because of License\n" + # "uncertainty. Contact alainplattner@gmail.com\n" + # "if you would like to use it.") + + + + + # Topo Correct + topoCorrectButton = tk.Button( + text="topo correct", fg="black", + command=lambda : [self.topoCorrect(proj), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + topoCorrectButton.config(height = 1, width = 2*halfwid) + topoCorrectButton.grid(row=15, column=rightcol, sticky='nsew',columnspan=colsp) + self.balloon.bind(topoCorrectButton, + "Reads a comma- or tab-separated file containing\n" + "either 3 columns (easting, northing, elevation)\n" + "or two columns (profile position, elevation).\n" + "All coordinates in meters.") + + + + startPickButton = tk.Button( + text="start pick", fg="black", + command=lambda : self.startPicking(proj,fig=fig,a=a,canvas=canvas)) + startPickButton.config(height = 1, width = halfwid) + startPickButton.grid(row=16, column=rightcol, sticky='nsew',columnspan=1) + self.balloon.bind(startPickButton, + "Start collecting location information\n" + "by clicking on the profile.") + + + stopPickButton = tk.Button( + text="stop pick", fg="black", + command=lambda : [self.stopPicking(proj,canvas), + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas)]) + + stopPickButton.config(height = 1, width = halfwid) + stopPickButton.grid(row=16, column=rightcol+1, sticky='nsew',columnspan=1) + self.balloon.bind(stopPickButton, + "Stop collecting location information\n" + "and save the locations you collected\n" + "in a text file.") + + # Save data + SaveButton = tk.Button( + text="save data", fg="black", + command=lambda : self.saveData(proj)) + SaveButton.config(height = 1, width = 2*halfwid) + SaveButton.grid(row=17, column=rightcol, sticky='nsew',columnspan=colsp) + self.balloon.bind(SaveButton, + 'saves the processed data including its history in a\n' + '.gpr file. The resulting file will contain absolute\n' + 'path names of the used data and topography files.\n' + 'Visualization settings such as "set x-range" or\n' + '"contrast" will not be saved.') + + + + # Print Figure + PrintButton = tk.Button( + text="print figure", fg="black", + command=lambda : self.printProfileFig(proj=proj,fig=fig)) + PrintButton.config(height = 1, width = 2*halfwid) + PrintButton.grid(row=18, column=rightcol, sticky='nsew',columnspan=colsp) + self.balloon.bind(PrintButton, + "Saves the current visible figure in a pdf with \n" + "chosen resolution. If there is a hyperbola on\n" + "the current figure, then the hyperbola will also\n" + "appear on the printed figure.") + + + # Export to VTK + VTKButton = tk.Button( + text="export to VTK", fg="black", + command = lambda : self.exportVTK(proj)) + VTKButton.config(height = 1, width = 2*halfwid) + VTKButton.grid(row=19, column=rightcol, sticky='nsew',columnspan=colsp) + self.balloon.bind(VTKButton, + "Exports the processed figure to a\n" + "VTK format, that can be read by\n" + "Paraview or similar 3D programs") + + + + + # Write script + HistButton = tk.Button( + text="write script", fg="black", + command=lambda : self.writeHistory(proj)) + HistButton.config(height = 1, width = 2*halfwid) + HistButton.grid(row=20, column=rightcol, sticky='nsew',columnspan=colsp) + self.balloon.bind(HistButton, + 'Writes a python script to reproduce the \n' + 'current status.\n' + '\n' + 'If the current data is from a .gpr file, \n' + 'then the python script will contain all \n' + 'steps going back to the raw data. \n' + '\n' + 'The script will not contain visualization \n' + 'settings such as x-range settings, unless \n' + 'the "print figure" command was used. ') + + + + + def undo(self,proj): + if self.picking: + self.picked=self.picked[0:-1,:] + else: + proj.undo() + + + def setYrng(self): + ylow = sd.askfloat("Input","Min Y value") + if ylow is not None: + yhigh = sd.askfloat("Input","Max Y value") + if yhigh is not None: + self.prevyrng=self.yrng + self.yrng=[ylow,yhigh] + + + def resetYrng(self,proj): + # Only needed in undo, and only if what you want to + # undo changed the y axis + if ("setVelocity" in proj.history[-1]) or ("topoCorrect" in proj.history[-1]): + self.yrng=self.prevyrng + + + def setAspect(self): + self.asp = sd.askfloat("Input","Plotting aspect ratio") + + + def setFullView(self,proj): + self.xrng=[np.min(proj.profilePos),np.max(proj.profilePos)] + if proj.velocity is None: + self.yrng=[np.min(proj.twtt),np.max(proj.twtt)] + elif proj.maxTopo is None: + self.yrng=[np.min(proj.depth),np.max(proj.depth)] + else: + self.yrng=[proj.minTopo-np.max(proj.depth),proj.maxTopo-np.min(proj.depth)] + + + def toggleGrid(self): + self.grid = not self.grid + + + def setXrng(self): + xlow = sd.askfloat("Input","Min X value") + if xlow is not None: + xhigh = sd.askfloat("Input","Max X value") + if xhigh is not None: + self.xrng=[xlow,xhigh] + + + def adjProfile(self,proj): + flipit = mesbox.askyesno("Question","Flip the profile (left to right)?") + if flipit: + proj.flipProfile() + minPos = sd.askfloat("Input","Start x coordinate") + if minPos is not None: + maxPos = sd.askfloat("Input","End x coordinate") + if maxPos is not None: + proj.adjProfile(minPos=minPos,maxPos=maxPos) + self.xrng=[minPos,maxPos] + + def pickZeroTime(self,proj): + trace_pour_afficher = sd.askfloat("Input","Number of trace to print") + proj.pick_time_zero(trace_pour_afficher) + + def setZeroTime(self,proj): + #newZeroTime = sd.askfloat("Input","New zero time") + #if newZeroTime is not None: + # proj.setZeroTime(newZeroTime=newZeroTime) + # trace_pour_afficher = sd.askfloat("Input","Trace number to print") + # proj.pick_time_zero(trace_pour_afficher) + # newZeroTime = sd.askfloat("Input","New zero time") + proj.setZeroTime() + + def bandpass(self, proj): + trace_pour_afficher = sd.askfloat("Input","Number of trace to print") + proj.bandpass(trace_pour_afficher) + + def pickFrequencies(self,proj): + trace_pour_afficher = sd.askfloat("Input","Number of trace to print") + proj.pick_frequencies(trace_pour_afficher) + + def dewow(self,proj): + window = sd.askinteger("Input","Dewow window width (number of samples)") + if window is not None: + proj.dewow(window=window) + + + def smooth(self,proj): + window = sd.askinteger("Input","Smoothing window width (number of samples)") + if window is not None: + proj.smooth(window=window) + + + def remMeanTrace(self,proj): + ntraces = sd.askinteger("Input","Remove mean over how many traces?") + if ntraces is not None: + proj.remMeanTrace(ntraces=ntraces) + + + def tpowGain(self,proj): + power = sd.askfloat("Input","Power for tpow gain?") + if power is not None: + proj.tpowGain(power=power) + + + def agcGain(self,proj): + window = sd.askinteger("Input","Window length for AGC?") + if window is not None: + proj.agcGain(window=window) + + def truncateY(self,proj): + maxY = sd.askfloat("Input","Truncate at what y value\n" + "(two-way travel time or depth)") + if maxY is not None: + proj.truncateY(maxY) + + def cut(self,proj): + minX = sd.askfloat("Input","Minimum profile position") + if minX is not None: + maxX = sd.askfloat("Input","Maximum profile position") + if maxX is not None: + proj.cut(minX,maxX) + + def setVelocity(self,proj): + velocity = sd.askfloat("Input","Radar wave velocity [m/ns]?") + if velocity is not None: + proj.setVelocity(velocity) + self.prevyrng=self.yrng + self.yrng=[0,np.max(proj.depth)] + + def antennaSep(self,proj): + if proj.velocity is None: + mesbox.showinfo("Antenna Sep Error","You have to set the velocity first") + proj.antennaSep() + + + def fkMigration(self,proj): + if proj.velocity is None: + mesbox.showinfo("Migration Error","You have to set the velocity first") + proj.fkMigration() + + + def profileSmooth(self,proj): + ntraces = sd.askinteger("Input","Smooth over how many traces (m)") + if ntraces is not None: + noversample = sd.askinteger("Input","Make how many copies of each trace (n).\nRecommended: Same as number of traces to be smoothed.") + if noversample is not None: + proj.profileSmooth(ntraces,noversample) + + + def topoCorrect(self,proj): + if proj.velocity is None: + mesbox.showinfo("Topo Correct Error","You have to set the velocity first") + return + topofile = fd.askopenfilename() + if topofile is not '': + out = self.getDelimiter() + proj.topoCorrect(topofile,self.delimiter) + self.prevyrng=self.yrng + self.yrng=[proj.minTopo-np.max(proj.depth),proj.maxTopo] + + + + + def startPicking(self,proj,fig,a,canvas): + self.picking = True + self.picked = np.asmatrix(np.empty((0,2))) + print("Picking mode on") + def addPoint(event): + self.picked = np.append(self.picked,np.asmatrix([event.xdata,event.ydata]),axis=0) + self.plotProfileData(proj,fig=fig,a=a,canvas=canvas) + print(self.picked) + self.pick_cid = canvas.mpl_connect('button_press_event', addPoint) + + + + def stopPicking(self,proj,canvas): + filename = fd.asksaveasfilename() + if filename is not '': + self.picking = False + canvas.mpl_disconnect(self.pick_cid) + print("Picking mode off") + np.savetxt(filename+'_profile.txt',self.picked,delimiter='\t') + print('saved picked file as "%s"' %(filename+'_profile.txt')) + # If we have 3D info, also plot it as 3D points + if proj.threeD is not None: + # First calculate along-track points + topoVal = proj.threeD[:,2] + npos = proj.threeD.shape[0] + steplen = np.sqrt( + np.power( proj.threeD[1:npos,0]-proj.threeD[0:npos-1,0] ,2.0) + + np.power( proj.threeD[1:npos,1]-proj.threeD[0:npos-1,1] ,2.0) + + np.power( proj.threeD[1:npos,2]-proj.threeD[0:npos-1,2] ,2.0) + ) + alongdist = np.cumsum(steplen) + topoPos = np.append(0,alongdist) + pick3D = np.zeros((self.picked.shape[0],3)) + # If profile is adjusted, need to start the picked at zero. + pickProfileShifted = self.picked[:,0] - np.min(proj.profilePos) + #for i in range(0,3): + for i in range(0,2): + pick3D[:,i] = interp.pchip_interpolate(topoPos, + proj.threeD[:,i], + pickProfileShifted).squeeze() + #self.picked[:,0]).squeeze() + + pick3D[:,2] = self.picked[:,1].squeeze() + + np.savetxt(filename+'_3D.txt',pick3D,delimiter='\t') + print('saved picked file as "%s"' %(filename+'_3D.txt')) + + + def loadData(self,proj): + filename = fd.askopenfilename( filetypes= (("All", "*.*"), + ("GPRPy (.gpr)", "*.gpr"), + ("Sensors and Software (.DT1)", "*.DT1"), + ("GSSI (.DZT)", "*.DZT"), + ("BSQ header","*.GPRhdr"), + ("MALA header","*.rad"))) + if filename: + marker_existance = sd.askstring("'yes' or 'no'","Marker file exists ?") + + if marker_existance == 'yes' : + marker = fd.askopenfilename() #ajout aline : fichier marqueur pour la position + proj.importdata(filename=filename, marker = marker) + else : + proj.importdata(filename=filename) + + self.xrng = [np.min(proj.profilePos),np.max(proj.profilePos)] + if proj.depth is None: + self.yrng = [0,np.max(proj.twtt)] + else: + if proj.maxTopo is None: + self.yrng = [0,np.max(proj.depth)] + else: + self.yrng = [proj.maxTopo-np.max(proj.depth), proj.maxTopo] + self.asp=None + # Just in case someone presses undo before changing yrange + self.prevyrng=self.yrng + print("Loaded " + filename) + + + + + def saveData(self,proj): + filename = fd.asksaveasfilename(defaultextension=".gpr") + if filename is not '': + proj.save(filename) + + + def exportVTK(self,proj): + outfile = fd.asksaveasfilename() + if outfile is not '': + #thickness = sd.askfloat("Input","Profile thickness [m]") + thickness = 0 + if self.asp is None: + aspect = 1.0 + else: + aspect = self.asp + + if proj.threeD is None: + gpyes = mesbox.askyesno("Question","Do you have topography data for this profile?") + if gpyes: + filename = fd.askopenfilename() + self.getDelimiter() + proj.exportVTK(outfile,gpsinfo=filename,thickness=thickness,delimiter=self.delimiter,aspect=aspect) + else: + proj.exportVTK(outfile,gpsinfo=proj.threeD,thickness=thickness,delimiter=self.delimiter,aspect=aspect) + print('... done with exporting to VTK.') + + def writeHistory(self,proj): + filename = fd.asksaveasfilename(defaultextension=".py") + if filename is not '': + proj.writeHistory(filename) + print("Wrote script to " + filename) + + + def plotProfileData(self,proj,fig,a,canvas): + # Clear cursor coordinate cid if if exists to avoid multiple instances + if 'self.cursor_cid' in locals(): + canvas.mpl_disconnect(self.cursor_cid) + dx=proj.profilePos[3]-proj.profilePos[2] + dt=proj.twtt[3]-proj.twtt[2] + a.clear() + stdcont = np.nanmax(np.abs(proj.data)[:]) + if proj.velocity is None: + a.imshow(proj.data,cmap=self.color.get(),extent=[min(proj.profilePos)-dx/2.0, + max(proj.profilePos)+dx/2.0, + max(proj.twtt)+dt/2.0, + min(proj.twtt)-dt/2.0], + aspect="auto", + vmin=-stdcont/self.contrast.get(), vmax=stdcont/self.contrast.get()) + a.set_ylim(self.yrng) + a.set_xlim(self.xrng) + a.set_ylabel("two-way travel time [ns]", fontsize=mpl.rcParams['font.size']) + a.invert_yaxis() + elif proj.maxTopo is None: + dy=dt*proj.velocity + a.imshow(proj.data,cmap=self.color.get(),extent=[min(proj.profilePos)-dx/2.0, + max(proj.profilePos)+dx/2.0, + max(proj.depth)+dy/2.0, + min(proj.depth)-dy/2.0], + aspect="auto", + vmin=-stdcont/self.contrast.get(), vmax=stdcont/self.contrast.get()) + a.set_ylabel("depth [m]", fontsize=mpl.rcParams['font.size']) + a.set_ylim(self.yrng) + a.set_xlim(self.xrng) + a.invert_yaxis() + else: + dy=dt*proj.velocity + a.imshow(proj.data,cmap=self.color.get(),extent=[min(proj.profilePos)-dx/2.0, + max(proj.profilePos)+dx/2.0, + proj.minTopo-max(proj.depth)-dy/2.0, + proj.maxTopo-min(proj.depth)+dy/2.0], + aspect="auto", + vmin=-stdcont/self.contrast.get(), vmax=stdcont/self.contrast.get()) + a.set_ylabel("elevation [m]", fontsize=mpl.rcParams['font.size']) + a.set_ylim(self.yrng) + a.set_xlim(self.xrng) + + a.get_xaxis().set_visible(True) + a.get_yaxis().set_visible(True) + a.set_xlabel("profile position [m]", fontsize=mpl.rcParams['font.size']) + a.xaxis.tick_top() + a.xaxis.set_label_position('top') + if self.asp is not None: + a.set_aspect(self.asp) + + # Set grid + a.grid(self.grid) + + # In case you are picking + if self.picking: + a.plot(self.picked[:,0],self.picked[:,1],'-x',color='yellow',linewidth=3*self.highfac) + a.plot(self.picked[:,0],self.picked[:,1],'-x',color='black',linewidth=2*self.highfac) + + # Allow for cursor coordinates being displayed + def moved(event): + if event.xdata is not None and event.ydata is not None: + canvas.get_tk_widget().itemconfigure(tag, text="(x = %5.5g, y = %5.5g)" % (event.xdata, event.ydata)) + + self.cursor_cid = canvas.mpl_connect('button_press_event', moved) + tag = canvas.get_tk_widget().create_text(20, 20, text="", anchor="nw") + + canvas.get_tk_widget().grid(row=2,column=0,columnspan=figcolsp, rowspan=figrowsp, sticky='nsew') + canvas.draw() + + + # Show hyperbola + def showHyp(self,proj,a): + x0 = sd.askfloat("Input","Hyperbola center on profile [m]") + if x0 is not None: + t0 = sd.askfloat("Input","Hyperbola apex location (two-way travel time [ns])") + if t0 is not None: + v = sd.askfloat("Input","Estimated velocity [m/ns]") + if v is not None: + y=proj.profilePos-x0 + d=v*t0/2.0 + k=np.sqrt(d**2 + np.power(y,2)) + t2=2*k/v + a.plot(proj.profilePos,t2,'--c',linewidth=3) + + + def printProfileFig(self,proj,fig): + figname = fd.asksaveasfilename(defaultextension=".pdf") + if figname is not '': + dpi = sd.askinteger("Input","Resolution in dots per inch? (Recommended: 600)") + if dpi is not None: + fig.savefig(figname, format='pdf', dpi=dpi) + # Put what you did in history + if self.asp is None: + histstr = "mygpr.printProfile('%s', color='%s', contrast=%g, yrng=[%g,%g], xrng=[%g,%g], dpi=%d)" %(figname,self.color.get(),self.contrast.get(),self.yrng[0],self.yrng[1],self.xrng[0],self.xrng[1],dpi) + else: + histstr = "mygpr.printProfile('%s', color='%s', contrast=%g, yrng=[%g,%g], xrng=[%g,%g], asp=%g, dpi=%d)" %(figname,self.color.get(),self.contrast.get(),self.yrng[0],self.yrng[1],self.xrng[0],self.xrng[1],self.asp,dpi) + proj.history.append(histstr) + print("Saved figure as %s" %(figname+'.pdf')) + + + + def getDelimiter(self): + commaQuery = tk.Toplevel(self.window) + commaQuery.title("Comma or tab separated?") + text = tk.Label(commaQuery,text="Is this a comma- or tab-separated file?",fg='red') + text.pack(padx=10,pady=10) + commaButton = tk.Button(commaQuery,text="comma",width=10, + command = lambda: [self.setComma(), + commaQuery.destroy()]) + commaButton.pack(side="left") + tabButton = tk.Button(commaQuery,text="tab",width=10, + command = lambda: [self.setTab(), + commaQuery.destroy()]) + tabButton.pack(side="right") + #self.window.frame().tansient(self.window) + #self.window.frame().grab_set() + self.window.wait_window(commaQuery) + def setComma(self): + self.delimiter = ',' + print("Delimiter set to comma") + def setTab(self): + self.delimiter = '\t' + print("Delimiter set to tab") + + +# root = tk.Tk() + +# for col in range(rightcol): +# root.columnconfigure(col, weight=1) +# for row in range(figrowsp): +# root.rowconfigure(row, weight=1) + +# app = GPRPyApp(root) + +# root.mainloop() + + + diff --git a/readmrq.py b/readmrq.py new file mode 100644 index 0000000..afabff8 --- /dev/null +++ b/readmrq.py @@ -0,0 +1,57 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Jul 1 16:34:37 2020 + +@author: aline +""" +import numpy as np +import matplotlib.pyplot as plt +from scipy import interpolate + + +def read_mrk(filename, nb_trace): + """ + @author : aline + + Parameters + ---------- + filename : TYPE string + DESCRIPTION. path+filename du fichier marqueur à lire + delta_x : TYPE float ou int + DESCRIPTION. pas en x (m) + nb_trace : TYPE int + DESCRIPTION. nombre de traces + + Returns + ------- + TYPE list + DESCRIPTION. positions (en m) pour chaque trace + + """ + + + fid = open(filename, 'r') + position = {} + xp = [] #liste des numeros de trace + yp = [] #liste des positions (en m) + + for line in fid : + num_trace = line.split()[0] + pos = line.split()[1] + + xp.append(int(num_trace)) + yp.append(int(pos)) + + fid.close() + + delta_x = (yp[-1]-yp[0])/(xp[-1]-xp[0]) + + + f = interpolate.interp1d(xp, yp, fill_value="extrapolate") + traces = np.linspace(0,nb_trace,nb_trace) + position = f(traces) + #plt.plot(position,traces) + #print('position',position[:nb_trace]) + + return position[:nb_trace] \ No newline at end of file