R

Conditional simulation

Example

bore <- read.table("borehole.txt", header=T)
v <- variogram(log(lambda)~1, ~x+y, bore)
m <- vgm(.003, "Exp", 30)
plot(v, model=m)
sim <- krige(formula=log(lambda)~1, locations=~x+y, model=m, data=bore, newdata=dom, nmax=15, beta=1.16, nsim=4)
levelplot(exp(z)~x+y|name, map.to.lev(sim, z=c(3:6)), aspect="iso")
levelplot(z~x+y|name, map.to.lev(sim, z=c(3:6)), aspect="iso")

3D SGSIM

library(gstat)
library(lattice)
range1d <- seq(from=1, to=100, length=50)
grid3d <- expand.grid(x=range1d, y=range1d, z=range1d)
coordinates(grid3d) <- ~x+y+z
bore3d <- expand.grid(x=c(1,100), y=50, z=range1d)
bore3d <- cbind(bore3d, lambda=3.2)
coordinates(bore3d) <- ~x+y+z
res3D <- krige(formula = log(lambda) ~ 1, bore3d, grid3d, model = vgm(0.03, "Exp", 30), beta=1.16, nmax=15, nsim=1, debug.level=-1)
levelplot(exp(sim1) ~ x + y | z, as.data.frame(res3D))
summary(exp(as.data.frame(res3D)$sim1))

Histogram

library(lattice)
histogram(meuse$lead, nint=12)

hist(mhw$grain, breaks=seq(2.6, 5.2, by=0.1), col="lightblue", border="red", main="Grain yield")

h <- hist(mhw$grain, breaks = seq(2.6, 5.2, by = 0.2), plot = F)
plot(h, col = heat.colors(length(h$mids))[length(h$count) -
     rank(h$count) + 1], ylim = c(0, max(h$count) + 5),
     main = "Frequency histogram, Mercer & Hall grain yield",
     sub = "Counts shown above bar, actual values shown with rug plot",
     xlab = "Grain yield, lbs. per plot")
text(h$mids, h$count, h$count, pos = 3)

hist(mhw$grain, breaks = 20, col = "lightblue", border = "red", freq = F, main = "Grain yield, lbs per plot")
rug(mhw$grain)
lines(density(mhw$grain))

Empirical cumulative distribution function

plot(ecdf(mhw$grain), pch = 1, xlab = "Mercer & Hall, Grain yield, lbs. per plot",
     ylab = "Cumulative proportion of plots", main = "Empirical CDF",
     sub = "Quantiles shown with vertical lines")

Cumulative distribution function

x<-rnorm(1000)
h<-hist(x)
h$counts<-cumsum(h$counts)
plot.histogram(h)

Graph

Line

lines(density(mhw$grain))

Straight line

abline(v = median(grain), col = "blue")

Title

title("Straw yield predicted by grain yield")

Plot

plot(mhw$grain, mhw$straw)
pts <- identify(mhw$grain, mhw$straw)

Variogram

3D

library(gstat)
demo(gstat3D)
 
# $Id: gstat3D.R,v 1.4 2006-02-10 19:05:02 edzer Exp $
# simple demo of 3D interpolation of 50 points with random normal values,
# randomly located in the unit cube
n <- 50
 
data3D <- data.frame(x = runif(n), y = runif(n), z = runif(n), v = rnorm(n))
coordinates(data3D) = ~x+y+z
 
range1D <- seq(from = 0, to = 1, length = 20)
grid3D <- expand.grid(x = range1D, y = range1D, z = range1D)
gridded(grid3D) = ~x+y+z
 
res3D <- krige(formula = v ~ 1, data3D, grid3D, model = vgm(1, "Exp", .2))
 
library(lattice)
 
levelplot(var1.pred ~ x + y | z, as.data.frame(res3D))
 

Show variogram

show.vgms()
show.vgms(models = "Exp", sill=0.03, range=30, max=100)

Grid

xy <- expand.grid(x=1:100, y=1:100)
xy <- expand.grid(x=1:100, y=1:100, z=1:100)

Utility commands

Interpolation

n <- seq(0.001, 0.2, by=0.001)
k <- (234.0*(n)^1.25+2094.0*(10*n)^3.88)*1e-18
f_kn <- approxfun(k, n)

Generate random numbers from uniform distribution

runif(12, -1, 1)

Round numbers to two decimal place

round(sample, 2)

Version of R

R.version

Find out if the method is available in any package

help.search("lqs")

Show current directly

getwd()

Change current directly

setwd("E:\tmp\user")

Show list of files on the current directly

list.files()

remove all variables

rm(list = ls())

Define as a spatial object

coordinates(sim) <- ~x+y

Min and max values in data

range(mhw$grain)

which.min(grain)
mhw[which.min(grain), "straw"]

Subset

subset(mhw, grain < 3)