forked from GeostatsGuy/2DayCourse_Exercises
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathvariogram_demo.R
138 lines (101 loc) · 7.96 KB
/
variogram_demo.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
# Variogram Analysis in R for Engineers and Geoscientists New to R
# Michael Pyrcz, University of Texas at Austin, Twitter @GeostatsGuy
# This will be used in my Introduction to Geostatistics Undergraduate Class
# It is assumed that students have no previous R experience.
# This utilizes the gstat library by Edzer Pedesma, appreciation to Dr. Pedesma for assistance.
# Load the required libraries, you may have to first go to "Tools/Install Packages..." to install these first
library(gstat) # geostatistical methods by Edzer Pebesma
library(sp) # spatial points addition to regular data frames
library(plyr) # splitting, applying and combining data by Hadley Wickham
library(ggplot2)
# Declare functions, this function completes standard normal transform on a dataset
nscore <- function(x) { # written by Ashton Shortridge, May/June, 2008
# Takes a vector of values x and calculates their normal scores. Returns
# a list with the scores and an ordered table of original values and
# scores, which is useful as a back-transform table. See backtr().
nscore <- qqnorm(x, plot.it = FALSE)$x # normal score
trn.table <- data.frame(x=sort(x),nscore=sort(nscore))
return (list(nscore=nscore, trn.table=trn.table))
}
# Set the working directory, I always like to do this so I don't lose files and to simplify subsequent read and writes
setwd("C:/PGE337")
# Read the data table from a comma delimited file - data on GitHub/GeostatsGuy/GeoDataSets
mydata = read.csv("2D_MV_200Wells.csv") # read in comma delimited data file
# Let's visualize the first several rows of our data so we can make sure we successfully loaded it
head(mydata) # show the first several rows of a data table in the console
# The columns are variables with variable names at the top and the rows are samples
# Convert the dataframe to a spatial points dataframe
class(mydata) # confirms that it is a dataframe
coordinates(mydata) = ~X+Y # indicate the X, Y spatial coordinates
summary(mydata) # confirms that it is now a spatial points dataframe
head(coordinates(mydata)) # check the first several coordinates
# Normal scores transform of the porosity data
npor.trn = nscore(mydata$porosity) # normal scores transform
mydata[["NPorosity"]]<-npor.trn$nscore # append the normal scores transform into the spatial data table
head(mydata)
# Visualize the dataset with a bubble plot
par(mfrow=c(2,2)) # set up a 2x2 matrix of plots
# Now let's check the porosity data distribution
hist(mydata$porosity,main="Porosity (%)",xlab="Porosity (%)",nclass = 15) # hist builds a regular frequency histogram
# ecdf makes a cdf object and plot command plots it
plot(ecdf(mydata$porosity),main="Porosity",xlab="Porosity (%",ylab="Cumulative Probability")
# Now let's check the normal score transform of the porosity data distribution
hist(mydata$NPorosity,main="N[Porosity (%)]",xlab="N[Porosity (%)]",nclass = 15) # hist builds a regular frequency histogram
# ecdf makes a cdf object and plot command plots it
plot(ecdf(mydata$NPorosity),main="N[Porosity]",xlab="N[Porosity (%)]",ylab="Cumulative Probability")
# We can plot the porosity data location map
par(mfrow=c(2,1)) # set up a 2x2 matrix of plots
spplot(mydata, "porosity", do.log = TRUE, # location map of porosity data
key.space=list(x=1.05,y=0.97,corner=c(0,1)),
scales=list(draw=T),xlab = "X (m)", ylab = "Y (m)",main ="Porosity (%)")
# We can also plot a simple bubble plot of the porosity data
bubble(mydata, "porosity", fill = FALSE, maxsize = 2, main ="Porosity (%)",identify = FALSE,xlab = "X (m)", ylab = "Y (m)")
# Calculate and fit an isotropic model to the normal score porosity variogram
por.vg.iso = variogram(NPorosity~1,mydata,cutoff = 3000,width =400,alpha = 35.0,tol.hor=90.0) # isotropic N[porosity]
plot(por.vg.iso,scaled = FALSE,var.lines = TRUE,main="Porosity Isotropic Variogram")
por.vm.iso <- vgm(psill = 0.6, "Exp", 800, anis = c(0, 0, 0, 1.0, 1.0),nugget=0.4)
por.vm.iso # check the variogram model parameters
plot(por.vg.iso,por.vm.iso,main="Porosity Isotropic Variogram") # use the built in gstat variogram plot
# Let's calculate variogram maps to identify possible directionality
par(mfrow=c(2,2)) # set up a 2x32 matrix of plots
plot(variogram(NPorosity~1,mydata, cutoff=1500, width=400, map=TRUE),main = "Semivariogram Map",max=1.0)
plot(variogram(NPorosity~1,mydata, cutoff=1500, width=400, map=TRUE),main = "Number of Pairs",np=TRUE)
por.vg.035 = variogram(NPorosity~1,mydata,cutoff = 3000,width =500,alpha = 35.0,tol.hor=22.5) # Calculate directional variograms of N[Porosity]
por.vg.125 = variogram(NPorosity~1,mydata,cutoff = 3000,width =500,alpha = 125.0,tol.hor=22.5)
plot(por.vg.035,main="Porosity Anisotropic 035 Variogram")
plot(por.vg.125,main="Porosity Anisotropic 125 Variogram")
# Anisotropy is parameterized as c(azimuth,dip,plunge,hratio,vratio) in #3D and c(azimuth,hratio) in 2D.
por.vm.ani <- vgm(psill = 0.6, "Exp", 800, anis = c(035, 0.5),nugget=0.4)
por.vm.ani # check the variogram model parameters
plot(por.vg.035,por.vm.ani,main="Porosity Anisotropic 035 Variogram") # use the built in gstat variogram plot
plot(por.vg.125,por.vm.ani,main="Porosity Anisotropic 125 Variogram") # use the built in gstat variogram plot
# Use auto-fitting to try to improve the variogram model, then check model and plot
por.vm.ani.auto <- fit.variogram(por.vg.iso,por.vm.ani,fit.sills = FALSE)
por.vm.ani.auto # check the autofit parameters and compare to our fit
plot(por.vg.035,por.vm.ani.auto,main="Porosity Anisotropic 035 Variogram")
plot(por.vg.125,por.vm.ani.auto,main="Porosity Anisotropic 125 Variogram")
# Let's make our own custom variogram plots
# The builtin variogram models do no fit into matrix plots, and are infelxible (e.g. no sill included)
par(mfrow=c(2,2)) # set up a 2x32 matrix of plots
name = c("Iso","035","125") # make name matrix
color = c("black","blue","red") # make color matrix
# Isotropic plot
plot(por.vg.iso$dist,por.vg.iso$gamma,main="Porosity Isotropic Variogram",xlab=" Lag Distance (m) ",ylab=" Semivariogram ", col=color[1],ylim=c(0,1.2))
abline(h = 1.0)
lines(por.vm.iso$model)
unit_vector = c(1,0,0) # unit vector doesn't matter since isotropic
vm.iso <- variogramLine(por.vm.iso,maxdist=3000,min=0.0001,n=100,dir=unit_vector,covariance=FALSE) # calculate variogram model
lines(vm.iso$dist,vm.iso$gamma,col=color[1]) # include variogram model
legend(2000,.8,name, cex=0.8, col=color,pch=c(21,21,21),lty=c(1,1,1)) # add legend
# Anisotropic plot
plot(por.vg.035$dist,por.vg.035$gamma,main="Porosity Anisotropic Variogram",xlab=" Lag Distance (m) ",ylab=" Semivariogram ", col=color[2],ylim=c(0,1.2))
points(por.vg.125$dist,por.vg.125$gamma,col=color[3])
abline(h = 1.0)
lines(por.vm.iso$model)
unit_vector = c(sin(35*pi/180),cos(35*pi/180),0) # unit vector for 035 azimuth
vm.ani.035 <- variogramLine(por.vm.ani,maxdist=3000,min=0.0001,n=100,dir=unit_vector,covariance=FALSE) # calculate 035 variogram model
lines(vm.ani.035$dist,vm.ani.035$gamma,col=color[2]) # include variogram model
unit_vector = c(sin(55*pi/180),-1*cos(35*pi/180),0) # unit vector for 125 azimuth
vm.ani.125 <- variogramLine(por.vm.ani,maxdist=3000,min=0.0001,n=100,dir=unit_vector,covariance=FALSE) # calculate 125 variogram model
lines(vm.ani.125$dist,vm.ani.125$gamma,col=color[3]) # include variogram model
legend(2000,.8,name, cex=0.8, col=color,pch=c(21,21,21),lty=c(1,1,1)) # add legend