################################################################################
### R BASICS WORKSHOP                                                        ###
### PRESENTATION 8: GRAPHICS                                                 ### 
###                                                                          ###
### Center for Conservation and Sustainable Development                      ###
### Missouri Botanical Garden                                                ###
### Website: rbasicsworkshop.weebly.com                                      ### 
################################################################################

# The material below is largely based on chapter 4 of "R for beginners". It is
# divided in seven sections:

# A) Introduction
# B) Managing graphics and graphical devices: opening multiple graphical devices 
# C) Managing graphics and graphical devices: partitioning a graphical device
# D) High level graphical functions (create a graph) 
# E) Low level graphical functions (affect an existing graph)
# F) Graphical parameters (determine graph format and can be modified using 
#    function "par")
# G) Graphical functions whose results can be assigned to an object (in contrast 
#    to "plot" or "points")
 

################################################################################
# A) Introduction
################################################################################

# R offers a wide variety of graphs. Type this code to get a sense of 
# what is possible:
demo(graphics)

#from the help page for function "image"
# A prettier display of the volcano
x <- 10*(1:nrow(volcano))
y <- 10*(1:ncol(volcano))
image(x, y, volcano, col = terrain.colors(100), axes = FALSE)
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))
box()
title(main = "Maunga Whau Volcano", font.main = 4)

# At the end of this section we will see some graphical functions whose results 
# can be assigned to objects (functions "hist", "boxplot", "barplot"). However, the 
# result of several graphical functions cannot be assigned to an object, and is 
# instead sent to a graphical device. Graphical devices are windows or files.


################################################################################
# B) Managing graphics and graphical devices: opening multiple graphical devices
################################################################################

# When executing a graphical function R opens a window to show the graph if 
# there is no open graphical device:
data(iris)
plot(iris$Sepal.Length, iris$Sepal.Width, pch=19) 
# you can use the menu under "File" to save the figure

#Graphical windows are called with function "dev.new":
dev.new()
plot(iris$Sepal.Length, iris$Petal.Length, pch=19)

# Graphical devices that are files are called with a function according to 
# the file type of interest. Before executing the code below, make sure the 
# working directory is where you want to store graphical files. Use functions 
# "getwd" and "setwd" as needed.

jpeg(file="SepalLenght_vs_SepalWidth.jpeg")
plot(iris$Sepal.Length, iris$Sepal.Width, pch=19)
dev.off() #closes graphical device

png(file="SepalLenght_vs_SepalWidth.png")
plot(iris$Sepal.Length, iris$Sepal.Width, pch=19)
dev.off()

pdf(file="SepalLenght_vs_SepalWidth.pdf")
plot(iris$Sepal.Length, iris$Sepal.Width, pch=19)
dev.off()

postscript(file="SepalLenght_vs_SepalWidth.ps")   # often used for publication
plot(iris$Sepal.Length, iris$Sepal.Width, pch=19)
dev.off()

postscript(file="SepalLenght_vs_SepalWidth.eps")  # often used for publication
plot(iris$Sepal.Length, iris$Sepal.Width, pch=19)
dev.off()

# A list of available graphical devices is obtained with this code:
?device

# You can have several graphical devices open at once. To determine the number of
# open devices at any one time use:
dev.list()

# To determine which graphical device is active at any one time use:
dev.cur()

# To change the active graphical device:
dev.set(2)

# To close a graphical device:
dev.off(3)

# To close the active graphical device:
dev.off()


################################################################################
# C) Managing graphics and graphical devices: partitioning a graphical device
################################################################################

# Function "layout" divides the graphical device into parts where graphs will be 
# displayed successively. Its main argument is a matrix with integer numbers 
# indicating the numbers of the sub-windows.
layout(matrix(1:4, 2, 2)) #see help for "layout"
layout.show(4)		#see help for "layout.show"
plot(iris$Sepal.Length, iris$Sepal.Width, pch=19, cex.lab=1.5, cex.axis=1.5, 
    xlab="Sepal length (cm)", ylab="Sepal width (cm)")
