############################################################
#                                                          #
#   R Code for demonstrations in course                    #
#                                                          #
#   "Using R for statistical data analysis and graphics"   #
#                                                          #
#   Section 4: Graphics                                    #
#                                                          #
############################################################

# 2009-10-12 A. Papritz

library( sp )

setwd( "~/Desktop" )

## 
 # load meuse data set and get some information about the data
 ##

data( meuse )
data( meuse.riv )

?meuse
?meuse.riv

str( meuse )
summary( meuse )



## 
 # Scatterplots
 ##

# of two vectors

plot( x = meuse[, "x"], y = meuse[, "y"], 
    xlab = "easting", ylab = "northing",
    main = "position of soil sampling locations"
)

# of a single vector

plot( x = meuse[,"zinc"], ylab = "zinc")

# of a two column ``matrix''

# plot( meuse[, c("x","y")], main = "sampling locations" )

# using a formula

plot( 
    zinc ~ dist.m, data = meuse, 
    main = "Zn vs. distance"
)


# Common arguments

# controlling the scales of axes

plot( 
    zinc ~ dist.m, data = meuse, 
    log = "y", 
)

plot( 
    zinc ~ dist.m, data = meuse, 
    xlim = c( 5, 2000 ), ylim = c( 100, 5000 )
)

plot( 
    y ~ x, data = meuse,
    asp = 1,
    main = "sampling locations" 
)


# controlling type, color and size of symbols

plot( 
    zinc ~ dist.m, data = meuse, log = "y",
    pch = 3, col = "red", cex = 3 
)

plot( 1:25, pch = 1:25, cex = 2 )


# "bubble plot": a more refined scatterplot

plot( 
    y ~ x, data = meuse, asp = 1, 
    col = as.numeric( meuse[, "ffreq"] ),
    cex = sqrt( meuse[, "zinc"] )/10
)



## 
 # Boxplots
 ##

# of a single variable

boxplot( meuse[, "zinc"], ylab = "zinc" )

# horizontal with notches and a logarithmic scale

boxplot( 
    meuse[, "zinc"], horizontal = TRUE, notch = TRUE,
    log = "x", xlab = "zinc [mg/kg]" 
)

# of several variables

boxplot( 
    meuse[, c("zinc", "lead", "copper", "cadmium")], 
    log = "y", 
    xlab = "metal content [mg/kg]"
)

# of values of a single variable for several groups

boxplot(
    zinc ~ ffreq, data = meuse,
    log = "y", notch = TRUE,
    names = c( "often", "intermediate", "rarely" ),
    xlab = "flooding frequency", ylab = "zinc [mg/kg]"
)




## 
 # Adding points and lines
 ##

# adding points

plot( 
    lead ~ dist.m, data = meuse, 
    ylim = c( 10, 1000 ),
    pch = 1 , log = "y"
)

points(
    x = meuse[, "dist.m"], y = meuse[, "copper"], 
    col = c( "red", "blue" ), pch = 2:18, cex = seq( 0.1, 2, by = 0.1 )
)

# adding lines

str( meuse.riv )

plot( y ~ x, meuse, asp = 1, pch = 16, cex = 1.5 )

lines( 
    x = meuse.riv[,1], y = meuse.riv[,2],
    lty = "solid", lwd = 2, col = "cyan" 
)

# keeping control of the sequence in which the different things are plotted

plot( y ~ x, meuse, asp = 1, type = "n" )

lines( 
    x = meuse.riv[,1], y = meuse.riv[,2], 
    lty = "solid", lwd = 2, col = "cyan" 
)

points( x = meuse[, "x"], y = meuse[, "y"], pch = 16, cex = 1.5 )

# adding straight lines

plot( lead ~ dist.m, meuse )

abline( h = c( 200, 500 ), col = c( "orange", "red"), lty = "dotted" ) 

abline( v = 350, col = 4, lty = 3, lwd = 2 ) 

abline( a = 500, b = -1, lty = "dashed", lwd = 2, col = "black" ) 

