RSA_geoarch.r

#get data from google sheet
# connect to google sheet
require(RCurl)
options(RCurlOptions = list(capath = system.file("CurlSSL", "cacert.pem", package = "RCurl"), ssl.verifypeer = FALSE))
#in google spreadsheet, go to file-> publish to web -> get link to publish to web -> get csv file
goog <- "https://docs.google.com/spreadsheet/pub?key=0As7CmPqGXTzldFRsVi1VZ2EyNXJ1ZEV5SG5GSExwRHc&single=true&gid=5&output=csv"
data <- read.csv(textConnection(getURL(goog)), stringsAsFactors = FALSE)

# extract just data for plotting: pH, SOM, CaCO3, MS-LF, MS-FD
plotting_data <- na.omit(data[,c('Sample.number',
                                 'Average.pH', 
                                 'mean.EC',
                                 'LOI.percent.organic.material..Average', 
                                 'LOI.Percent.carbonate.content..Average',
                                 'MS.LF.reading.average.mass.susceptibility',
                                 'MS.FD'
)])
# rename the columns of the plotting data with shorter names and measurement units
names(plotting_data) <- c("depth (m)", "pH", "EC", "Soil organic \nmatter \n(% mass)", "Carbonates \n(% mass)", "Magnetic \nsusceptibility \n(SI units)", "Frequency \ndependency \n(%)")
# replace missing data in google sheet with 0 so it can be handled by the plotting function
plotting_data[plotting_data == '#DIV/0!']  <- 0      
# convert all the numbers to numbers (some columns are in character type)
plotting_data <- as.data.frame(sapply(plotting_data, as.numeric, as.character))
# create vector for depth of samples
depth <- na.omit(as.numeric(plotting_data$`depth (m)`))
depth <- 1:20
# create data frame of eveything except depth
plotting_data <- plotting_data[,-(names(plotting_data) == 'depth (m)')]

# setup for plotting
install.packages("rioja")
library(rioja)

#### plain plot 

x <- strat.plot(plotting_data,
                yvar = depth, 
                y.rev = TRUE, 
                ylabel = "Depth below surface (m)", 
                col.line = "black", # try others that you like
                lwd.line = 1       # ditto
)
dev.off() # clear that plot

#### biplots: for plotting two of the variables to explore their relationship
library(scales)
library(ggplot2)
library(stringr)
# change what's in the "": lookback up to line 25 for the proper names
x <- as.numeric(as.character(plotting_data [,which(colnames(plotting_data)=="Magnetic \nsusceptibility \n(SI units)")]))
y <- as.numeric(as.character(plotting_data [,which(colnames(plotting_data )=="Soil organic \nmatter \n(% mass)")]))
#
#
# change the xlab and ylab below
tmp <- data.frame(cbind(x,y))
ggplot(tmp, aes(x, y, label=depth)) +                            # make the basic plot object
  scale_x_continuous(limit=c(min(tmp$x),max(tmp$x)), 
                     breaks=round(fivenum(tmp$x),1)) +           # set the locations of the x-axis labels
  # change x label to match the data you've selected:
  xlab(paste("Magnetic sus (r = ", round(cor(x,y),4),", p-value = ", 
             round(anova(lm(y~x))$'Pr(>F)'[1],4),")",sep="")) +  # paste in the r and p values with the axis label
  scale_y_continuous(limit=c(min(tmp$y),max(tmp$y)),
                     breaks=round(fivenum(tmp$y),1)) +           # set the locations of the y-axis labels
  # change y label to match the data you've selected:
  ylab("Soil Organic Matter (% mass)") +                   
  geom_text() +                 # specify that the data points are the sample names
  geom_rug(size=0.1) +         # specify that we want the rug plot
  geom_smooth(method="lm", se=FALSE, colour="grey80", aes(group=1)) + 
  # specify a line of best fit calculated by a linear model
  # se is the standard error area, change to TRUE to display it
  theme(panel.background = element_blank(),          # suppress default background
        panel.grid.major = element_blank(),              # suppress default major gridlines
        panel.grid.minor = element_blank(),              # suppress default minor gridlines
        axis.ticks = element_blank(),                    # suppress tick marks
        axis.title.x=element_text(size=12),              # increase axis title size slightly  
        axis.title.y=element_text(size=12, angle=90),    # increase axis title size slightly and rotate
        axis.text.x=element_text(size=12),               # increase size of numbers on x-axis
        axis.text.y=element_text(size=12),               # increase size of numbers on y-axis
        aspect.ratio=1)                                  # force the plot to be square

#get correlations of all variables
cors <- cor(plotting_data)