MT_WGS, #Spatial dataframe in WGS coordinates to make map

v.MT, #variogram for MT dataset using all data

rv.MT, #Cressie's robust variogram for MT dataset using all data

Vcut, #variogram cutoff in m

Vbins, #Number of variogram bins

o, #index for the operator to be left out

LowLim = 2, #Lower limit for number of wells drilled by an operator. Need at least LowLim wells to make a plot for the operator.

HistSep = 10, #x-axis bar separartion on histogram

Histylim = 300, #y-axis upper limit on histogram

RegName, #Region name for figure

res = 300, #png resolution

height = 8, #png height in inches

width = 8, #png width in inches

MaxLagDist, #maximum lag separation distance for which to compute jackknife distance metrics

SensitivityDist = NA, #If supplied, will provide results of the metrics in the specified variogram lag distance interval for MaxLagDist

plt = TRUE # should diagnostic plots be made for each operator?

){

#Unique operators

UniOps = unique(MT$Operator)

#Only plot the operator if they have more than LowLim well. It doesn't make sense otherwise.

if(nrow(MT[MT$Operator == UniOps[o],]) >= LowLim){

#Operator left-out variograms

v.o = variogram(Qs~1, MT[-which(MT$Operator == UniOps[o]),], cutoff=Vcut, width = Vcut/Vbins)

rv.o = variogram(Qs~1, MT[-which(MT$Operator == UniOps[o]),], cutoff=Vcut, width = Vcut/Vbins, cressie = TRUE)

if (plt){

png(paste0(RegName, '_OperatorESDA_', o, '.png'), res = res, height = height, width = width, units = 'in')

#Histograms Overlain

#layout(rbind(c(1,4), c(2,3)))

#Histogram showing difference between current operator (red) and all data (black)

#hist(MT$Qs, col = 'black', main = UniOps[o], xlab = expression(paste('Heat Flow (mW/m'^2,')')), ylab = 'Frequency', ylim = c(0,Histylim), xlim = c(0,round(max(MT$Qs + HistSep/2),-1)), breaks = seq(round(min(MT$Qs - HistSep/2),-1),round(max(MT$Qs + HistSep/2),-1),HistSep))

#par(new=TRUE)

#hist(MT$Qs[-which(MT$Operator == UniOps[o])], border = 'red', axes = FALSE, xlab = '', main = '', ylab = '', ylim = c(0,Histylim), xlim = c(0,round(max(MT$Qs + HistSep/2),-1)), breaks = seq(round(min(MT$Qs - HistSep/2),-1),round(max(MT$Qs + HistSep/2),-1),HistSep))

#legend('topright', legend = c('All data', 'Operator Removed'), col = c('black', 'red'), pch = 15)

#Histograms on separate plots

#layout(rbind(c(1,5), c(2,5), c(3,4), c(3,4)))

#hist(MT$Qs, col = 'black', main = UniOps[o], xlab = expression(paste('Heat Flow (mW/m'^2,')')), ylab = 'Frequency', ylim = c(0,Histylim), xlim = c(0,round(max(MT$Qs + HistSep/2),-1)), breaks = seq(round(min(MT$Qs - HistSep/2),-1),round(max(MT$Qs + HistSep/2),-1),HistSep))

#par(new=TRUE)

#hist(MT$Qs[-which(MT$Operator == UniOps[o])], border = 'red', axes = FALSE, xlab = '', main = '', ylab = '', ylim = c(0,Histylim), xlim = c(0,round(max(MT$Qs + HistSep/2),-1)), breaks = seq(round(min(MT$Qs - HistSep/2),-1),round(max(MT$Qs + HistSep/2),-1),HistSep))

#legend('topright', legend = c('All data', 'Operator Removed'), col = c('black', 'red'), pch = 15)

#Boxplot

layout(rbind(c(1,4), c(2,3)))

boxplot(horizontal = TRUE, x = MT$Qs, main = UniOps[o], xlab = expression(paste('Heat Flow (mW/m'^2,')')), ylim = c(0,round(max(MT$Qs + HistSep/2),-1)), at = 1, xlim = c(0,2))

par(new = TRUE)

#All red

#plot(x = MT$Qs[which(MT$Operator == UniOps[o])], y = jitter(rep(1, length(MT$Qs[which(MT$Operator == UniOps[o])])), amount = 0.4), col = 'red', xlab = '', ylab = '', axes = FALSE, xlim = c(0,round(max(MT$Qs + HistSep/2),-1)), ylim = c(0,2))

#Colored by value of heat flow

plot(x = MT$Qs[which(MT$Operator == UniOps[o])], y = jitter(rep(1, length(MT$Qs[which(MT$Operator == UniOps[o])])), amount = 0.4), col = colFun(MT$Qs[which(MT$Operator == UniOps[o])]), xlab = '', ylab = '', axes = FALSE, xlim = c(0,round(max(MT$Qs + HistSep/2),-1)), ylim = c(0,2))

legend('topright', legend = c(paste("All Wells: N wells =", nrow(MT)), paste("Operator's Wells: N wells =", nrow(MT[which(MT$Operator == UniOps[o]),]))), col = c('black', 'red'), lty = c(1,NA), pch = c(NA,1))

#Variogram of region with and without operator

par(mar = c(5,5,2,1))

plot(v.MT$dist, v.MT$gamma, ylim = c(0,round(max(v.MT$gamma),-1)), xlim = c(0,Vcut), xlab = 'Separation Distance (m)', ylab = expression(paste('Semivariance (mW/m'^2,')'^2)), main = 'MOM Semi-variogram')

par(new=TRUE)

plot(v.o$dist, v.o$gamma, ylim = c(0,round(max(v.MT$gamma),-1)), xlim = c(0,Vcut), xlab = '', ylab = '', col = 'red')

legend('bottomright', legend = c('All Wells', 'Operator Removed'), col = c('black', 'red'), pch = 1)

#Robust variogram of region with and without operator

plot(rv.MT$dist, rv.MT$gamma, ylim = c(0,round(max(v.MT$gamma),-1)), xlim = c(0,Vcut), xlab = 'Separation Distance (m)', ylab = expression(paste('Semivariance (mW/m'^2,')'^2)), main = 'Robust Semi-variogram')

par(new=TRUE)

plot(rv.o$dist, rv.o$gamma, ylim = c(0,round(max(v.MT$gamma),-1)), xlim = c(0,Vcut), xlab = '', ylab = '', col = 'red')

legend('bottomright', legend = c('All Wells', 'Operator Removed'), col = c('black', 'red'), pch = 1)

#Map

plot(MT_WGS, pch = 16, cex = 0.4, col ='white')

plot(Counties[which(Counties$STATEFP == 42 | Counties$STATEFP == 36 | Counties$STATEFP == 54 | Counties$STATEFP == 51| Counties$STATEFP == 24| Counties$STATEFP == 21),], add=TRUE, border = 'grey')

plot(NY, lwd = 2, add=TRUE)

plot(PA, lwd = 2, add=TRUE)

plot(WV, lwd = 2, add=TRUE)

plot(MD, lwd = 2, add=TRUE)

plot(KY, lwd = 2, add=TRUE)

plot(VA, lwd = 2, add=TRUE)

north.arrow(-83.5, 37.8, 0.05, lab = 'N', cex.lab = 1.5, col='black', cex = 0.7)

degAxis(side = 2, seq(34, 46, 1), cex.axis = 1.5)

degAxis(side = 2, seq(34, 46, 1), labels = FALSE)

degAxis(side = 4, seq(34, 46, 1), labels = FALSE)

degAxis(side = 1, seq(-70, -86, -1), cex.axis = 1.5)

degAxis(side = 3, seq(-70, -86, -1), labels = FALSE)

degAxis(side = 1, seq(-70, -86, -1), labels = FALSE)

plot(MT_WGS, pch = 16, cex = 0.4, add = T)

#All red

#plot(MT_WGS[which(MT$Operator == UniOps[o]),], pch = 16, cex = 0.4, add = T, col = 'red')

#Colored by heat flow value. Deeper on top of shallower where there are overlaps in space.

plot(MT_WGS[which(MT$Operator == UniOps[o]),][order(MT_WGS[which(MT$Operator == UniOps[o]),]$WellDepth, decreasing = FALSE),], pch = 16, cex = 0.4, add = T, col = colFun(MT_WGS$Qs[which(MT_WGS$Operator == UniOps[o])][order(MT_WGS[which(MT$Operator == UniOps[o]),]$WellDepth, decreasing = FALSE)]))

legend('topleft', legend = c('Other Wells', "Operator's Wells"), col = c('black', 'red'), pch = 16, cex = 1)

dev.off()

}

#Number of operator wells

Nop = length(which(MT$Operator == UniOps[o]))

#Check if a sensitivity analysis on maximum variogram lag distance is requested.

if (!is.na(SensitivityDist)){

#Run sensitivity analysis

#Determine the distances to evaluate

Dists = seq(SensitivityDist, MaxLagDist, SensitivityDist)

#Compute metrics for each of the distances

#Get the lag that should be selected for cumulative sums

DistLags = DiffMOMW = DiffRobustW = vector('numeric', length(Dists))

for (d in 1:length(Dists)){

DistLags[d] = max(which(v.MT$dist <= Dists[d]))

#Weighting by the separation distance lags

#Note: not using different weights for MOM and Robust because want to compare the same number of lags for both

weights = v.MT$np[which(v.MT$dist <= Dists[d])]*(v.MT$dist[which(v.MT$dist <= Dists[d])]^(-2))

weights = weights/sum(weights)

DiffMOMW[d] = sum((v.MT$gamma - v.o$gamma)[v.MT$dist <= Dists[d]]*weights)

DiffRobustW[d] = sum((rv.MT$gamma - rv.o$gamma)[v.MT$dist <= Dists[d]]*weights)

}

#Return the sum of the difference between the variograms with and without operator over the spatial area that they operate

#Using only the full dataset to define the lags computed because when an operator is left out is may result in changing the lag distance slightly.

#Weighting by the difference per well

DiffMOM = cumsum((v.MT$gamma - v.o$gamma))[DistLags]/Nop

DiffRobust = cumsum((rv.MT$gamma - rv.o$gamma))[DistLags]/Nop

#This metric is not recommended because it can make increases in semi-variances seem like decreases.

#DiffMOM_Jack = sum((v.MT$gamma - v.o$gamma*(nrow(MT) - Nop)/nrow(MT))[v.MT$dist <= MaxLagDist])

#DiffRobust_Jack = sum((rv.MT$gamma - rv.o$gamma*(nrow(MT) - Nop)/nrow(MT))[rv.MT$dist <= MaxLagDist])

#This metric is better because it adjusts for number of points for each operator, but it may not be meaningful.

DiffMOM_Jack = cumsum((v.MT$gamma - v.o$gamma))[DistLags]*(nrow(MT) - Nop)/nrow(MT)

DiffRobust_Jack = cumsum((rv.MT$gamma - rv.o$gamma))[DistLags]*(nrow(MT) - Nop)/nrow(MT)

#Return the mean of the full dataset and the operator's dataset

DiffMeans = mean(MT$Qs) - mean(MT$Qs[which(MT$Operator == UniOps[o])])

t_diff = t.test(x = MT$Qs[-which(MT$Operator == UniOps[o])], y = MT$Qs[which(MT$Operator == UniOps[o])], alternative = 'two.sided', mu = 0)$p.value

w_diff = wilcox.test(x = MT$Qs[-which(MT$Operator == UniOps[o])], y = MT$Qs[which(MT$Operator == UniOps[o])], alternative = 'two.sided', mu = 0, correct = TRUE)$p.value

#Return the mean and scaled mean of the operator-left-out dataset

DiffMeanRmOp = mean(MT$Qs) - mean(MT$Qs[-which(MT$Operator == UniOps[o])])

ScaledDiffMeanRmOp = mean(MT$Qs) - ((nrow(MT) - Nop)*mean(MT$Qs[-which(MT$Operator == UniOps[o])]))/nrow(MT)

#2-sample t-test for difference in means. These are large samples. Reporting wilcoxon alternative anyway.

t_RmOp = t.test(x = MT$Qs, y = MT$Qs[-which(MT$Operator == UniOps[o])], alternative = 'two.sided', mu = 0)$p.value

w_RmOp = wilcox.test(x = MT$Qs, y = MT$Qs[-which(MT$Operator == UniOps[o])], alternative = 'two.sided', mu = 0, correct = TRUE)$p.value

lst = list(Op = UniOps[o], DiffMat = data.frame(DiffMOM = DiffMOM, DiffMOMW = DiffMOMW, DiffMOM_Jack = DiffMOM_Jack, DiffRobust = DiffRobust, DiffRobustW = DiffRobustW, DiffRobust_Jack = DiffRobust_Jack), DiffMeans = DiffMeans, pt_DiffMean = t_diff, pw_DiffMean = w_diff, DiffMeanRmOp = DiffMeanRmOp, pt_DiffRmOp = t_RmOp, pw_DiffRmOp = w_RmOp, ScaledDiffMeanRmOp = ScaledDiffMeanRmOp, NumWells = Nop)

}else{

#Just report one result supplied by MaxLagDist

#Return the sum of the difference between the variograms with and without operator over the spatial area that they operate

#Using only the full dataset to define the lags computed because when an operator is left out is may result in changing the lag distance slightly.

#Weighting by the difference per well

DiffMOM = sum((v.MT$gamma - v.o$gamma)[v.MT$dist <= MaxLagDist])/Nop

DiffRobust = sum((rv.MT$gamma - rv.o$gamma)[rv.MT$dist <= MaxLagDist])/Nop

#Weighting by the separation distance lags

weights = v.MT$np[which(v.MT$dist <= MaxLagDist)]*(v.MT$dist[which(v.MT$dist <= MaxLagDist)]^(-2))

weights = weights/sum(weights)

DiffMOMW = sum((v.MT$gamma - v.o$gamma)[v.MT$dist <= MaxLagDist]*weights)

DiffRobustW = sum((rv.MT$gamma - rv.o$gamma)[rv.MT$dist <= MaxLagDist]*weights)

#This metric is not recommended because it can make increases in semi-variances seem like decreases.

#DiffMOM_Jack = sum((v.MT$gamma - v.o$gamma*(nrow(MT) - Nop)/nrow(MT))[v.MT$dist <= MaxLagDist])

#DiffRobust_Jack = sum((rv.MT$gamma - rv.o$gamma*(nrow(MT) - Nop)/nrow(MT))[rv.MT$dist <= MaxLagDist])

#This metric is better because it adjusts for number of points for each operator, but it may not be meaningful.

DiffMOM_Jack = sum((v.MT$gamma - v.o$gamma)[v.MT$dist <= MaxLagDist])*(nrow(MT) - Nop)/nrow(MT)

DiffRobust_Jack = sum((rv.MT$gamma - rv.o$gamma)[rv.MT$dist <= MaxLagDist])*(nrow(MT) - Nop)/nrow(MT)

#Return the mean of the full dataset and the operator's dataset

DiffMeans = mean(MT$Qs) - mean(MT$Qs[which(MT$Operator == UniOps[o])])

t_diff = t.test(x = MT$Qs[-which(MT$Operator == UniOps[o])], y = MT$Qs[which(MT$Operator == UniOps[o])], alternative = 'two.sided', mu = 0)$p.value

w_diff = wilcox.test(x = MT$Qs[-which(MT$Operator == UniOps[o])], y = MT$Qs[which(MT$Operator == UniOps[o])], alternative = 'two.sided', mu = 0, correct = TRUE)$p.value

#Return the mean and scaled mean of the operator-left-out dataset

DiffMeanRmOp = mean(MT$Qs) - mean(MT$Qs[-which(MT$Operator == UniOps[o])])

ScaledDiffMeanRmOp = mean(MT$Qs) - ((nrow(MT) - Nop)*mean(MT$Qs[-which(MT$Operator == UniOps[o])]))/nrow(MT)

#2-sample t-test for difference in means. These are large samples. Reporting wilcoxon alternative anyway.

t_RmOp = t.test(x = MT$Qs, y = MT$Qs[-which(MT$Operator == UniOps[o])], alternative = 'two.sided', mu = 0)$p.value

w_RmOp = wilcox.test(x = MT$Qs, y = MT$Qs[-which(MT$Operator == UniOps[o])], alternative = 'two.sided', mu = 0, correct = TRUE)$p.value

lst = data.frame(Op = UniOps[o], OpNum = o, DiffMOM = DiffMOM, DiffMOMW = DiffMOMW, DiffMOM_Jack = DiffMOM_Jack, DiffRobust = DiffRobust, DiffRobustW = DiffRobustW, DiffRobust_Jack = DiffRobust_Jack, DiffMeans = DiffMeans, pt_DiffMean = t_diff, pw_DiffMean = w_diff, DiffMeanRmOp = DiffMeanRmOp, pt_DiffRmOp = t_RmOp, pw_DiffRmOp = w_RmOp, ScaledDiffMeanRmOp = ScaledDiffMeanRmOp, NumWells = Nop, stringsAsFactors = FALSE)

}

return(lst)

}

){

DiffMOM = DiffMOMW = DiffRobust = DiffRobustW = matrix(0, nrow = nrow(OpDiag), ncol = nrow(OpDiag[1,2][[1]]))

for (i in 1:nrow(OpDiag)){

DiffMOM[i,] = OpDiag[i,2][[1]]$DiffMOM

DiffMOMW[i,] = OpDiag[i,2][[1]]$DiffMOMW

DiffRobust[i,] = OpDiag[i,2][[1]]$DiffRobust

DiffRobustW[i,] = OpDiag[i,2][[1]]$DiffRobustW

}

RankOpDiag_MOMmat = apply(X = DiffMOM, MARGIN = 2, FUN = order, decreasing = TRUE)

RankOpDiag_MOMWmat = apply(X = DiffMOMW, MARGIN = 2, FUN = order, decreasing = TRUE)

RankOpDiag_Robustmat = apply(X = DiffRobust, MARGIN = 2, FUN = order, decreasing = TRUE)

RankOpDiag_RobustWmat = apply(X = DiffRobustW, MARGIN = 2, FUN = order, decreasing = TRUE)

return(list(RankOpDiag_MOMmat = RankOpDiag_MOMmat, RankOpDiag_MOMWmat = RankOpDiag_MOMWmat, RankOpDiag_Robustmat = RankOpDiag_Robustmat, RankOpDiag_RobustWmat = RankOpDiag_RobustWmat, DiffMOM = DiffMOM, DiffMOMW = DiffMOMW, DiffRobust = DiffRobust, DiffRobustW =DiffRobustW))

res = 300, #png resolution

height = 10, #png height in inches

width = 7, #png width in inches

xlim = c(1,60), #x axis limits of plot

NumOps = 100, #y axis number of operators to show

NumOpsStep = 20, #y axis tick mark locations for operators

xStep = 2.5, #x-axis step size

PlotName, #name to be added to the plot file

cols = NA #color palette function for the plots

){

if(is.na(cols)){

cols = colorRampPalette(brewer.pal(7, name = 'PuOr'))

}

png(paste0(PlotName, '_OpDiagnostics_DistanceSensitivity.png'), res = res, width = width, height = height, units = 'in')

layout(cbind(c(1,2,3,4)))

par(mar = c(2,5,3,1))

for (i in 1:nrow(OpDiagRanks$RankOpDiag_MOMmat)){

if (i == 1){

plot(seq(xStep,max(xlim),xStep), which(OpDiagRanks$RankOpDiag_MOMmat == i) - (seq(0,ncol(OpDiagRanks$RankOpDiag_MOMmat)-1,1)*nrow(OpDiagRanks$RankOpDiag_MOMmat)),

type = 'o', lty = 1, pch = 16,

ylim = rev(c(0,nrow(OpDiagRanks$RankOpDiag_MOMmat))), xlim = xlim,

col = cols(nrow(OpDiagRanks$RankOpDiag_MOMmat))[which(OpDiagRanks$RankOpDiag_MOMmat[,1] == i)],

ylab = 'Operator Rank', xlab = '', main = 'MOM Semi-variogram, Weighted by Number of Wells',

cex.lab = 1.5, cex.main = 2, axes = FALSE)

box()

axis(side = 2, at = seq(0,NumOps,NumOpsStep), labels = TRUE, cex.axis = 1.5)

axis(side = 4, at = seq(0,NumOps,NumOpsStep), labels = FALSE)

axis(side = 1, at = seq(xStep,max(xlim),xStep), labels = FALSE)

}else{

plot(seq(xStep,max(xlim),xStep), which(OpDiagRanks$RankOpDiag_MOMmat == i) - (seq(0,ncol(OpDiagRanks$RankOpDiag_MOMmat)-1,1)*nrow(OpDiagRanks$RankOpDiag_MOMmat)),

type = 'o', lty = 1, pch = 16,

ylim = rev(c(0,nrow(OpDiagRanks$RankOpDiag_MOMmat))), xlim = xlim,

col = cols(nrow(OpDiagRanks$RankOpDiag_MOMmat))[which(OpDiagRanks$RankOpDiag_MOMmat[,1] == i)],

axes = FALSE, xlab = '', ylab = '')

}

par(new = T)

}

par(new = FALSE)

for (i in 1:nrow(OpDiagRanks$RankOpDiag_MOMmat)){

if (i == 1){

plot(seq(xStep,max(xlim),xStep), which(OpDiagRanks$RankOpDiag_Robustmat == i) - (seq(0,ncol(OpDiagRanks$RankOpDiag_Robustmat)-1,1)*nrow(OpDiagRanks$RankOpDiag_Robustmat)),

type = 'o', lty = 1, pch = 16,

ylim = rev(c(0,nrow(OpDiagRanks$RankOpDiag_Robustmat))), xlim = xlim,

col = cols(nrow(OpDiagRanks$RankOpDiag_Robustmat))[which(OpDiagRanks$RankOpDiag_Robustmat[,1] == i)],

ylab = 'Operator Rank', xlab = '', main = 'Robust Semi-variogram, Weighted by Number of Wells',

cex.lab = 1.5, cex.main = 2, axes = FALSE)

box()

axis(side = 2, at = seq(0,NumOps,NumOpsStep), labels = TRUE, cex.axis = 1.5)

axis(side = 4, at = seq(0,NumOps,NumOpsStep), labels = FALSE)

axis(side = 1, at = seq(xStep,max(xlim),xStep), labels = FALSE)

}else{

plot(seq(xStep,max(xlim),xStep), which(OpDiagRanks$RankOpDiag_Robustmat == i) - (seq(0,ncol(OpDiagRanks$RankOpDiag_Robustmat)-1,1)*nrow(OpDiagRanks$RankOpDiag_Robustmat)),

type = 'o', lty = 1, pch = 16,

ylim = rev(c(0,nrow(OpDiagRanks$RankOpDiag_Robustmat))), xlim = xlim,

col = cols(nrow(OpDiagRanks$RankOpDiag_Robustmat))[which(OpDiagRanks$RankOpDiag_Robustmat[,1] == i)],

axes = FALSE, xlab = '', ylab = '')

}

par(new = T)

}

par(new = FALSE)

for (i in 1:nrow(OpDiagRanks$RankOpDiag_MOMmat)){

if (i == 1){

plot(seq(xStep,max(xlim),xStep), which(OpDiagRanks$RankOpDiag_MOMWmat == i) - (seq(0,ncol(OpDiagRanks$RankOpDiag_MOMWmat)-1,1)*nrow(OpDiagRanks$RankOpDiag_MOMWmat)),

type = 'o', lty = 1, pch = 16,

ylim = rev(c(0,nrow(OpDiagRanks$RankOpDiag_MOMWmat))), xlim = xlim,

col = cols(nrow(OpDiagRanks$RankOpDiag_MOMWmat))[which(OpDiagRanks$RankOpDiag_MOMWmat[,1] == i)],

ylab = 'Operator Rank', xlab = '', main = 'MOM Semi-variogram, Weighted by Separation Distance',

cex.lab = 1.5, cex.main = 2, axes = FALSE)

box()

axis(side = 2, at = seq(0,NumOps,NumOpsStep), labels = TRUE, cex.axis = 1.5)

axis(side = 4, at = seq(0,NumOps,NumOpsStep), labels = FALSE)

axis(side = 1, at = seq(xStep,max(xlim),xStep), labels = FALSE)

}else{

plot(seq(xStep,max(xlim),xStep), which(OpDiagRanks$RankOpDiag_MOMWmat == i) - (seq(0,ncol(OpDiagRanks$RankOpDiag_MOMWmat)-1,1)*nrow(OpDiagRanks$RankOpDiag_MOMWmat)),

type = 'o', lty = 1, pch = 16,

ylim = rev(c(0,nrow(OpDiagRanks$RankOpDiag_MOMWmat))), xlim = xlim,

col = cols(nrow(OpDiagRanks$RankOpDiag_MOMWmat))[which(OpDiagRanks$RankOpDiag_MOMWmat[,1] == i)],

axes = FALSE, xlab = '', ylab = '')

}

par(new = T)

}

par(new = FALSE)

for (i in 1:nrow(OpDiagRanks$RankOpDiag_MOMmat)){

if (i == 1){

plot(seq(xStep,max(xlim),xStep), which(OpDiagRanks$RankOpDiag_RobustWmat == i) - (seq(0,ncol(OpDiagRanks$RankOpDiag_RobustWmat)-1,1)*nrow(OpDiagRanks$RankOpDiag_RobustWmat)),

type = 'o', lty = 1, pch = 16,

ylim = rev(c(0,nrow(OpDiagRanks$RankOpDiag_RobustWmat))), xlim = xlim,

col = cols(nrow(OpDiagRanks$RankOpDiag_RobustWmat))[which(OpDiagRanks$RankOpDiag_RobustWmat[,1] == i)],

ylab = 'Operator Rank', xlab = 'Separation Distance (km)', main = 'Robust Semi-variogram, Weighted by Separation Distance',

cex.lab = 1.5, cex.main = 2, axes = FALSE)

box()

axis(side = 2, at = seq(0,NumOps,NumOpsStep), labels = TRUE, cex.axis = 1.5)

axis(side = 4, at = seq(0,NumOps,NumOpsStep), labels = FALSE)

axis(side = 1, at = seq(2.5,57.5,5), labels = FALSE)

axis(side = 1, at = seq(5,60,5), labels = seq(5,60,5), cex.axis = 1.5)

}else{

plot(seq(xStep,max(xlim),xStep), which(OpDiagRanks$RankOpDiag_RobustWmat == i) - (seq(0,ncol(OpDiagRanks$RankOpDiag_RobustWmat)-1,1)*nrow(OpDiagRanks$RankOpDiag_RobustWmat)),

type = 'o', lty = 1, pch = 16,

ylim = rev(c(0,nrow(OpDiagRanks$RankOpDiag_RobustWmat))), xlim = xlim,

col = cols(nrow(OpDiagRanks$RankOpDiag_RobustWmat))[which(OpDiagRanks$RankOpDiag_RobustWmat[,1] == i)],

axes = FALSE, xlab = '', ylab = '')

}

par(new = T)

}

dev.off()

OpDiagRanks, #Output from the GetOperatorRanks function

col, #column to plot in OpDiagRanks

res = 300, #png resolution

heightPanel = 5, #png height in inches

widthPanel = 10, #png width in inches

heightParAx = 5, #png height in inches

widthParAx = 10, #png width in inches

NumOps = 100, #y axis number of operators to show

NumOpsStep = 20, #y axis tick mark locations for operators

PlotName, #name to be added to the plot file

cols = NA, #color palette function for the parallel axis plots

xStep #separation distance step size (km) for variogram. Used to place captions on plots

){

#Plot for variance represented by robust variogram difference, weighted by number of points the operator has

# Colored by the bias, represented by the p-value for difference in center for operator data vs. rest of data.

png(paste0(PlotName, '_OpDiagnostics_PanelPlot.png'), res = res, width = widthPanel, height = heightPanel, units = 'in')

par(mar = c(5,5,5,3))

layout(rbind(c(1,2)))

plot(y = OpDiagRanks$DiffMOM[,col], x = OpDiagRanks$DiffMOMW[,col],

col = colFun(as.numeric(OpDiag[,"pw_DiffMean"])),

main = 'Difference in Variograms: \n All Data - Operator Removed \n 0 - 30 km Separation', ylab = 'MOM, weight = 1/number operator wells', xlab = 'MOM, weight ~ 1/(distance lag)^2')

legend('topleft', title = 'p-value', legend = c('< 0.05', '< 0.10', '< 0.15', '< 0.20', '> 0.20'), col = colFun(c(0,0.05,0.10,0.15,0.20,0.21)), pch = 1)

plot(y = OpDiagRanks$DiffMOM[,col], x = OpDiagRanks$DiffRobust[,col],

col = colFun(as.numeric(OpDiag[,"pw_DiffMean"])),

main = 'Difference in Variograms: \n All Data - Operator Removed \n 0 - 30 km Separation', ylab = 'MOM, weight = 1/number operator wells', xlab = 'Robust, weight = 1/number operator wells')

legend('topleft', title = 'p-value', legend = c('< 0.05', '< 0.10', '< 0.15', '< 0.20', '> 0.20'), col = colFun(c(0,0.05,0.10,0.15,0.20,0.21)), pch = 1)

dev.off()

#plot(OpDiagRanks$DiffRobust, OpDiagRanks$DiffRobustW, col = colFun(OpDiagRanks$pw_DiffMean), main = 'Difference in MOM Semi-variance: 0 to 20 km Separation', xlab = 'Difference in MOM semi-variance: 0 - 20 km Separation', ylab = 'Difference in Robust semi-variance: 0 - 20 km Separation')

#plot(OpDiagRanks$DiffMOMW, OpDiagRanks$DiffRobustW, col = colFun(OpDiagRanks$pw_DiffMean), main = 'Difference in MOM Semi-variance: 0 to 20 km Separation', xlab = 'Difference in MOM semi-variance: 0 - 20 km Separation', ylab = 'Difference in Robust semi-variance: 0 - 20 km Separation')

#Univariate plots

#plot(seq(1,length(OpDiagRanks$DiffMOM),1), abs(sort(as.numeric(OpDiagRanks$DiffMOM))), main = 'Difference in MOM Semi-variance: 0 to 20 km Separation', xlab = 'Sorted Operator ID', ylab = 'Difference')

#plot(seq(1,length(OpDiagRanks$DiffMOM),1), abs(sort(OpDiagRanks$DiffRobust)), main = 'Difference in Robust Semi-variance: 0 to 20 km Separation', xlab = 'Sorted Operator ID', ylab = 'Difference')

#Make a parallel axis plot of the rank of the operator for each metric

#Make a matrix for the ranks

RankMat = cbind(OpDiagRanks$RankOpDiag_MOMmat[,col],

OpDiagRanks$RankOpDiag_Robustmat[,col],

OpDiagRanks$RankOpDiag_MOMWmat[,col],

OpDiagRanks$RankOpDiag_RobustWmat[,col],

order(as.numeric(OpDiag[,"pw_DiffMean"]), decreasing = FALSE))

if(is.na(cols)){

cols = colorRampPalette(brewer.pal(7, name = 'PuOr'))

}

png(paste0(PlotName, '_OpDiagnostics_ParallelAxis_Ranks.png'), res = res, width = widthParAx, height = heightParAx, units = 'in')

par(mar = c(5,5,1,3))

for (i in 1:nrow(RankMat)){

if (i == 1){

plot(seq(1,ncol(RankMat),1), which(RankMat == i) - (seq(0,ncol(RankMat)-1,1)*nrow(RankMat)),

type = 'o', lty = 1, pch = 16,

ylim = rev(c(0,nrow(RankMat))), xlim = c(1,ncol(RankMat)),

col = cols(nrow(RankMat))[which(RankMat[,1] == i)],

ylab = 'Rank of Operator for Metric', xlab = '', cex.lab = 1.5, axes = FALSE)

box()

axis(side = 2, at = seq(0,NumOps,NumOpsStep), labels = TRUE, cex.axis = 1.5)

axis(side = 1, at = seq(1,ncol(RankMat),1), line = 2, tick = FALSE, padj = -0.3,

labels = c(paste0('MOM Variogram \n 0 - ', col*xStep, ' km \n w = N wells'),

paste0('Robust Variogram \n 0 - ', col*xStep, ' km \n w = N wells'),

paste0('MOM Variogram \n 0 - ', col*xStep, ' km \n w = spatial lag'),

paste0('Robust Variogram \n 0 - ', col*xStep, ' km \n w = spatial lag'),

'p-value \n Wilcoxon Rank Sum \n'))

}else{

plot(seq(1,ncol(RankMat),1), which(RankMat == i) - (seq(0,ncol(RankMat)-1,1)*nrow(RankMat)),

type = 'o', lty = 1, pch = 16,

ylim = rev(c(0,nrow(RankMat))), xlim = c(1,ncol(RankMat)), col = cols(nrow(RankMat))[which(RankMat[,1] == i)],

axes = FALSE, xlab = '', ylab = '')

}

par(new = T)

}

dev.off()

#Make a matrix for the ranks

NormMat = cbind((OpDiagRanks$DiffMOM[,col] - min(OpDiagRanks$DiffMOM[,col])) / (max(OpDiagRanks$DiffMOM[,col]) - min(OpDiagRanks$DiffMOM[,col])),

(OpDiagRanks$DiffRobust[,col] - min(OpDiagRanks$DiffRobust[,col])) / (max(OpDiagRanks$DiffRobust[,col]) - min(OpDiagRanks$DiffRobust[,col])),

(OpDiagRanks$DiffMOMW[,col] - min(OpDiagRanks$DiffMOMW[,col])) / (max(OpDiagRanks$DiffMOMW[,col]) - min(OpDiagRanks$DiffMOMW[,col])),

(OpDiagRanks$DiffRobustW[,col] - min(OpDiagRanks$DiffRobustW[,col])) / (max(OpDiagRanks$DiffRobustW[,col]) - min(OpDiagRanks$DiffRobustW[,col])),

as.numeric(OpDiag[,"pw_DiffMean"]))

png(paste0(PlotName, '_OpDiagnostics_ParallelAxis_Normalized.png'), res = res, width = widthParAx, height = heightParAx, units = 'in')

par(mar = c(5,5,1,3))

for (i in 1:nrow(NormMat)){

if (i == 1){

plot(seq(1,ncol(NormMat),1), NormMat[i,],

type = 'o', lty = 1, pch = 16,

ylim = c(0,1), xlim = c(1,ncol(NormMat)),

col = cols(nrow(RankMat))[which(RankMat[,1] == i)],

ylab = 'Normalized Value of Metric', xlab = '', cex.lab = 1.5, axes = FALSE)

box()

axis(side = 2, at = seq(0,1,.1), labels = TRUE, cex.axis = 1.5)

axis(side = 1, at = seq(1,ncol(RankMat),1), line = 2, tick = FALSE, padj = -0.3,

labels = c(paste0('MOM Variogram \n 0 - ', col*xStep, ' km \n w = N wells'),

paste0('Robust Variogram \n 0 - ', col*xStep, ' km \n w = N wells'),

paste0('MOM Variogram \n 0 - ', col*xStep, ' km \n w = spatial lag'),

paste0('Robust Variogram \n 0 - ', col*xStep, ' km \n w = spatial lag'),

'p-value \n Wilcoxon Rank Sum \n'))

}else{

plot(seq(1,ncol(NormMat),1), NormMat[i,],

type = 'o', lty = 1, pch = 16,

ylim = c(0,1), xlim = c(1,ncol(NormMat)), col = cols(nrow(RankMat))[which(RankMat[,1] == i)],

axes = FALSE, xlab = '', ylab = '')

}

par(new = T)

}

dev.off()

OpRanks, #One metric's dataframe output from the GetOperatorRanks function

col, #column of OpRanks to select

rank #Rank to select. 1 is highest value in col

){

RankOp = OpRanks[rank,col]

as.numeric(strsplit(names(OpDiag[,1][RankOp]), split = '.', fixed = TRUE)[[1]][2])

OpRanks, #Output from the GetOperatorRanks function

col, #column of OpRanks to select

N, #Ranks to select. The top N are selected for each metric

ID = TRUE, #Should plot ID number be returned? If FALSE, returns name of operator

WellData #Needs a column named Operator to return the name of the operator.

){

TopN = vector('numeric')

for (i in 1:N){

a = ReturnOperatorPlotNum(OpDiag, OpRanks$RankOpDiag_MOMmat, col = 12, rank = i)

b = ReturnOperatorPlotNum(OpDiag, OpRanks$RankOpDiag_MOMWmat, col = 12, rank = i)

c = ReturnOperatorPlotNum(OpDiag, OpRanks$RankOpDiag_Robustmat, col = 12, rank = i)

d = ReturnOperatorPlotNum(OpDiag, OpRanks$RankOpDiag_RobustWmat, col = 12, rank = i)

TopN = unique(c(a,b,c,d,TopN))

}

if(ID){

#PlotID numbers

return(sort(TopN))

}else{

#Operator Names

return(unique(WellData$Operator)[sort(TopN)])

}