plot(iris$Sepal.Length, iris$Petal.Length, pch=19, cex.lab=1.5, cex.axis=1.5, 
    xlab="Sepal length (cm)", ylab="Petal length (cm)")
plot(iris$Sepal.Length, iris$Petal.Width, pch=19, cex.lab=1.5, cex.axis=1.5, 
    xlab="Sepal length (cm)", ylab="Petal width (cm)") 
boxplot(iris$Sepal.Length ~ iris$Species, cex.lab=1.5, cex.axis=1.5,
	xlab="Species", ylab="Sepal Length (cm)",
	names=c(expression(italic("I. set.")), expression(italic("I. ver.")), expression(italic("I. vir."))))

layout(matrix(1:6, 3, 2))
layout.show(6)
plot(iris$Sepal.Length, iris$Sepal.Width, pch=19)
plot(iris$Sepal.Length, iris$Petal.Length, pch=19)
plot(iris$Sepal.Length, iris$Petal.Width, pch=19)
plot(iris$Sepal.Width, iris$Petal.Length, pch=19)
plot(iris$Sepal.Width, iris$Petal.Width, pch=19)
plot(iris$Petal.Length, iris$Petal.Width, pch=19)

# By default layout() divides the graphical device into regular heights and 
# widths, use arguments "widths" and "heights" to modify this division

layout(matrix(1:6, 3, 2))
layout.show(6)

layout(matrix(1:6, 2, 3))
layout.show(6)

m <- matrix(c(1:3, 3), 2, 2)
layout(m)
layout.show(3)

m <- matrix(1:4, 2, 2)
layout(m, widths=c(1, 3), heights=c(3, 1))
layout.show(4)

m <- matrix(c(1,1,2,1),2,2)
layout(m, widths=c(2, 1), heights=c(1, 2))
layout.show(2)

layout(matrix(1, 1, 1))
layout.show(1)


################################################################################
# D) High level graphical functions (create a graph)
################################################################################

#a few examples, but see pages 40-41 in "R for beginners"
iris[1:5,]
plot(iris$Sepal.Length, iris$Sepal.Width)
plot(iris$Petal.Length, iris$Petal.Width, xlab="Sepal length (cm)", 
    ylab="Petal Width (cm)", cex.axis=1.5, cex.lab=1.5, bty="n", pch=19)
sunflowerplot(iris$Sepal.Length, iris$Sepal.Width)
boxplot(iris$Sepal.Length ~ iris$Species)
boxplot(iris$Sepal.Length ~ iris$Species, names=expression(italic("Iris setosa"), 
    italic("Iris versicolor"), italic("Iris virginica")), 
    ylab="Sepal length (cm)", xlab="", cex.axis=1.5, cex.lab=1.5)
coplot(iris$Petal.Length ~ iris$Petal.Width | iris$Sepal.Length, overlap=0, pch=19)
pairs(iris)
hist(iris$Sepal.Length, breaks=25)


################################################################################
# E) Low level graphical functions (modify an existing graph)
################################################################################

# Below are a few examples with functions "points" and "legend"; see pages 41-42 
# in "R for beginners" for more low level graphical functions

plot(iris$Petal.Length, iris$Petal.Width, xlab="Sepal length (cm)", 
    ylab="Petal Width (cm)", cex=1.3, cex.axis=1.5, cex.lab=1.5, bty="n")
    
points(iris$Petal.Length[iris$Species=="setosa"], 
    iris$Petal.Width[iris$Species=="setosa"], cex=1.3, pch=19, col="red")
points(iris$Petal.Length[iris$Species=="versicolor"], 
    iris$Petal.Width[iris$Species=="versicolor"], cex=1.3, pch=19, col="blue")
points(iris$Petal.Length[iris$Species=="virginica"], 
    iris$Petal.Width[iris$Species=="virginica"], cex=1.3, pch=19, col="green")
legend("topleft", c("Iris setosa", "I. versicolor", "I. virginica"), 
    pch=19, col=c("red", "blue", "green"), cex=1.3)