# adding line segments

plot( y ~ x, meuse, asp = 1, pch = 16, cex = 1.5 )

t.xr <- range( meuse[, "x"] )
t.yr <- range( meuse[, "y"] )

t.xr
t.yr

cbind(rep(t.xr[1],2), t.yr )
cbind( rep(t.xr[2],2), t.yr[2:1] )

segments( rep(t.xr[1],2), t.yr, rep(t.xr[2],2), t.yr[2:1], col = "red", lwd = 5 )

# adding polygons

plot( y ~ x, meuse, asp = 1, pch = 16, cex = 1.5 )

polygon( 
    x = meuse.riv[,1], y = meuse.riv[,2], 
    lwd = 2, border = "blue", col = "cyan", 
    angle = 135, density = 10 
)





## 
 # Adding text
 ##

# adding labels to points

plot( y ~ x, meuse, asp = 1, pch = 16, cex = 1.5 )
text( 
    x = meuse[,"x"], y = meuse[,"y"], 
    labels = meuse[, "landuse"], pos = 4 
)

# or 

plot( y ~ x, meuse, asp = 1, type = "n" )
text( meuse[, c("x","y")], labels = meuse[, "landuse"] )

# adding text within the window of the plot

plot( y ~ x, meuse, asp = 1, pch = 16, cex = 1.5 )

t.xr <- range( meuse[, "x"] )
t.yr <- range( meuse[, "y"] )

segments( rep(t.xr[1],2), t.yr, rep(t.xr[2],2), t.yr[2:1], col = "red", lwd = 5 )

text( 177800, 331500, cex = 2, 
    labels = "this is non-sense", pos = 4, col = "red", srt = 90
)

# adding a legend

plot( 
    y ~ x, meuse, asp = 1, 
    col = as.numeric( meuse[, "ffreq"] ),
    cex = sqrt( meuse[, "zinc"] )/10
)

legend( 
    x = "topleft", 
    legend = c( "flooding frequency", "often", "intermediate", "rarely" ),
    pch = c( NA, rep( 1, 3 ) ), col = c( NA, "black", "red", "green" )
)

t.zinc <- c( 100, 200, 500, 1000, 2000 )
legend( 
    x = 182200, y = 329600, xjust = 1, yjust = 0,
    x.intersp = 1.5, y.intersp = 1.7,
    legend = c( "zinc [mg/kg]", as.character( t.zinc ) ),
    pch = 1, pt.cex = c( NA, sqrt( t.zinc )/10 ),
    bty = "n"
)




## 
 # Interacting with graphical displays
 ##

# Identify points in a plot

plot( 
    y ~ x, meuse, asp = 1, 
    col = as.numeric( meuse[, "ffreq"] ),
    cex = sqrt( meuse[, "lead"] )/10
)

identify( x = meuse[,"x"], y = meuse[,"y"] )

identify( 
    x = meuse[,"x"], y = meuse[,"y"], 
    labels = meuse[,"lead"], col = "blue" 
)

identify(
    meuse[, c("x","y")], 
    labels = meuse[,"landuse"]
)

# read off the coordinates of points in a plot

plot( 
    y ~ x, meuse, asp = 1, 
    col = as.numeric( meuse[, "ffreq"] ),
    cex = sqrt( meuse[, "lead"] )/10
)

locator()

locator( 3, type = "p", pch = 3, col = "blue"  )

plot( 
    y ~ x, meuse, asp = 1, 
    col = as.numeric( meuse[, "ffreq"] ),
    cex = sqrt( meuse[, "lead"] )/10
)

t.polygon <- locator(, type = "p", pch = 3, col = "blue")
t.polygon

polygon( t.polygon , border ="cyan" )




## 
 # Controlling graphics devices
 ##

dev.new( 
    height = 10, width = 10, 
    pointsize = 14 
)
pairs( meuse[, c("zinc", "lead", "copper", "cadmium")] )

# generating a pdf file

