#################################################
#################################################
#################################################
###Python code used throughout thesis titled:
###'Electrical resistivity and surface-wave methods
###for the detection of shallow objects in 
###glaciolacustrine clay'
###by Ian Robert Deniset
###Completed in 2020
#################################################
#################################################
#################################################

#################################################
#################################################
###STANDARD MODULES
#################################################
#numerical
import numpy as np
from scipy.interpolate import griddata 

#plotting
import matplotlib.pyplot as plt
from matplotlib import gridspec
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from matplotlib.ticker import (MultipleLocator, FormatStrFormatter, AutoMinorLocator)
import matplotlib.colors
from matplotlib import cm

#file manipulation
import os
from pathlib import Path
import sys

#suppress warnings
import warnings
warnings.filterwarnings('ignore')


#################################################
#################################################
###SEISMIC CODE
#################################################

##MODULES for SEISMIC
from obspy import read

##################
#FILE LOADING CODE
##################

##################
#Load modeled seismic gathers

def loadModelDAT(file):
	'''
	Takes a ParkSeis '.dat' file generated during modeling and returns an Obspy stream object and raw traces
	'''

	#load in raw .dat file
	strm = read(file)
	
	#generator function here to save memory for big files; basically list comprehension 
	trcs = np.stack(t.data for t in strm.traces)
	
	return trcs, strm


##################
#Load seismic data exported from ParkSEIS as a TXT file

def loadParkSeisDataTEXT(file):
	'''
	Takes in a seismic data file exported from ParkSEIS as a .TXT file and returns
	model parameters and traces loaded into a numpy array of shape (trNum, sampNum)
	'''
	
	#load in the text file using np.genfromtext
	tr, t, amp = np.genfromtxt(file, unpack=True)
	
	#get model parameters and load into dict
	nTrace = int(np.unique(tr)[-1])
	npts = int(tr.size/nTrace)
	dt = t[5]-t[4]
	
	modParamsTXT = {
		'numTrc' : nTrace,
		'numSamps' : npts,
		'dt' : dt
	}
	
	#load amplitudes into array to return as 'traces'
	txtTraces = np.full((nTrace, npts), np.nan) #create empty array of proper shape and fill with np.nan values
	
	for trace in np.unique(tr): #load into created empty array
		trIdx = int(trace)-1
		txtTraces[trIdx] = amp[tr==trace]
		
	return modParamsTXT, txtTraces


##################
#Load a dispersion curve

def loadTheorDC(fname, asDict=False, dictKey=''):
	'''
	Take in a file path to a ParkSeis dispersion curve and extracts the curve data into
	an array or dictionary
	'''

	numPts, modeNum = getTheorDCParams(fname) #read file and get info

	M = np.genfromtxt(fname, skip_header=1, usecols=(1,2), max_rows=numPts)

	if asDict == True: #for loading singular experimental curves from inversions
		mDict = {} 
		mDict[dictKey] = M.T
		return mDict
	else:
		return M.T


##################
#Load all dispersion curves for a single model at once

def loadTheorDC_PATH(path):
	'''
	Takes a path to a folder containing all theoretical dispersion curves generated for a given model in ParkSeis and returns
	a dictionary 
	'''

	dcCurves = {}

	for root, dirs, files in os.walk(path):
		fNum = 0
		for f in files:

			#this method allows for more than 9 modes and also loading in AM0 curves
			if f.endswith('.DC'):
				key = f.split('(')[-1].split(')')[0] #this is to get the mode number label from the bracket value in file name

				dcCurves[key] = loadTheorDC(os.path.join(path, f))

	return dcCurves


##################
#Read a *.LYR and return parameters

def readLYRFile(fname):
	''''
	Takes a layer file and returns all parameter values for all layers

	For use in plotting of the layer file
	'''
	#get number of layers in model; note returned value is number of distinct layers PLUS one (i.e. includes half-space)
	with open(fname, 'r') as lyrFile:
		lines = lyrFile.readlines()
		numLayers = int(lines[2][14:]) #number of layers is on 3rd line starting after the 14 col
		paramLines = lines[3:(3+numLayers)] #gets all lines containing model parameters
	
	
	#load in the model parameters
	z, h, vs, vp, pr, rho, qs, qp, vsU, vsL = np.genfromtxt(fname, skip_header=4, usecols=(1,2,3,4,5,6,7,8,9,10), 
												  max_rows=numLayers, unpack=True)
	
	#calculate parameters needed
	hsADD = 1.5 #depth value to add to half-space for plotting
	hsDepth = z[numLayers-2] #get depth to half-space
	z[-1] = hsDepth + hsADD #sets depth of half-space base
	h[-1] = z[-1] - z[-2] #sets thickness of half-space
	
	dz = 0.0001
	zProf = np.arange(0, hsDepth+hsADD, dz) #calculate depth vector
	vsProf = np.full_like(zProf, np.nan)
	vpProf = np.full_like(zProf, np.nan)
	prProf = np.full_like(zProf, np.nan)
	rhoProf = np.full_like(zProf, np.nan)
	qsProf = np.full_like(zProf, np.nan)
	qpProf = np.full_like(zProf, np.nan)
	vsUProf = np.full_like(zProf, np.nan)
	vsLProf = np.full_like(zProf, np.nan)
	
	for layer in range(numLayers):
		upper = z[layer] - h[layer]
		lower = z[layer]
		
		vsProf[(zProf >= upper) & (zProf <= lower)] = vs[layer]
		vpProf[(zProf >= upper) & (zProf <= lower)] = vp[layer] 
		prProf[(zProf >= upper) & (zProf <= lower)] = pr[layer] 
		rhoProf[(zProf >= upper) & (zProf <= lower)] = rho[layer]
		qsProf[(zProf >= upper) & (zProf <= lower)] = qs[layer] 
		qpProf[(zProf >= upper) & (zProf <= lower)] = qp[layer]
		vsUProf[(zProf >= upper) & (zProf <= lower)] = vsU[layer]
		vsLProf[(zProf >= upper) & (zProf <= lower)] = vsL[layer]

	layerFileValues = {
	'Depth':zProf,
	'Vs':vsProf,
	'Vp':vpProf,
	'Poissons Ratio':prProf,
	'Density':rhoProf,
	'Qs':qsProf,
	'Qp':qpProf,
	'Vs Upper':vsUProf,
	'Vs Lower':vsLProf
	}

	return layerFileValues


##################
#PLOTTING CODE
##################

##################
#Plot modeled gather as wiggles 

def plotModelData(data, startTime, endTime, X1, dx, percentile=99.99, gain=1, traceNorm=None, title='Shot Record', ax=None):
	'''
	Takes either a model file or Obspy stream and plots the data

	 X1: Offset from source to first receiver
	 dx: receiver spacing
	'''

	if data.endswith('.TXT'):
		#case if model file is exported as TXT file
		modParams, traces = loadParkSeisDataTEXT(data)

	elif isinstance(data, str):
		#case if data is a file path; loads stream and parameters then stacks traces into array
		traces, stream = loadModelDAT(data)
		modParams = getModelParams(stream)

	else:
		#case if data is a Obspy stream object
		traces = np.stack(t.data for t in data.traces)
		modParams = getModelParams(data)

	#parameters for plotting
	perc = np.percentile(traces, percentile)
	gain = gain
	nTrace = modParams['numTrc']
	npts = modParams['numSamps']
	dt = modParams['dt']
	dx = dx
	srcOffset = X1
	start = int(startTime/dt)
	end = int(endTime/dt)

	xOffset = np.arange(0, (nTrace)*dx, dx)
	t = np.arange(0, npts*dt, dt)


	#plot
	if ax != None:
		ax = ax
		msg = 'Master Plot'
	else:
		fig, ax = plt.subplots(1,1, figsize=(15,15))
		msg = 'No Master Plot'

	trNum = 1

	for xpos, tr in zip(xOffset, traces):
		#optional normalization
		if traceNorm == 'Per Trace': amp = gain * tr / np.nanmax(tr) + (xpos+srcOffset) 
		else: amp = gain * tr / perc + (xpos+srcOffset)

		ax.plot(amp[start:end], t[start:end], c='k', lw=0.5)
		ax.fill_betweenx(t[start:end], amp[start:end], 
			(xpos+srcOffset), amp[start:end] > (xpos+srcOffset), color='r', lw=0, alpha=0.7)
		ax.fill_betweenx(t[start:end], amp[start:end], 
			(xpos+srcOffset), amp[start:end] < (xpos+srcOffset), color='blue', lw=0, alpha=0.7)
		ax.text(xpos+srcOffset, 0, str(trNum), horizontalalignment='center', fontsize=9)

		trNum += 1

	ax.set_ylabel('Time [s]')
	ax.set_xlabel('Offset from Source [m]')
	ax.set_title(title)
	ax.invert_yaxis()
	if msg == 'No Master Plot':
		plt.show()


	return