legend("bottomright", c(expression(italic("Iris setosa")), 
    expression(italic("Iris versicolor")), expression(italic("Iris virginica"))), 
    pch=19, col=c("red", "blue", "green"), cex=1.3)
legend(1.3, 1.9, c(expression(italic("Iris setosa")), 
    expression(italic("Iris versicolor")), expression(italic("Iris virginica"))), 
    pch=19, col=c("red", "blue", "green"), cex=1.3)


################################################################################
# F) Graphical parameters (can be modified using function "par")
################################################################################

# You will take full control of graphs by getting to know function "par" and the 
# associated help page. The following three blocks of code illustrate the use of 
# function "par" as well as three low level graphical functions: "axis", "mtext" 
# and "legend":

par(mar = c(5, 4, 4, 2) + 0.1) # These are default values, see help for "par"
hist(iris$Sepal.Length, breaks=seq(0,8,0.5), col="red2", border="red4", 
    xlab="", ylab="", main="", xaxt="n", yaxt="n", ylim=c(0,40))
par(new=T)
hist(iris$Petal.Length, breaks=seq(0,8,0.5), density=30, xlab="", ylab="", 
    main="", xaxt="n", yaxt="n", ylim=c(0,40), col="blue")
axis(1, at=seq(0,8,1), cex.axis=2, lwd=2)
mtext("Length (cm)", side=1, line=3, cex=2)
axis(2, cex.axis=1.5, cex.axis=2, lwd=2)
mtext("Specimens", side=2, line=3, cex=2)
legend("topright", c("Sepals", "Petals"), fill=c("red2", "blue"), 
    density=c(1000, 30), pt.cex=2, border=c("red4", "blue"), cex=2)
# There are problems, see them? In the next block of code the problems
# are solved using argument "mar" of function "par" and argument "ylim"
# of function "hist".

par(mar = c(5, 5, 4, 2) + 0.1)
hist(iris$Sepal.Length, breaks=seq(0,8,0.5), col="red2", border="red4", 
    xlab="", ylab="", main="", xaxt="n", yaxt="n", ylim=c(0,42))
par(new=T)
hist(iris$Petal.Length, breaks=seq(0,8,0.5), density=30, xlab="", ylab="", 
    main="", xaxt="n", yaxt="n", ylim=c(0,42), col="blue")
axis(1, at=seq(0,8,1), cex.axis=2, lwd=2)
mtext("Length (cm)", side=1, line=3, cex=2)
axis(2, cex.axis=1.5, cex.axis=2, lwd=2)
mtext("Specimens", side=2, line=3, cex=2)
legend("topright", c("Sepals", "Petals"), fill=c("red2", "blue"), 
    density=c(1000, 30), pt.cex=2, border=c("red4", "blue"), cex=2)

par(mfcol = c(1,2))
hist(iris$Sepal.Length, breaks=seq(0,8,0.5), col="gray", border="black", 
    ylim=c(0,40), xlab="Length (cm)", ylab="Specimens", cex.lab=1.5, 
    cex.axis=1.5, main="Sepals", cex.main=1.5)
hist(iris$Petal.Length, breaks=seq(0,8,0.5), col="gray", border="black", 
    ylim=c(0,40), xlab="Length (cm)", ylab="Specimens", cex.lab=1.5, 
    cex.axis=1.5, main="Petals", cex.main=1.5)


################################################################################
# G) Graphical functions whose results can be assigned to an object (in contrast
# to "plot" or points)
################################################################################

par(mfcol = c(1,1))
h.Sepal.Length <- hist(iris$Sepal.Length, breaks=seq(0,8,0.5))
h.Sepal.Length
class(h.Sepal.Length)
attributes(h.Sepal.Length)
str(h.Sepal.Length)
h.Sepal.Length$mids
h.Sepal.Length$counts

bp.Sepal.Length <- boxplot(iris$Sepal.Length ~ iris$Species)
bp.Sepal.Length
class(bp.Sepal.Length)
attributes(bp.Sepal.Length)
str(bp.Sepal.Length)
bp.Sepal.Length[[1]]
bp.Sepal.Length[[1]][1,]
bp.Sepal.Length[[2]]