pdf( file = "plot.pdf", height = 10, width = 10, pointsize = 14 )

plot( zinc ~ dist.m, data = meuse )

# list open devices

dev.list()

# list active device

dev.cur()

dev.off()

dev.list()
dev.cur()


graphics.off()

dev.list()





## 
 # Querying and setting default parameters to control appearance of a
 # graphic
 ##


# querying and changing default parameters

par( "pch", "lty", "cex" )

par()

# changing default parameters

par( pch = 3, lty = "dotted", col = "red" )

plot( y ~ x, data = meuse )
lines( x = meuse.riv[,1], y = meuse.riv[,2] )

# storing old values of parameters

t.oldpar <- par( pch = 5, col = "magenta", lty = "dotdash" )

t.oldpar

par( "pch", "col", "lty" )

plot( y ~ x, data = meuse )
lines( x = meuse.riv[,1], y = meuse.riv[,2] )

# restoring old parameters values

par( t.oldpar )

plot( y ~ x, data = meuse )
lines( x = meuse.riv[,1], y = meuse.riv[,2] )

dev.off()


## 
 # Multiple plots in one graphic
 ##

# by row 

par( mfrow = c( 2, 2 ) )

plot( 
    y ~ x, data = meuse, 
    cex = sqrt( meuse[,"zinc"])/10, asp = 1 
)

plot( zinc ~ dist.m, data = meuse )
hist( meuse[, "zinc"] )

boxplot( 
    meuse[, c("lead", "copper", "cadmium")], 
    log = "y", lab = "metal content [mg/kg]"
)
  
# by column

dev.new()
par( mfcol = c( 2, 2 ) )

plot( 
    y ~ x, data = meuse, 
    cex = sqrt( meuse[,"zinc"])/10, asp = 1 
)

plot( zinc ~ dist.m, data = meuse )
hist( meuse[, "zinc"] )

boxplot( 
    meuse[, c("lead", "copper", "cadmium")], 
    log = "y", lab = "metal content [mg/kg]"
)

# par( "cex" )

graphics.off()


# controlling the figure region

par( fig = c( 0, 0.5, 0, 0.7 ) )
plot( 1:10 )

par( fig = c( 0.5, 1, 0.2, 1 ) )
plot( 1:20, col = "red" )


par( fig = c( 0, 0.5, 0, 0.7 ) )
plot( 1:10 )

par( 
    fig = c( 0.5, 1, 0.2, 1 ), 
    new = TRUE 
)
plot( 1:20, col = "red" )

# controlling the plot regions

dev.new()
par( 
    fig = c( 0, 0.5, 0, 0.7 ),
    plt = c( 0, 1, 0.2, 1 )
)
plot( 1:10 )

par( 
    fig = c( 0.5, 1, 0.2, 1 ), 
    plt = c( 0, 1, 0, 1 ),
    new = TRUE 
)
plot( 1:20, col = "red" )

graphics.off()


# setting the width of the outer margins

dev.new()

plot( 1:10 )
par( "oma")

dev.new()
par( oma = c( 5, 5, 5, 5 ) )
plot( 1:10 )

dev.off()

# setting the width of the plot margins

dev.new()

par( mar = c( 10, 10, 0, 0 ) )
plot( 1:10 )

graphics.off()



## 
 # Fine-tuning the plotting of axes
 ##

# adding customized axes

plot( zinc~dist.m, data = meuse, log = "y" )

dev.new()

plot( 
    log10( zinc)~dist.m, data = meuse,
    yaxt = "n",
    ylab = "zinc [mg/kg]"
)

axis( 3, at = seq( 100, 1000, by = 100 ) )

t.yticks <- c( (1:9)* 100, (1:9) * 1000 )
t.yticks 

axis( 2, at = log10( t.yticks ), labels = t.yticks )

graphics.off()

# more control on axis labelling

dev.new()
plot( zinc~dist.m, data = meuse, log = "y" )

dev.new()

par( las = 2 )

plot( zinc~dist.m, data = meuse, log = "y" )