##################
#Calculate and plot dispersion image for (modeled) gather

def dcModelPhaseShift(data, minVel, maxVel, minFreq, maxFreq, X1, dx, dv=0.1, padLen=None, title='Dispersion Image', 
	overLayDC=None, ax=None):
	'''Takes either a gather file or an Obspy stream along with min/max velocity values to check
	and min/max frequency values to check between; output is a dispersion curve image.
	
	After Park et al. 1998

	dv: velocity step used in image calculation
	padLen: number of samples to pad spatial axis with
	overLayDC: dictionary containing dispersion curves to plot on top of image
	'''
	
	################
	#Calculate Image
	################

	if isinstance(data, str):
		#load in file to get raw data matrix and stream then extract parameters
		gather, gatherStream = loadModelDAT(data)
		params = getModelParams(gatherStream)
		
	else:
		#extract traces from stream and get acquisition parameters
		gather = np.stack(t.data for t in data.traces)
		params = getModelParams(data)
		

	#pad spatial axis if specified
	if padLen != None:
		gather = np.pad(gather, ((0,0),(padLen,padLen)), 'constant', constant_values=(0,0))

	#calculate offset vector
	xx = np.arange(X1, X1 + params['numTrc']*dx, dx)
	
	#compute fft and associated freqs
	G = np.fft.fft(gather)
	freqs = np.fft.fftfreq(gather[0].size, params['dt'])
	
	#select only positive frequencies from fft ouput; i.e. first half
	Gp = G[:,:freqs.size//2]
	freqsp = freqs[:freqs.size//2]
	
	#select frequencies
	df = freqs[10]-freqs[9]
	fMax = int(maxFreq/df)
	fMin = int(minFreq/df)

	
	#set up velocities to test
	testVels = np.arange(minVel, maxVel, dv)
	
	#create empty array to hold transformed data and mode picks
	V = np.zeros((len(freqsp[fMin:fMax]), len(testVels)))
	numF = V.shape[0]
	M0 = np.zeros((len(freqsp[fMin:fMax])))
	
	######TRANSFORM
	#run through freqs first
	for f in range(len(freqsp[fMin:fMax])):
		#then run through each test velocity
		for v in range(len(testVels)):
			V[f,v] = np.abs( np.sum( Gp[:,f+fMin]/np.abs(Gp[:,f+fMin]) * np.exp(1j*2*np.pi*freqsp[f+fMin]*xx /testVels[v]) ) )

	#normalize by the number of traces in the gather (as suggested by Olafsdottir 2018)
	Vnorm = V/params['numTrc']  

	#calculate summation for all frequencies along single velocity value
	vSum = Vnorm.sum(axis=0) / numF  #normalize by number of frequency values
	vSum_array = np.asarray([testVels, vSum])      
		
	
	###########
	#PLOT
	###########
	if ax != None:
		ax = ax
		msg = 'Master Plot'
	else:
		fig, ax = plt.subplots(1,1, figsize=(15,8))
		msg = 'No Master Plot'


	majorXLocator = MultipleLocator(10) #sets the major tick interval to 10
	majorXFormatter = FormatStrFormatter('%d')
	minorXLocator = MultipleLocator(5) #sets the minor tick interval to every 5
	majorYLocator = MultipleLocator(50) #sets the major tick interval to every 50
	majorYFormatter = FormatStrFormatter('%d')
	minorYLocator = MultipleLocator(25) #sets the minor tick interval to every 25

	#plot theoretical dispersion curve if specified; must be first to not override extents
	if overLayDC != None:
		#loop through dictionary of DC curves and plot
		for dc in overLayDC:
			if overLayDC[dc][0].size > 1: #ensures mode curve exists in file
				dcF, dcC = overLayDC[dc][0,:], overLayDC[dc][1,:] #get mode frequency and phase vels
				if dc == 'Experimental DC':
					ax.plot(dcF, dcC, lw=3, c='m', label=dc)
				elif dc == 'Inverted Model DC':
					ax.plot(dcF, dcC, lw=3, label=dc, ls=':', c='k')
				else:
					ax.plot(dcF, dcC, lw=3, alpha=1, label=dc)

	#plot calculated dispersion image
	dispC = ax.imshow(Vnorm.T, aspect='auto', interpolation='none', extent=[fMin*df,fMax*df,maxVel,minVel])
	ax.invert_yaxis()
	ax.set_xlabel('Frequency [Hz]')
	ax.set_ylabel('Phase Velocity [m/s]')
	ax.xaxis.set_major_locator(majorXLocator)
	ax.xaxis.set_major_formatter(majorXFormatter)
	ax.xaxis.set_minor_locator(minorXLocator)
	ax.yaxis.set_major_locator(majorYLocator)
	ax.yaxis.set_major_formatter(majorYFormatter)
	ax.yaxis.set_minor_locator(minorYLocator)
	ax.set_title(title)

				
	ax.legend()
	ax.grid(which='both', linestyle='--', alpha=0.75)
	
	if msg == 'No Master Plot':
		fig.colorbar(dispC, ax = ax, shrink = 0.7)
		plt.show()

	return


##################
#Plot a MODEL *.LYR file generated in ParkSEIS

def plotLYRFile(fname, title='Layer Model Values', ax=None):
	'''
	Takes in a path to a ParkSeis layer file and outputs a visual plot of the model parameters
	'''
	
	lyrFileVals = readLYRFile(fname)
	z = lyrFileVals['Depth']
	vs = lyrFileVals['Vs']
	vp = lyrFileVals['Vp']
	rho = lyrFileVals['Density']

	#PLOT
	if ax != None:
		ax = ax
		msg = 'Master Plot'
	else:
		fig, ax = plt.subplots(1,3, figsize=(6,10))
		msg = 'No Master Plot'


	majorYLocator = MultipleLocator(1) #sets the major tick interval to every meter
	majorYFormatter = FormatStrFormatter('%d')
	minorYLocator = MultipleLocator(0.25) #sets the minor tick interval to every 0.25 meter

	ax[0].plot(vs, z, c='r', label='Vs')
	ax[0].set_xlabel('Vs [m/s]')
	ax[0].set_ylabel('Depth [m]')
	ax[0].locator_params(axis='x', tight=True, nbins=12)
	
	ax[1].plot(vp, z, c='b', label='Vp')
	ax[1].set_xlabel('Vp [m/s]')
	ax[1].locator_params(axis='x', tight=True, nbins=8)
	
	ax[2].plot(rho, z, c='k', label='Density')
	ax[2].set_xlabel('Density [kg/m$^3$]')
	ax[2].set_xlim(1,3)
	ax[2].locator_params(axis='x', tight=True, nbins=4)

	for i in range(len(ax)):
		ax[i].invert_yaxis()
		ax[i].yaxis.set_major_locator(majorYLocator)
		ax[i].yaxis.set_major_formatter(majorYFormatter)
		ax[i].yaxis.set_minor_locator(minorYLocator)
		ax[i].grid(which='both', linestyle='--', alpha=0.75)
		ax[i].tick_params(axis='x', labelrotation=90)
		if i > 0:
			ax[i].set_yticklabels([])
			ax[i].tick_params(axis='y', which='both', left='off')
	
	
	if msg == 'No Master Plot':
		fig.suptitle(title, y=0.92)
		plt.subplots_adjust(wspace=0.05, hspace=0)
		plt.show()

	
	return


##################
#Plot an INVERSION *.LYR file from ParkSEIS

def plotLYRFile_INV(invFile, originalFile=None, title='Vs Profile', ax=None, vBounds=None, zBounds=None, legendFS=10, 
		save=False, fs=None, **kwargs):
	'''
	Takes at min one layer file which is output from inversion and plots it
	Optional second argument is the original model layer file which can be plotted for comparison
	'''

	#parse kwargs
	saveTitle = kwargs.get('saveTitle',None)

	#take file and turn into dict if only one is passed
	if isinstance(invFile, dict) == False:
		fileDict = {'':invFile}
	else:
		fileDict = invFile

	#load original file if passed
	if originalFile != None:
		origFileVals = readLYRFile(originalFile)
		z_orig = origFileVals['Depth']
		vs_orig = origFileVals['Vs']
		vp_orig = origFileVals['Vp']
		rho_orig = origFileVals['Density']

	#PLOT
	if fs != None: #set global font size for thesis plots
		plt.rcParams.update({'font.size': fs}) #global font size

	if ax != None:
		ax = ax
		msg = 'Master Plot'
	else:
		fig, ax = plt.subplots(1,1, figsize=(5,10))
		msg = 'No Master Plot'

	#plot original file first
	if originalFile != None:
		ax.plot(vs_orig, z_orig, c='k', lw=2, ls='--', label='Original')

	#plot all files in dict
	for fKey in fileDict:
		#load layer files
		invFileVals = readLYRFile(fileDict[fKey])
		z_inv = invFileVals['Depth']
		vs_inv = invFileVals['Vs']
		vp_inv = invFileVals['Vp']
		rho_inv = invFileVals['Density']


		majorYLocator = MultipleLocator(1) #sets the major tick interval to every meter
		majorYFormatter = FormatStrFormatter('%d')
		minorYLocator = MultipleLocator(0.25) #sets the minor tick interval to every 0.25 meter

		majorXLocator = MultipleLocator(100) #sets the major tick interval to every meter
		majorXFormatter = FormatStrFormatter('%d')
		minorXLocator = MultipleLocator(50) #sets the minor tick interval to every 0.25 meter

		ax.plot(vs_inv, z_inv, lw=2, label='Inverted: \n %s' %fKey)

	#format
	ax.set_xlabel('Vs [m/s]')
	ax.set_ylabel('Depth [m]')
	ax.invert_yaxis()
	#yticks
	ax.yaxis.set_major_locator(majorYLocator)
	ax.yaxis.set_major_formatter(majorYFormatter)
	ax.yaxis.set_minor_locator(minorYLocator)
	#xticks
	ax.xaxis.set_major_locator(majorXLocator)
	ax.xaxis.set_major_formatter(majorXFormatter)
	ax.xaxis.set_minor_locator(minorXLocator)
	
	ax.grid(which='both', linestyle='--', alpha=0.75)
	ax.set_axisbelow(True)
	ax.tick_params(axis='x', labelrotation=90)
	ax.legend(fontsize=legendFS)
	ax.set_title(title)

	#apply x and y limits if specified
	if vBounds != None:
		vB = np.asarray(vBounds)
		ax.set_xlim(vB[0], vB[1])
	if zBounds != None:
		zB = np.asarray(zBounds)
		ax.set_ylim(zB[0], zB[1])
		ax.invert_yaxis() #re invert after setting bounds

	
	if msg == 'No Master Plot':
		#fig.suptitle(title, y=0.92)
		plt.subplots_adjust(wspace=0.05, hspace=0)
		plt.show()

	if save:      
		plt.savefig(saveTitle, dpi=100, bbox_inches='tight')
		plt.show()


	return


##################
#Plot a 2D pseudo-inverse

def approxINV_2D(path, dx, zConv=0.3, vsConv=0.92, legend=True, returnDC=False, flip=False):
	'''
	Take a path to a folder containing experimental DC for a 2D profile and
	plots them as an approximate inversion

	Optionally returns all experimental dispersion curves in a sorted dictionary if desired

	path: path to extracted dispersion curves
	dx: spacing between each measurement
	zConv: wavelength conversion for depth
	vsConv: relation between measured phase velocity and shear velocity

	flip: option to flip plotting order (i.e. reverse profile) 

	**********NOTE************
	THIS ASSUMES SHOTS ARE NAMED IN SEQUENCE WITHIN THE FOLDER WITH
	THEIR RELATIVE POSITION AT THE END OF THE FILE NAME
	eg:
	Line01_Shot01(Model).DC
	Line01_Shot02(Model).DC
	...
	etc.
	'''

	##########
	#LOAD THE EXP. DC'S
	##########

	expDCs = {}

	for root, dirs, files in os.walk(path):
		for f in files:
			if f.endswith('.DC'):

				key = f.replace('.DC', '')

				expDCs[key] = loadTheorDC(os.path.join(path, f))

	#get reversed order dict keys for plotting to match DC data
	dcKeys_sorted = sorted(expDCs, reverse=True)
	
	##########
	#LOOP THROUGH AND PLOT
	##########

	#first calculate offset position vector
	xx = np.arange(0, dx*len(expDCs), dx)

	#format settings
	majorYLocator = MultipleLocator(1) #sets the major tick interval to every meter
	majorYFormatter = FormatStrFormatter('%d')
	minorYLocator = MultipleLocator(0.25) #sets the minor tick interval to every 0.25 meter

	majorXLocator = MultipleLocator(5) #sets the major tick interval to every meter
	majorXFormatter = FormatStrFormatter('%d')
	minorXLocator = MultipleLocator(1) #sets the minor tick interval to every 0.25 meter


	#set up and plot
	fig, ax = plt.subplots(1,2, figsize=(15,6), gridspec_kw={'width_ratios': [3, 1]})

	for (xPos, line, sortedKey) in zip(xx, expDCs, dcKeys_sorted):
		if flip:
			f, c = expDCs[sortedKey][0], expDCs[sortedKey][1]
		else:	
			f, c = expDCs[line][0], expDCs[line][1]

		#convert from frequency to wavelength
		lamb = c/f

		#calculate aprox depths and Vs values
		z = lamb*zConv
		vs = c/vsConv

		#create x positional vector for scatter 
		xcord = np.full_like(z, xPos)

		#set marker size base on spacing
		mSize = int(np.nanmax(xx)*4.5)
		ax[0].scatter(xcord, z, c=vs, cmap='plasma', s=mSize)
		ax[1].plot(vs, z, marker='o', fillstyle='none', label=line)

	####ax formatting
	def axFormat(ax):
		ax.set_ylabel('Depth [m]', fontsize=16)
		ax.tick_params(axis='both', labelsize=16)
		ax.invert_yaxis()
		#yticks
		ax.yaxis.set_major_locator(majorYLocator)
		ax.yaxis.set_major_formatter(majorYFormatter)
		ax.yaxis.set_minor_locator(minorYLocator)
		
		ax.set_axisbelow(True)
		ax.grid(which='both', linestyle='--', alpha=0.75)

	axFormat(ax=ax[0])
	axFormat(ax=ax[1])

	ax[0].set_xlabel('Offset [m]', fontsize=16)
	
	ax[1].set_xlabel('Vs [m/s]', fontsize=16)
	if legend:
		ax[1].legend()

	if returnDC:
		return expDCs
	else: 
		return


##################
#Plot a 2D velocity section created through inversion in ParkSEIS

def plotInverted2DSection(file, vLimits=None, zLimits=None, flip=False, ax=None, cBar='viridis'):
	'''
	Takes in a *.TXT file exported from ParkSEIS that contains a 2D inverted model
	and plots the section

	vLimits: list of bounding velocity values for plotting
	zLimits: list of depth limits for plotting

	flip: option to reverse profile direction
	'''
	
	#######
	#READ AND MANIPULATE DATA
	#######
	#read file to get number of data lines (rows)
	with open(file, 'r') as textFile:
		lines = textFile.readlines()
		dataLinesNum = len([l for l in lines if l.startswith(' ')]) #list comp
		
	#read file into arrays using genfromtext
	x, z, vs = np.genfromtxt(file, usecols=(0,1,2), max_rows=dataLinesNum, unpack=True)
	z = -1*z #get rid of negative sign
	
	#get some basic info
	uniqueXpos = np.unique(x).size #number of x locations (i.e. shots)
	numLay = dataLinesNum/uniqueXpos #number of layers in inverted model
	
	#create new arrays which also contain 'top layer' measurements for plotting
	epsilon = 0.001 #depth offset of top and bottom layer of adjacent rows
	
	xLay, zLay, vsLay = np.array([]), np.array([]), np.array([])
	
	for row in range(dataLinesNum): #loop through and append
		if row%numLay == 0: #for first data point at each new x location
			xLay, zLay, vsLay = np.append(xLay, x[row]), np.append(zLay, 0), np.append(vsLay, vs[row])
		else:
			xLay, zLay, vsLay = np.append(xLay, x[row]), np.append(zLay, z[row-1]+epsilon), np.append(vsLay, vs[row])
		
		xLay, zLay, vsLay = np.append(xLay, x[row]), np.append(zLay, z[row]), np.append(vsLay, vs[row])
		
	#######
	#PLOT
	#######
	
	if ax != None:
		ax = ax
		msg = 'Master Plot'
	else:
		fig, ax = plt.subplots(1,1,figsize=(15,8))
		msg = 'No Master Plot'
	
	#flip is desired
	if flip:
		xLay = xLay[::-1]

	if vLimits != None: #allow control of colorbar extents
		if np.nanmin(vsLay) < vLimits[0] and np.nanmax(vsLay) > vLimits[1]: ext = 'both'
		if np.nanmin(vsLay) < vLimits[0] and np.nanmax(vsLay) < vLimits[1]: ext = 'min'
		if np.nanmin(vsLay) > vLimits[0] and np.nanmax(vsLay) > vLimits[1]: ext = 'max'
		if np.nanmin(vsLay) > vLimits[0] and np.nanmax(vsLay) < vLimits[1]: ext = 'neither'

		#modify color bar accordingly
		levs = np.linspace(vLimits[0],vLimits[1], 50) #get levels for cbar
		ctr = ax.tricontourf(xLay, zLay, vsLay, levels=levs, extend=ext, cmap=cBar)
		ctr.set_clim(vLimits[0],vLimits[1])

	else:
		levs = np.linspace(min(vsLay), max(vsLay), 50)
		ctr = ax.tricontourf(xLay,zLay,vsLay, levels=levs, cmap=cBar)
	
	#plot original data location
	ax.scatter(xLay, zLay, c='k', marker='+', s=35, alpha = 0.55)
	
	#formatting
	if zLimits != None:
		ax.set_ylim(zLimits[0], zLimits[1])
	else:
		ax.set_ylim(min(zLay), max(zLay))
	ax.invert_yaxis()
	ax.set_ylabel('Depth [m]')
	
	ax.set_xlim(min(xLay), max(xLay))
	ax.set_xlabel('Position Along Line [m]')

	ax.set_title('Inverted Shear Velocity Profile')
	
	if msg == 'No Master Plot':
		cbar = fig.colorbar(ctr)
		cbar.set_label('Inverted Shear Velocity [m/s]')
		plt.show()

		return
	
	else:
		return ctr



##################
#Plot backscattered analysis section

def plotBSAsection(file, time=None, clim=None):
    '''
    Takes a BSA section from ParkSEIS that has been exported as a TXT file and plots it
    
    NOTE: this is a 'dumb' function as it will only work for the geometry I designed it for
    i.e. dSR of 2*dx, dx of 1.5m, and SR of 9 m with 24 geophones
    
    time: list of time range to plot
    clim: list containing min and max amplitudes to plot within (i.e. gain)
    '''
    
    #read data from the file
    x, t, amp = np.genfromtxt(file, usecols=(0,1,2), unpack=True)
    
    #get various parameters for calculations and loading
    numTr = np.unique(x).size #get the number of unique x positions or traces
    n = int(x.size / numTr) #number of sample in each trace
    dt = t[6] - t[5] #sample interval
    maxT = n * dt #max time
    
    #create empty array to store data
    bsaSection = np.full((numTr, n), np.nan)
    
    #loop through data and store
    trace = 0
    tSamp = 0
    for row in range(x.size): 
        bsaSection[trace, tSamp] = amp[row]
        tSamp += 1
        if tSamp%n == 0 and tSamp > 0:
            trace += 1
            tSamp = 0
            
    ##################################
    #PLOT
    ##################################
    fig, ax = plt.subplots(1,1, figsize=(15,8))
    
    majorXLocator = MultipleLocator(10) #sets the major tick interval to 10
    majorXFormatter = FormatStrFormatter('%d')
    minorXLocator = MultipleLocator(2) #sets the minor tick interval to every 5
    
    if time != None:
        tMin, tMax = int(time[0]/dt), int(time[1]/dt)
        eTop, eBttm = time[0], time[1]
    else:
        tMin, tMax = 0, (n-1)
        eTop, eBttm = 0, maxT
    ax.imshow(np.fliplr(bsaSection[:,tMin:tMax].T), aspect='auto', interpolation='bicubic', 
    	extent=[0-17.25,45+17.25+9,eBttm,eTop], clim=clim)
    ax.set_xlabel('Surface Location [m]', fontsize=20)
    ax.xaxis.tick_top()
    ax.xaxis.set_label_position('top') 
    ax.set_ylabel('Time [s]', fontsize=20)
    ax.tick_params(axis='both', labelsize=16)
    ax.grid(which='both', linestyle='--', alpha=0.75)
    
    ax.xaxis.set_major_locator(majorXLocator)
    ax.xaxis.set_major_formatter(majorXFormatter)
    ax.xaxis.set_minor_locator(minorXLocator)

    return


##################
#SUMMARY PLOTS
##################

##################
#Plot a summary image of layer model, dispersion image, and dispersion curves

def plotModelAnalysis(modPath, fmin=2, fmax=100, vmin=50, vmax=1500, dv=0.1, dx=1.5, X1=9, printParams=False,
	fs=12, save=False):
	'''
	Takes a path to a model folder where results are stored in folders of:
		'DC' - all dispersion curves and associated .LYR file
		'seis' - all seismic model files generated (i.e. .dat files)
	'''
	
	#get DC folders, seismic paths, and LYR paths first
	dcFolder, seisFiles, lyrFiles = getModelFolderPaths(modPath)
	
	###
	#MAIN LOOP THROUGH MODEL FILES
	###
	
	for dcF, seisF, lyrF in zip(dcFolder, seisFiles, lyrFiles):

		if printParams: #print out the model parameters
			printLYRfile(lyrF)

		#LOAD SEISMIC
		modTrcs, modStrm = loadModelDAT(seisF)
		
		#LOAD IN EXPERIMENTAL DC FILES
		dcDict = loadTheorDC_PATH(dcF)
		
		###
		#PLOT
		###
		plt.rcParams.update({'font.size': fs}) #global font size

		fig = plt.figure(figsize=(15,5))
		gs = gridspec.GridSpec(1,5, width_ratios=[0.4,0.4,0.4,0.3,3.25], wspace=0.05)
		
		#get plot title
		plotTitle = seisF.split('\\')[-1] #gets title for model
		
		#LAYER PLOT
		lyrAX = [plt.subplot(gs[0,0]), plt.subplot(gs[0,1]), plt.subplot(gs[0,2])] #create list of axes to pass to function
		plotLYRFile(lyrF, ax=lyrAX)
		lyrAX[0].set_xlim(50,350)
		lyrAX[0].locator_params(axis='x', tight=True, nbins=8)
		lyrAX[0].tick_params(axis='x', which='major', labelsize=fs*0.75)
		lyrAX[1].tick_params(axis='x', which='major', labelsize=fs*0.75)
		lyrAX[2].set_xlabel('$\\rho$ [kg/m$^3$]')
		lyrAX[2].tick_params(axis='x', which='major', labelsize=fs*0.75)
		
		#BLANK SPACER PLOT
		axSPACE = plt.subplot(gs[0,3])
		axSPACE.axis('off')

		#DISP IMAGE PLOT
		dcAX = plt.subplot(gs[0,4])
		dcModelPhaseShift(modStrm, minVel=vmin, maxVel=vmax, dv=dv, minFreq=fmin, maxFreq=fmax, padLen=500,
			X1=X1, dx=dx, overLayDC=dcDict, velSum=True, ax=dcAX)
		dcAX.set_title('')
		dcAX.get_legend().remove() #removes legend; over crowded and can't see response
				
		#Set title and save
		if save == False:
			supT = fig.suptitle(plotTitle, x=0.445, y=0.99, fontsize=14)

		if save:
			 figTitle = plotTitle.replace('.dat','')       
			 plt.savefig(figTitle, dpi=100, bbox_inches='tight')      
	   
	
		plt.show()
		
	return


##################
#Plot an inversion summary image

def plotINVSummary(seisFile, expDC, invLYRFile, invDC, oneD=True, origLYRFile=None, fmin=2, fmax=100, 
	vmin=50, vmax=1500, dv=0.1, dx=1.5, X1=9, vBounds=None, zBounds=None, save=False, fs=None):
	'''
	Takes at minimum a original seismic file, extracted experimental DC, and a inverted model contained in a .LYR file and plots
	a summary figure.  

	Optional input is a original .LYR for comparison which can either be the known true model or initial model for inversion

	INPUTS:
	seisFile: .dat file (vertical comp.) containing seismic record that was inverted
	expDC: the extracted fundamental mode file from the DC image generated from the provided seisFile
	invDC: file pointing to the theoretical DC calculated from the inverted model
	invLYRFile: .LYR file containing the resultant layered earth model generated during inversion
	origLYRFile: the original or initial .LYR file containing the layered earth model
	'''
	
	#load traces and stream from provided seismic file
	t, s = loadModelDAT(seisFile)

	#load the extracted experimental DC curve provided into dict
	DC_exp = loadTheorDC(expDC, asDict=True, dictKey='Experimental DC')

	#if provided inverted DC is single file, load the modeled DC curve provided into dict
	if isinstance(invDC, dict) == False:
		DC_inv = loadTheorDC(invDC, asDict=True, dictKey='Inverted Model DC')
	elif isinstance(invDC, dict):
		DC_inv = invDC #store provided dict into var to be merged

	#merge two DC into single dict for plotting
	DC_plotDict = {**DC_exp, **DC_inv}
	
	###
	#PLOT SUMMARY FIG
	###

	if fs != None: #set global font size 
		plt.rcParams.update({'font.size': fs}) #global font size
	
		
	fig = plt.figure(figsize=(14,5.5))
	gs = gridspec.GridSpec(1,2, width_ratios=[0.2,0.9])
		
	#get plot title
	plotTitle = 'Inversion Summary for ' + seisFile.split('\\')[-1] #gets title for model
		
	#LAYER PLOT
	lyrAX = plt.subplot(gs[0,0])
	plotLYRFile_INV(invFile=invLYRFile, originalFile=origLYRFile, vBounds=vBounds, zBounds=zBounds, ax=lyrAX)
	if fs != None:
		lyrAX.legend(fontsize=fs*.75)
		lyrAX.set_title('')

		
	#DISP IMAGE PLOT
	dcAX = plt.subplot(gs[0,1])
	dcModelPhaseShift(s, minVel=vmin, maxVel=vmax, minFreq=fmin, dv=dv, maxFreq=fmax, padLen=1000,
		X1=X1, dx=dx, overLayDC=DC_plotDict, ax=dcAX)
	if fs != None:
		dcAX.legend(fontsize=fs*.75)
		dcAX.set_title('')
	if save == False:	
		fig.colorbar(dcAX.images[0], ax=dcAX, shrink=0.7)
				
	#Set title and save

	if save == False: supT = fig.suptitle(plotTitle, x=0.445, y=0.97, fontsize=14)

	if save:
		figTitle = 'INVsummary_' + seisFile.split('\\')[-1]
		figTitle = figTitle.replace('.dat','')        
		plt.savefig(figTitle, dpi=100, bbox_inches='tight')      
	   
	
	plt.show()

	#PRINT INV SUMMARY
	if isinstance(invLYRFile, dict) == False:
		oneORtwo = oneD #if the inverted layer files are from a 1D or 2D inversion
		printLYRfile(invLYRFile, inv=oneORtwo)
	
	return


##################
#Plot a 2D velocity profile inversion summary

def inversionSUMMARY_2D(file, vLimits=None, zLimits=None, cBar='viridis', legend=True, flip=False, returnMODLyrs=False):
	'''
	Takes in a .TXT file from ParkSEIS that contains 2D inverted model and plots
	the 2D section and all the individual layer file inversions which make it up

	Optionally returns dictionary containing all layer files associated with the inversion

	vLimits: list of bounding velocity values for plotting
	zLimits: list of depth limits for plotting

	flip: option to reverse profile direction
	'''

	########
	#First get individual layer files
	########
	XRLmodFile = file.replace('(2DVs)','')
	XRLmodFile = XRLmodFile.replace('.TXT', '(Model).XRL')

	#read the file
	with open(XRLmodFile, 'r') as XRLfile:
		lines = XRLfile.readlines()
		numLayerFiles = len(lines)

	#load layer file paths into array
	modLYRfiles = np.genfromtxt(XRLmodFile, usecols=(2), max_rows=numLayerFiles, unpack=True, dtype='str') #load file paths
	LYRFilesXpos = np.genfromtxt(XRLmodFile, usecols=(0), max_rows=numLayerFiles, unpack=True, dtype='str') #load shot numbers
	LYRFilesXpos = [float(xPos.replace(',', '')) for xPos in LYRFilesXpos] #convert shot numbers to numbers

	modelLYRdict = dict(zip(LYRFilesXpos, modLYRfiles)) #store layer file paths and associated shot nums in dict

	#########
	#PLOT
	#########

	fig, ax = plt.subplots(1,2, figsize=(20,8), gridspec_kw={'width_ratios': [3, 1], 'wspace': 0.29})

	#plot 2D section
	ctr = plotInverted2DSection(file, vLimits=vLimits, zLimits=zLimits, cBar=cBar, flip=flip, ax=ax[0])
	cbar = fig.colorbar(ctr, ax=ax[0], shrink=0.9, fraction=0.05, aspect=30, format='%d')
	cbar.set_label('Inverted Shear Velocity [m/s]', fontsize=16)
	cbar.ax.tick_params(labelsize=14)

	#plot layer files
	plotLYRFile_INV(invFile=modelLYRdict, ax=ax[1])
	ax[1].set_ylim(ax[0].get_ylim()) #sets y limits of plots to be the same for alllll the asssthetics
	if legend == False:
		ax[1].get_legend().remove()

	for item in ([ax[0].title, ax[0].xaxis.label, ax[0].yaxis.label] + ax[0].get_xticklabels() + ax[0].get_yticklabels()):
		item.set_fontsize(20)
	for item in ([ax[1].title, ax[1].xaxis.label, ax[1].yaxis.label] + ax[1].get_xticklabels() + ax[1].get_yticklabels()):
		item.set_fontsize(20)

	plt.tight_layout()

	if returnMODLyrs:
		return modelLYRdict
	else:
		return

##################
#UTILITY CODE
##################

##################
#Determine properties of dispersion curve file

def getTheorDCParams(fname):
	'''
	Takes in a file path to a ParkSeis generate dispersion curve and extracts the number of points and mode number
	'''
	with open(fname, 'r') as dispFile:
		lines = dispFile.readlines()

		#get number of points on DC 
		numPts = int(lines[0][7:]) #number of points begins after the 7th character

		#get the mode number
		if fname.endswith('(Model).DC'):
			modeNum = 0
		else:
			modeNum = int(fname[-5:-4])	#get mode number from provided file name

	return numPts, modeNum


##################
#Load all relevant folder/file paths for multiple models at once

def getModelFolderPaths(path):

	'''
	Loop through a specified path and return modeled DC folder paths, paths to modeled seismic files, and paths to model .LYR files
	
	NOTE: 'path' specified must contain a folder called 'DC' holding dispersion curves and a LYR file, as well as a 'seis' folder
	containing all modeled seismic files
	'''

	dcModelFolders = []
	seisModelFolders = []
	seisFiles = []
	lyrFiles = []
	
	#get folder paths fopr DC curves and seismic files
	for root, dirs, files in os.walk(path):
		if root.endswith('DC'):
			#print('Dispersion curves found in folder: ', root)
			dcModelFolders.append(root)
		if root.endswith('seis'):
			#print('Seismic folder found in folder: ', root)
			seisModelFolders.append(root)
			
	#get seismic file paths
	for p in seisModelFolders:
		for entry in os.listdir(p):
			if entry.endswith('v.dat'): #SELECT ONLY VERTICAL COMP OUTPUT
				seisFiles.append(os.path.join(p,entry))

	#get layer file paths
	for p in dcModelFolders:
		for entry in os.listdir(p):
			if entry.endswith('.LYR'):
				lyrFiles.append(os.path.join(p,entry))
				
	if len(dcModelFolders) != len(seisFiles):
		raise Exception('Caution! The number of model dispersion curve folders is not the same as the number of model seismic files!')
	if len(dcModelFolders) != len(lyrFiles):
		raise Exception('Caution! The number of model dispersion curve folders is not the same as the number of model .LYR files!')
				
	return dcModelFolders, seisFiles, lyrFiles


##################
#Print out parameters of a ParkSEIS *.LYR file

def printLYRfile(fname, inv=False):
	'''
	Takes in a layer file from ParkSEIS and prints out the table of model parameters
	'''
	with open(fname, 'r') as lyrFile:
		lines = lyrFile.readlines()
		numLayers = int(lines[2][14:]) #number of layers is on 3rd line starting after the 14 col
		paramLines = lines[1:(4+numLayers)] #gets all lines containing model parameters

		#for printing inversion file parameters
		if inv:
			#get inv parameters
			sectionLineNums = []
			for i, l in enumerate(lines):
				if l.startswith('------------------------------------------------------------'):
					sectionLineNums.append(i)
	
			if sectionLineNums: #if no section lines exists, still print results
				invSecStart, invSecEnd = sectionLineNums[0]+1, sectionLineNums[1]-2 #gets the line range for inversion parameters
				matchSecStart, matchSecEnd = sectionLineNums[1]+1, sectionLineNums[2] #same but for match info
	
				invParams = lines[invSecStart:invSecEnd]
				matchParams = lines[matchSecStart:matchSecEnd]

		
	print('---- PARAMETERS FOR MODEL %s: \n' %fname)
	for l in paramLines:
		print(l)

	if inv and (sectionLineNums): #again, make sure sections actually exist so can still print without error
		print('----INVERSION PARAMETERS')
		for l in invParams:
			print(l)
		print('----INVERSION MATCH')
		for l in matchParams:
			print(l)

	return


#################################################
#################################################
###RESISTIVITY CODE
#################################################

##################
#FILE LOADING CODE
##################

##################
#Load Raw Field Data

def loadProSysData(fname):
	''' 
	Takes in a file path to a *.csv file exported from ProSys
	Returns:
	Tx1: Transmitter 1 location
	Tx2: Transmitter 2 location
	Rx1: Rec. location 1
	Rx2: Rec. location 2
	a: base a spacing
	n: all n values for observations
	xLoc: lateral location for each measurement 
	rho: apparent resistivity for each measurement 
	'''

	Tx1, Tx2, Rx1, Rx2, rho = np.genfromtxt(fname, delimiter=',', skip_header=1,
		usecols=(2,3,4,5,6), unpack=True)

	#find base 'a' spacing
	aSpacing = Tx2[0] - Tx1[0]

	#calculate lateral location of each measurement
	xLoc = ((Rx1 - Tx2) / 2) + Tx2

	#calculate 'n' value for each measurement
	n = np.zeros_like(Tx1)

	for obs in range(Tx1.size):
		n[obs] = ((Rx1[obs] - Tx1[obs]) / aSpacing) - 1

	return Tx1, Tx2, Rx1, Rx2, rho, xLoc, n, aSpacing 



##################
#Load Modeled *.dat Pseudosection Data

def loadRES2DModData(fname):
	'''
	Takes in a *.dat file created from resd2dmod and returns
	xLoc: lateral position of measurement
	n: n value for measurement 
	rho: apparent resistivity for measurement
	aSpac: unit electrode spacing
	'''

	#get the number of inversion model blocks
	with open(fname, 'r') as datModFile:
		lines = datModFile.readlines()
		numObs = int(lines[3])	#info we want is on line 4 and is from location 20 to the end
		aSpac = int(lines[1][3]) #a spacing on line 2

	xLoc, n, rho = np.genfromtxt(fname, skip_header=6, usecols=(0,2,3), max_rows=numObs, unpack=True)

	return xLoc, n, rho, aSpac


##################
#Load inverted Data

def loadRes2DINVResults(fname):
	''' 
	Takes in a file path to a *.xyz file from res2dinv
	Returns:
	x: central location of model blocks
	z: this is the depth value for the center of the blocks
	rho: this is the inverted resistivity
	'''

	#get the number of inversion model blocks
	with open(fname, 'r') as xyzINVFile:
		lines = xyzINVFile.readlines()
		invBlockNum = int(lines[1][21:])  #info we want is on line 2 and is from location 20 to the end

	#extract info from file
	x, z, rho = np.genfromtxt(fname, skip_header=5, usecols=(0,1,2), max_rows=invBlockNum, unpack=True)

	return x, z, rho


##################
#PLOTTING CODE
##################

##################
#Plot Raw Field Data

def plotRawData(fname, sliceLoc, nSliceVal, plotTitle, gridDimension=500, limits=None, 
	showPoints=False, save=False, fs=12):
	'''
	Takes in a file path to a *.csv (or modeled *.dat file) file exported from ProSys
	Plots raw data as a pseudosection

	fname: file to be plotted
	sliceLoc: lateral location along profile to plot 1D sounding curve
	nSliceVal: 'n' value to plot a slice at
	gridDimension: the dimension of square grid to interpolate data onto 
	limits: set resistivity limits to contour within
	showPoints: show locations of true data points on plotted section
	fs: font size used during plotting
	'''

	#####################
	#Load the file
	if fname.endswith('.csv'):
		Tx1, Tx2, Rx1, Rx2, rho, xLoc, n, aSpacing = loadProSysData(fname)
	elif fname.lower().endswith('.dat'):
		xLoc, n, rho, aSpacing = loadRES2DModData(fname)


	#####################
	#Grid the Data

	#get max values
	max_x0, max_y0, max_z0 = np.amax(xLoc), np.amax(n), np.amax(rho)

	#get lateral step 
	dx = (max_x0 + aSpacing + 1) / gridDimension

	#establish grid
	xi = np.linspace(0, max_x0 + aSpacing + 1, gridDimension)
	yi = np.linspace(0, max_y0, gridDimension)

	#finally, grid the data
	zi = griddata((xLoc, n), rho, (xi[None,:], yi[:,None]), method='cubic')


	#####################
	#Plot the gridded data

	#get depth slice
	numNLevels = int(np.amax(yi)) #number of n levels
	dy = numNLevels / yi.size #y step of grid
	nSliceLoc = int(nSliceVal/dy) #array index for desired n level
	nSlice = zi[nSliceLoc,:] #get values at index

	#get 1D profile
	oneDProf = zi[:, int(sliceLoc/dx)] #get values at xloc

	#init figure
	fig = plt.figure(figsize=(17,6))
	gs = gridspec.GridSpec(2,3, width_ratios=[1,5,0.15], height_ratios=[3,1], hspace=.42, wspace=0.28)

	#main pseudosection
	ax0 = plt.subplot(gs[0,1])
	if limits != None: #allow control of colorbar extents
		if np.nanmin(zi) < limits[0] and np.nanmax(zi) > limits[1]: ext = 'both'
		if np.nanmin(zi) < limits[0] and np.nanmax(zi) < limits[1]: ext = 'min'
		if np.nanmin(zi) > limits[0] and np.nanmax(zi) > limits[1]: ext = 'max'
		if np.nanmin(zi) > limits[0] and np.nanmax(zi) < limits[1]: ext = 'neither'

		#modify color bar accordingly
		ctr = ax0.contourf(xi,yi,zi,np.arange(limits[0],limits[1],0.5), extend=ext)
		ctr.set_clim(limits[0],limits[1])
	else:
		ctr = ax0.contourf(xi,yi,zi)
	#draw contours
	ax0.contour(xi,yi,zi, 15, colors='gray', linewidths=0.75)

	if showPoints == True:
		ax0.scatter(xLoc, n, c='k', marker='+', s=35, alpha = 0.55)

	#main plot formatting 
	ax0.invert_yaxis()
	ax0.axvline(x=sliceLoc)
	ax0.axhline(y=nSliceVal)
	ax0.set_ylabel('n', fontsize=fs)
	ax0.set_title('Apparent Resistivity Pseudosection', fontsize=fs)
	ax0.tick_params(axis='both', which='major', labelsize=fs*0.75)
	ax0.tick_params(axis='both', which='minor', labelsize=fs*0.75)

	#plot cBar
	ax3 = plt.subplot(gs[:,2])
	cbar = fig.colorbar(ctr, cax=ax3)
	cbar.set_label('Apparent Resistivity [$\Omega .m$]', fontsize=fs)
	ax3.tick_params(axis='both', which='major', labelsize=fs*0.75)
	ax3.tick_params(axis='both', which='minor', labelsize=fs*0.75)

	#plot 1D profile at specified x location
	ax1 = plt.subplot(gs[:,0])
	ax1.plot(oneDProf, yi)
	ax1.scatter(rho[xLoc == sliceLoc], n[xLoc == sliceLoc])
	ax1.invert_yaxis()
	ax1.grid(which='both')
	ax1.set_ylabel('n', fontsize=fs)
	ax1.set_xlabel('$\\rho _a$ [$\Omega .m$]', fontsize=fs)
	ax1.set_title('$\\rho _a$ at %d m' %sliceLoc, fontsize=fs)
	ax1.tick_params(axis='both', which='major', labelsize=fs*0.75)
	ax1.tick_params(axis='both', which='minor', labelsize=fs*0.75)

	#plot n slice at specified level
	ax2 = plt.subplot(gs[1,1])
	ax2.plot(xi, nSlice)
	ax2.scatter(xLoc[n == nSliceVal], rho[n == nSliceVal])
	ax2.set_ylim(np.nanmin(nSlice)*0.8, np.nanmax(nSlice)*1.2)
	ax2.set_xlim(np.nanmin(xi), np.nanmax(xi))
	ax2.grid(which='both')
	ax2.set_ylabel('$\\rho _a$ [$\Omega .m$]', fontsize=fs)
	ax2.set_xlabel('Profile Distance [m]', fontsize=fs)
	ax2.set_title('$\\rho _a$ at n = %d' %nSliceVal, fontsize=fs)
	ax2.tick_params(axis='both', which='major', labelsize=fs*0.75)
	ax2.tick_params(axis='both', which='minor', labelsize=fs*0.75)

	if save ==False:
		supT = fig.suptitle(plotTitle, x=0.445, y=1.04, fontsize=15)

	figTitle = plotTitle.replace('.dat','')
	figTitle = plotTitle.replace('.csv','')

	if save:
		plt.savefig(figTitle, dpi=100, bbox_inches='tight')

	return


##################
#Plot Inverted Data

def plotINVResults(fname, sliceLoc, depthSlice, plotTitle, limits=None, AOI=None, 
	bounds=None, bCols=None, save=False, fs=12):
	'''
	Takes raw output from res2dinv and plots using tricontour (Delaunay triangulation)
	NOTE: this really only works for inversions with a significant number of model blocks and 
	extended models.  Use interpolation for more sparse data

	limits: specify the min and max resistivity value to allow in the colorbar
	AOI: this controls the plot extent and allows 'zooming' into an area of interest (x1,x2,z1,z2)
	bounds: list of n values between which resistivity values will be grouped
	bCols: list of n-1 colors which will represent the bounded values specified in bounds
	'''

	#####################
	#Load the file
	x, z, rho = loadRes2DINVResults(fname)


	#####################
	#Plot the inverted data

	fig = plt.figure(figsize=(17,6))

	gs = gridspec.GridSpec(2,3, width_ratios=[1,5,0.15], height_ratios=[3,1], hspace=.42,wspace=0.28)

	#main plot
	ax0 = plt.subplot(gs[0,1])

	#plot limited countours if specified
	if limits != None: #allow control of colorbar extents
		if np.nanmin(rho) < limits[0] and np.nanmax(rho) > limits[1]: ext = 'both'
		if np.nanmin(rho) < limits[0] and np.nanmax(rho) < limits[1]: ext = 'min'
		if np.nanmin(rho) > limits[0] and np.nanmax(rho) > limits[1]: ext = 'max'
		if np.nanmin(rho) > limits[0] and np.nanmax(rho) < limits[1]: ext = 'neither'

		#modify color bar accordingly
		levs = np.linspace(limits[0],limits[1], 20) #get levels for cbar
		ctr = ax0.tricontourf(x, z, rho, levels=levs, extend=ext)
		ctr.set_clim(limits[0],limits[1])
	
	#bounded resistivity range plot
	elif (bounds != None) and (bCols != None):
		cmap = matplotlib.colors.ListedColormap(bCols)
		norm = matplotlib.colors.BoundaryNorm(bounds, len(bCols), clip=True)

		levs = bounds #set levels for cbar
		ctr = ax0.tricontourf(x, z, rho, norm=norm, cmap=cmap, levels=levs, alpha=0.9)

	#plot regular contours if not specified
	else:
		levs = 10 #set levels for cbar
		ctr = ax0.tricontourf(x, z, rho, levels=levs)

	#draw contours
	ax0.tricontour(x, z, rho, colors='gray', linewidths=0.75, levels=levs)


	#if area of interest if specified
	if AOI:
		aoi_x0, aoi_x1, aoi_z0, aoi_z1 = AOI[0], AOI[1], AOI[2], AOI[3]

	if AOI:
		ax0.set_xlim(aoi_x0, aoi_x1)
		ax0.set_ylim(aoi_z0, aoi_z1)
	ax0.invert_yaxis()
	ax0.axvline(x=sliceLoc)
	ax0.axhline(y=depthSlice)
	ax0.set_ylabel('Depth [m]', fontsize=fs)
	ax0.set_title('Inverted Resistivity Profile', fontsize=fs)
	ax0.tick_params(axis='both', which='major', labelsize=fs*0.75)
	ax0.tick_params(axis='both', which='minor', labelsize=fs*0.75)

	#color bar for main plot
	ax3 = plt.subplot(gs[:,2])
	cbar = fig.colorbar(ctr, cax=ax3)
	cbar.set_label('Resistivity [$\Omega .m$]', fontsize=fs)
	ax3.tick_params(axis='both', which='major', labelsize=fs*0.75)
	ax3.tick_params(axis='both', which='minor', labelsize=fs*0.75)

	#depth profile
	ax1 = plt.subplot(gs[:,0])
	profileX, profileY = rho[x == find_nearest(x, sliceLoc)], np.unique(z)
	ax1.scatter(profileX, profileY)
	ax1.plot(profileX, profileY)
	if AOI:
		ax1.set_ylim(aoi_z0, aoi_z1)
	ax1.invert_yaxis()
	ax1.grid(which='both')
	ax1.set_ylabel('Depth [m]', fontsize=fs)
	ax1.set_xlabel('$\\rho$ [$\Omega .m$]', fontsize=fs)
	ax1.set_title('$\\rho$ at %d m' %sliceLoc, fontsize=fs)
	if limits != None:
		ax1.set_xlim(limits[0], limits[1])
	ax1.tick_params(axis='both', which='major', labelsize=fs*0.75)
	ax1.tick_params(axis='both', which='minor', labelsize=fs*0.75)

	#depth slice
	ax2 = plt.subplot(gs[1,1])
	depthX, depthY = np.unique(x), rho[z == find_nearest(z, depthSlice)]
	ax2.scatter(depthX, depthY)
	ax2.plot(depthX, depthY)
	if AOI:
		ax2.set_xlim(aoi_x0, aoi_x1)
	else:
		ax2.set_xlim([np.amin(x), np.amax(x)])
	ax2.grid(which='both')
	ax2.set_ylabel('$\\rho$ [$\Omega .m$]', fontsize=fs)
	ax2.set_xlabel('Profile Distance [m]', fontsize=fs)
	ax2.set_title('Inverted Resistivity at z = %d' %depthSlice, fontsize=fs)
	ax2.tick_params(axis='both', which='major', labelsize=fs*0.75)
	ax2.tick_params(axis='both', which='minor', labelsize=fs*0.75)
	

	if save ==False:
		supT = fig.suptitle(plotTitle, x=0.445, y=1.04, fontsize=15)

	figTitle = plotTitle.replace('.dat','')
	figTitle = plotTitle.replace('.xyz','')

	if save:
		plt.savefig(figTitle, dpi=100, bbox_inches='tight')

	return


##################
#Plot Percent Change in Raw Data

def plotRAWDataPercentChange(baseFile, secondFile, sliceLoc, nSliceVal, plotTitle, gridDimension=500, 
	showPoints=False, method='absolute', save=False, fs=12):
	'''
	Takes raw output from either 2 proSys *.csv files or 2 *.dat files, 
	calculates the difference, and plots the percent change from the base survey

	baseFile: this is the initial or base survey
	secondFile: this is the second or follow up survey
	method: 'absolute' (default) or 'signed'
	'''

	#####################
	#Load base survey file
	if baseFile.endswith('.csv'):
		Tx1_B, Tx2_B, Rx1_B, Rx2_B, rho_B, xLoc_B, n_B, aSpacing_B = loadProSysData(baseFile)
	elif baseFile.endswith('.dat'):
		xLoc_B, n_B, rho_B, aSpacing_B = loadRES2DModData(baseFile)

	#####################
	#Load second survey file
	if secondFile.endswith('.csv'):
		Tx1_2, Tx2_2, Rx1_2, Rx2_2, rho_2, xLoc_2, n_2, aSpacing_2 = loadProSysData(secondFile)
	elif secondFile.endswith('.dat'):
		xLoc_2, n_2, rho_2, aSpacing_2 = loadRES2DModData(secondFile)


	#####################
	#Calculate difference between two surveys and set color map
	if method == 'absolute':
		diffAbs = np.absolute(rho_B - rho_2)
		percentChange = np.absolute((diffAbs/rho_B)*100)
		cMap = 'viridis'
	if method == 'signed':
		diff_signed = rho_B - rho_2
		percentChange = ((diff_signed/rho_B)*100) * (-1.0) #negative one so that positive percents are increases
		cMap = 'seismic'


	#####################
	#Grid the Data

	#get max values
	max_x0, max_y0, max_z0 = np.amax(xLoc_B), np.amax(n_B), np.amax(percentChange)

	#get lateral step 
	dx = (max_x0 + aSpacing_B + 1) / gridDimension

	#establish grid
	xi = np.linspace(0, max_x0 + aSpacing_B + 1, gridDimension)
	yi = np.linspace(0, max_y0, gridDimension)

	#finally, grid the data
	zi = griddata((xLoc_B, n_B), percentChange, (xi[None,:], yi[:,None]), method='cubic')
	#get min and max values for plotting color bars etc.
	if method == 'absolute':
		minVal, maxVal = 0.0, np.nanmax(zi)
		zi[(zi != np.nan) & (zi < 0)] = 0 #set interpolated values below 0 to 0
	if method == 'signed':
		pChangeMax = np.nanmax(np.absolute(zi)) #get max percent change value for setting cbar limits (i.e centered)
		minVal, maxVal = (-1.0*pChangeMax), pChangeMax


	#####################
	#Plot the gridded data

	#get depth slice
	numNLevels = int(np.amax(yi)) #number of n levels
	dy = numNLevels / yi.size #y step of grid
	nSliceLoc = int(nSliceVal/dy) #array index for desired n level
	nSlice = zi[nSliceLoc,:] #get values at index

	#get 1D profile
	oneDProf = zi[:, int(sliceLoc/dx)] #get values at xloc

	#init figure
	fig = plt.figure(figsize=(17,6))
	gs = gridspec.GridSpec(2,3, width_ratios=[1,5,0.15], height_ratios=[3,1], hspace=.42, wspace=0.28)

	#main pseudosection
	ax0 = plt.subplot(gs[0,1])
	ctr = ax0.contourf(xi,yi,zi, cmap=cMap, levels=np.linspace(minVal, maxVal, 50))
	
	#draw contours
	ax0.contour(xi,yi,zi, 15, colors='gray', linewidths=0.75)

	if showPoints == True:
		ax0.scatter(xLoc_B, n_B, c='k', marker='+', s=35, alpha = 0.55)

	#main plot formatting 
	ax0.invert_yaxis()
	ax0.axvline(x=sliceLoc)
	ax0.axhline(y=nSliceVal)
	ax0.set_ylabel('n', fontsize=fs)
	ax0.set_title('Percent Change in Apparent Resistivity', fontsize=fs)
	ax0.tick_params(axis='both', which='major', labelsize=fs*0.75)
	ax0.tick_params(axis='both', which='minor', labelsize=fs*0.75)

	#plot cBar
	ax3 = plt.subplot(gs[:,2])
	cbar = fig.colorbar(ctr, cax=ax3, ticks=np.linspace(minVal, maxVal, 11), format='%.1f') #11 labels on color bar
	cbar.set_label('% Change in Apparent Resistivity', fontsize=fs)
	ax3.tick_params(axis='both', which='major', labelsize=fs*0.75)
	ax3.tick_params(axis='both', which='minor', labelsize=fs*0.75)

	#plot 1D profile at specified x location
	ax1 = plt.subplot(gs[:,0])
	ax1.plot(oneDProf, yi)
	ax1.scatter(percentChange[xLoc_B == sliceLoc], n_B[xLoc_B == sliceLoc])
	ax1.invert_yaxis()
	ax1.grid(which='both')
	ax1.set_ylabel('n', fontsize=fs)
	ax1.set_xlabel('% Change', fontsize=fs)
	ax1.set_title('%% Change at %d' %sliceLoc, fontsize=fs)
	ax1.tick_params(axis='both', which='major', labelsize=fs*0.75)
	ax1.tick_params(axis='both', which='minor', labelsize=fs*0.75)

	#plot n slice at specified level
	ax2 = plt.subplot(gs[1,1])
	ax2.plot(xi, nSlice)
	ax2.scatter(xLoc_B[n_B == nSliceVal], percentChange[n_B == nSliceVal])
	ax2.set_ylim(np.nanmin(nSlice)*0.8, np.nanmax(nSlice)*1.2)
	ax2.set_xlim(np.nanmin(xi), np.nanmax(xi))
	ax2.grid(which='both')
	ax2.set_ylabel('% Change', fontsize=fs)
	ax2.set_xlabel('Profile Distance [m]', fontsize=fs)
	ax2.set_title('%% Change at n = %d' %nSliceVal, fontsize=fs)
	ax2.tick_params(axis='both', which='major', labelsize=fs*0.75)
	ax2.tick_params(axis='both', which='minor', labelsize=fs*0.75)

	if save == False:
		supT = fig.suptitle(plotTitle, x=0.445, y=1.04, fontsize=15)

	figTitle = plotTitle.replace('.dat','')
	figTitle = plotTitle.replace('.csv','')

	if save:
		plt.savefig(figTitle, dpi=150, bbox_inches='tight')

	return


##################
#Plot Percent Change in Inverted Data

def plotINVPercentChange(baseFile, secondFile, sliceLoc, depthSlice, plotTitle, method='absolute',
	save=False, fs=12):
	'''
	Takes raw output from 2 resd2dinv *.xyz files, calculated the difference, and plots the percent
	change from the base survey

	NOTE: can also pass a float value as second file to subtract out a background value

	baseFile: this is the initial or base survey
	secondFile: this is the second or follow up survey
	method: 'absolute' (default) or 'signed'
	'''

	#####################
	#Load base survey file
	x_B, z_B, rho_B = loadRes2DINVResults(baseFile)

	#test if second file is background value or *.xyz file
	if isinstance(secondFile, str):
		#####################
		#Load second survey file
		x_2, z_2, rho_2 = loadRes2DINVResults(secondFile)
	else:
		x_2, z_2 = x_B, z_B
		rho_2 = np.full_like(rho_B, secondFile) #creates array like the base file and fills with specified value


	#####################
	#Calculate difference between two surveys
	if method == 'absolute':
		diffAbs = np.absolute(rho_B - rho_2)
		percentChange = np.absolute((diffAbs/rho_B)*100)
		minVal, maxVal = 0.0, np.nanmax(percentChange)
		cMap = 'viridis'
	if method == 'signed':
		diff_signed = rho_B - rho_2
		percentChange = ((diff_signed/rho_B)*100) * (-1.0) #negative one so that positive percents are increases
		pChangeMax = np.nanmax(np.absolute(percentChange)) #get max percent change value for setting cbar limits (i.e centered)
		minVal, maxVal = (-1.0*pChangeMax), pChangeMax
		cMap = 'seismic'


	#####################
	#Plot the percent change

	fig = plt.figure(figsize=(17,6))

	gs = gridspec.GridSpec(2,3, width_ratios=[1,5,0.15], height_ratios=[3,1], hspace=.42,wspace=0.28)

	#main plot
	ax0 = plt.subplot(gs[0,1])
	ctr = ax0.tricontourf(x_2, z_2, percentChange, levels=np.linspace(minVal, maxVal, 50), cmap=cMap) #levels to control 0 point
	ax0.invert_yaxis()
	ax0.axvline(x=sliceLoc)
	ax0.axhline(y=depthSlice)
	ax0.set_ylabel('Depth [m]', fontsize=fs)
	ax0.set_title('Percent Change in Inverted Resistivity', fontsize=fs)
	ax0.tick_params(axis='both', which='major', labelsize=fs*0.75)
	ax0.tick_params(axis='both', which='minor', labelsize=fs*0.75)

	#color bar for main plot
	ax3 = plt.subplot(gs[:,2])
	cbar = fig.colorbar(ctr, cax=ax3, ticks=np.linspace(minVal, maxVal, 11)) #11 labels on color bar
	cbar.set_label('% Change in Resistivity', fontsize=fs)
	ax3.tick_params(axis='both', which='major', labelsize=fs*0.75)
	ax3.tick_params(axis='both', which='minor', labelsize=fs*0.75)

	#depth profile
	ax1 = plt.subplot(gs[:,0])
	profileX, profileY = percentChange[x_2 == find_nearest(x_2, sliceLoc)], np.unique(z_2)
	ax1.scatter(profileX, profileY)
	ax1.plot(profileX, profileY)
	ax1.invert_yaxis()
	ax1.grid(which='both')
	ax1.set_ylabel('Depth [m]', fontsize=fs)
	ax1.set_xlabel('% Change', fontsize=fs)
	ax1.set_title('%% Change at %d m' %sliceLoc, fontsize=fs)
	ax1.tick_params(axis='both', which='major', labelsize=fs*0.75)
	ax1.tick_params(axis='both', which='minor', labelsize=fs*0.75)

	#depth slice
	ax2 = plt.subplot(gs[1,1])
	depthX, depthY = np.unique(x_2), percentChange[z_2 == find_nearest(z_2, depthSlice)]
	ax2.scatter(depthX, depthY)
	ax2.plot(depthX, depthY)
	ax2.grid(which='both')
	ax2.set_ylabel('% Change', fontsize=fs)
	ax2.set_xlabel('Profile Distance [m]', fontsize=fs)
	ax2.set_title('%% Change at z = %d' %depthSlice, fontsize=fs)
	ax2.set_xlim([np.amin(x_2), np.amax(x_2)])
	ax2.tick_params(axis='both', which='major', labelsize=fs*0.75)
	ax2.tick_params(axis='both', which='minor', labelsize=fs*0.75)

	if save ==False:
		supT = fig.suptitle(plotTitle, x=0.445, y=1.04, fontsize=14)

	figTitle = plotTitle.replace('.xyz','')

	if save:
		plt.savefig(figTitle, dpi=150, bbox_inches='tight')

	return



##################
#UTILITY CODE
##################

def find_nearest(array, value):
	#function for finding nearest value to depth and slice locations
	array = np.asarray(array)
	idx = (np.abs(array - value)).argmin()
	return array[idx]