graphics.off()




## 
 # More on color and size
 ##

# setting color of text annotations

par( 
    col.main = "magenta", cex.main = 3,
    col.axis = "red", cex.axis = 0.7,
    col.lab = "green", cex.lab = 1.5
)

plot( zinc~dist.m, main = "many colors", data = meuse)

# superseding on the fly the defaults set previously by par

plot( zinc~dist.m, main = "many colors", data = meuse, col.main = "blue")


# settting the background color

par( bg = "darkblue" )

plot( zinc~dist.m, main = "many colors", data = meuse, log = "y" )

# ... and foreground color

par( fg = "yellow" )

plot( zinc~dist.m, main = "many colors", data = meuse, log = "y" )

graphics.off()


# Managing color scales

# built-in colors

colors()

# default color scale

palette()

plot( 1:16, col = 1:16, pch = 16, cex = 3)

# rainbow color scale

palette( rainbow( 10 ) )
palette()

plot( 
    zinc ~ dist.m, meuse, pch = 16, cex = 3,
    col= as.numeric( cut(meuse[, "zinc"], 10) ) 
)


# other built-in color scales

palette( terrain.colors( 16 ) )
plot( 1:16, col = 1:16, pch = 16, cex = 3 )

palette( heat.colors( 16 ) )
plot( 1:16, col = 1:16, pch = 16, cex = 3 )

# restoring the default color scale

palette("default" )
palette()



## 
 # Further high-level plotting functions
 ##


# barplot

# of a vector

barplot( 1:10 )

# of a matrix (stacked)

barplot( VADeaths, ylab = "deaths per 1000 persons" )

# of a matrix (non-stacked)

barplot( 
    VADeaths, beside = TRUE,
    ylab = "deaths per 1000 persons",
    legend.text = rownames( VADeaths )
)



# scatterplot matrix

pairs( 
    meuse[, c("lead", "zinc", "copper", "cadmium")], 
    upper.panel = NULL 
)
  

# multiple scatterplots in one plot

matplot(
    x = meuse[, "dist.m"], 
    y = meuse[, c("lead", "zinc", "copper", "cadmium")],
    pch = 1:4, col = 1:4, log = "y",
    xlab = "distance to river", 
    ylab = "metal content [mg/kg]"
)

legend( 
    "topright", pch = 1:4, col = 1:6,
    legend = c("lead", "zinc", "copper", "cadmium"),
    horiz = TRUE
)


# 3D-scatterplot

library( scatterplot3d )

scatterplot3d(
    x = meuse[,"x"], y = meuse[,"y"],
    z = log10( meuse[,"zinc"]), type = "h",
    highlight.3d = TRUE
)



# image and contour

x <- 10*(1:nrow(volcano))
y <- 10*(1:ncol(volcano))

?volcano

image(
    x = x, y = y, z = volcano,
    col = terrain.colors( 100 ), 
    xaxt = "n", yaxt = "n",
    main = "Maunga Whau Volcano"
)

# add contour lines

contour(
    x, y, volcano, 
    levels = seq(90, 200, by = 5),
    add = TRUE, col = "peru"
)
    
axis( 1, at = seq( 100, 800, by = 100 ) )
axis( 2, at = seq( 100, 600, by = 100 ) )



# persp

z <- 2 * volcano
x <- 10 * (1:nrow(z))
y <- 10 * (1:ncol(z))

par( bg = "slategray" )

persp(
    x, y, z, 
    theta = 135, phi = 30, 
    col = "green3", scale = FALSE,
    ltheta = -120, 
    shade = 0.75, border = NA, box = FALSE
)







#############################
#                           #
#   R Code used on slides   #
#                           #
#############################

# overview

d.sport <- read.table(
    "http://stat.ethz.ch/Teaching/Datasets/WBL/sport.dat", header=TRUE
)
data( sleep )
meuse <- read.table(
    "ftp://uwis-ftp.ethz.ch/ITES/papritz/UsingR/data/meuse.txt", header=TRUE
)
meuse.riv <- read.table(
    "ftp://uwis-ftp.ethz.ch/ITES/papritz/UsingR/data/meuse_riv.txt", header=TRUE
)

t.y1 <- sleep[sleep[,"group"]==1,"extra"]
t.y2 <- sleep[sleep[,"group"]==2,"extra"]

plot(d.sport[,"kugel"],~d.sport[,"speer"],xlab="ball~push",ylab="javelin",pch=7)
plot(punkte~kugel+speer, data=d.sport)
pairs(d.sport)
boxplot(t.y1,t.y2,ylab="extra")

# scatterplot & boxplot

plot(x=meuse[,"x"], y=meuse[,"y"],  xlab="easting", ylab="northing", main="sampling locations")
plot(meuse[,"zinc"], ylab = "zinc")
plot(y~x, data=meuse, asp=1)
plot(zinc~dist, data=meuse)

boxplot(x=meuse[,"zinc"], notch=TRUE, horizontal=TRUE, log="x", xlab="zinc content")
boxplot(meuse[,c("zinc","lead")], horizontal=TRUE, log="x")
boxplot(zinc~ffreq, data=meuse)

# low-level functions
    
plot(lead~dist.m, data=meuse, ylim=c(10,1000), log="y")
points(meuse[,c("dist.m","copper")], col="red", pch=3)

plot(y~x, meuse, asp=1)
lines(x=meuse.riv[,1], y=meuse.riv[,2], lty="dotted", lwd=2.5, col= "cyan")

plot(lead~dist.m, meuse, asp=1)
abline(h=c(200,500), lty="dotted", col=c("orange","red"))
abline(a=300, b=-1, lty=3, lwd=3)

plot(y~x, data=meuse)
t.xr<-range(meuse[,"x"])
t.yr<-range(meuse[,"y"])
segments(x0=rep(t.xr[1],2), y0=t.yr, x1=rep(t.xr[2],2), y1=t.yr[2:1] )

plot(y~x, meuse, asp=1)
polygon(x=meuse.riv[,1], y=meuse.riv[,2], border="blue", col="cyan", angle=135, density=10)


plot(y~x, data=meuse, asp=1, type="n")
text(x=meuse[,"x"], y=meuse[,"y"], cex=0.7, labels=meuse[,"landuse"])

plot(y~x, data=meuse, asp=1,
    col=as.numeric(meuse[,"ffreq"]),
    cex=sqrt(meuse[,"zinc"])/10)
legend(x="topleft", pch=c(NA,rep(1,3)),
  col=c("black","black","red","green"),
  legend=c("flooding", "often", "intermediate","rarely"))


plot(y~x, data=meuse, asp=1,
    cex=sqrt(meuse[,"zinc"])/10))
identify(x=meuse[,"x"], y=meuse[,"y"], labels=meuse[,"zinc"])


plot(y~x, data=meuse, asp=1)
polygon(locator())

# par

par(mfrow=c(1,2))
plot(lead~dist.m, data=meuse)
plot(zinc~dist.m, data=meuse )










# # Variants of scatterplots
# 
# # add a smooth curve
# 
# scatter.smooth( 
#     x = meuse[, "dist.m"], y = meuse[, "zinc"], 
#     span = 0.75, 
#     xlab = "distance", ylab = "zinc",
#     main = "Zn vs. distance"
# )
# 
# # scatterplots of several y- against one x-variable
# 
# matplot( 
#     x = meuse[, "dist.m"], y = meuse[, c("zinc", "lead", "copper")],
#     log = "y", pch = 1:3,
#     ylab = "zinc, lead, copper", xlab = "distance"
# )
# 
# # conditioning plots
# 
# coplot(
#     zinc ~ dist.m | ffreq, 
#     data = meuse, 
#     rows = 1, columns = 3, 
#     col = as.numeric( meuse[, "soil"] ),
#     pch = as.numeric( meuse[, "soil"] )
# )
# 
#