TOPIC: POTM-02 (June): Multiple y-axes

POTM-02 (June): Multiple y-axes 1 year 2 weeks ago #22

Kim_Cressman
Offline
Moderator
Posts: 19
Karma: 2

Hello everyone! Here's the June Plot of the Month.

This is a plot that shows three different parameters, each on a different scale. I used it to look at data over a tidal cycle – phosphate from our ISCO sampler, plotted alongside salinity and depth from our YSI. You can see each of these parameters has its own y-axis on the left. Please forgive where the legend is; I haven't figured out all the subtleties of legend placement. But hopefully the point I've gotten to will at least be helpful for others!

I'd also love to see how this type of plot works out with other people's data, if anyone wants to have a go at it. See how chlorophyll relates to pH and DO over a tidal cycle; see how NH4, PO4 and chl all look on the same graph; whatever you're interested in. You can probably even play with the margins and make more than 3 axes.

Step 1: Download some data. The CDMO’s Advanced Query System has a Merge Query function that will package water quality, weather, and nutrient data together. I didn't download any weather data, but you could definitely do something with it. Nutrient data will show up in the row of WQ/MET data that’s closest to your sample time. So a nutrient sample taken on 3/14 at 15:16 is going to be on the same line as your WQ readings at 3/14 15:15. The WQ and NUT samples each get their own DateTimeStamp column.

Slight tangent about the Merge Query function and column names: You can combine many different sites into one file through the Merge Query function. The column headers aren’t as simple as they are from the CDMO’s other export tools. Instead of PO4F, for example, it comes back as GNDBCNUT_PO4F. Useful when you want data from multiple sites, right? I could have gotten GNDBLNUT_PO4F and GNDPCNUT_PO4F in the same file, and still been able to tell what was what.

To maintain sanity, I only worked with one site in each data file, and simplified column names by going into the data file and turning GNDBCNUT_PO4F into PO4. GNDBCWQ_Depth became Depth. GNDBCWQ_Sal became just Sal. You see where I’m going with this. I wanted to make a script that I could apply to multiple different files as long as I formatted the files the same way. But you could certainly handle this differently.

Step 2: Filter data in your file to get rid of what you don't need. Since I was playing with nutrient data, I deleted the vast majority of YSI readings, because they didn't have associated nutrient samples. I also ended up deleting rows with grab samples because the timestamps made my x-axis messy – so I got rid of everything that didn’t correspond to an ISCO nutrient data point. Then I saved it as a csv.

Step 3: Play in R! I’ve heavily annotated the script I used, because I’m still a beginner and I like to know why I need to do certain things. This is all based on a blog post I found at: http://www.r-bloggers.com/multiple-y-axis-in-a-r-plot/

Now on to the script!

#set working directory
#this is where R will look for files
#it's also where R will save files you create in this session
setwd('C:/Users/kimberly.cressman/Desktop/Modified after emergency backup/Phosphate Publication/data')
#read in the data file
#don't forget R is case sensitive
dat <- read.csv('BC_WQandNUT_2012.csv')
#make the date/time column POSIXct, so it will be treated
#like an actual date/time and not a factor
dat$DateTime <- as.POSIXct(dat$DateTimeStamp, format = "%m/%d/%Y %H:%M")
#subset into months of interest
Sept.dat <- subset(dat, dat$DateTime >= '2021-09-01' & dat$DateTime <= '2021-09-30')
Oct.dat <- subset(dat, dat$DateTime >= '2021-10-01' & dat$DateTime <= '2021-10-31')
Nov.dat <- subset(dat, dat$DateTime >= '2021-11-01' & dat$DateTime <= '2021-11-30')
#start a pdf file for the graphs to go into
pdf(file='BCISCOplots.pdf')
#make wide margins for multiple y axes
#also increase the bottom margin so the x-axis labels can be angled
par(mar=c(7, 12, 4, 4) + 0.1)

That's all set-up so far. Now, I wanted to be able to copy and paste the graph bits for different months without having to go find everywhere I had put 'Sept.dat' and change it to 'Oct.dat' or whatever other month I wanted. So I made 'graph.dat' the main dataset used for the plots, and I assign different subsets of the data to it for whichever month I'm interested in at the time. This example uses my September data, which is the plot shown above.

###
###
#CHANGE ME
graph.dat <- Sept.dat
###
###
#plot PO4
#we are NOT plotting axes with the normal plot command; nor are we including axis labels
#this ylim command inside the plot command tells the plot to be as tall as the max of what's in your data
#missing data points screw up that 'max' function, so tell R to remove them with na.rm=T
plot(PO4~DateTime, data=graph.dat, axes=F, ylim=c(0, max(graph.dat$PO4, na.rm = T)), 
     xlab="", ylab="", type="l", col="black", main="")
axis(2, at=c(0,1000),lwd=2, labels = F)
axis(2, ylim=c(0,max(graph.dat$PO4, na.rm = T)), col="black",lwd=2)
mtext(2,text="PO4 conc (mg/L)",line=2)
#we just did two lines calling the y-axis because R will automatically
#make tick marks at certain intervals, and R wants to stop the axis at a tick mark
#but we want the line to extend above it, so we're making a second axis
#without tick marks. We're telling it to go from 0-1000, but it will
#really stop at the max of your data (because that's where the graph stops)
#plot depth
par(new=T)
plot(Depth~DateTime, data=graph.dat, axes=F, ylim=c(0,max(graph.dat$Depth, na.rm = T)), xlab="", ylab="", 
     type="l",lty=2, main="",lwd=2)
axis(2, at=c(0,1000),line = 3.5, lwd = 2, labels = F)
axis(2, ylim=c(0,max(graph.dat$Depth, na.rm = T)),lwd=2,line=3.5)
mtext(2,text="Sonde Depth (m)",line=5.5)
#plot salinity
par(new=T)
plot(Sal~DateTime, data=graph.dat, axes=F, ylim=c(0,max(graph.dat$Sal, na.rm = T)), xlab="", ylab="", 
     type="l",lty=3, main="",lwd=2)
axis(2, at=c(0,1000), lwd=2, line = 7, labels = F)
axis(2, ylim=c(0,max(graph.dat$Sal, na.rm = T)),lwd=2,line=7)
mtext(2,text="Salinity (psu)",line=9)
#plot x axis
#make and format a text string of date/time in the dataset
#because you can't just rotate text by 45 degrees through any of the default graphing commands
axisdate <- seq(min(graph.dat$DateTime), max(graph.dat$DateTime), by = '2 hours')
axisdate <- format(axisdate, '%m/%d %H:%M')
#creating the axis twice again, like we did for the y-axes above
#one will have tick marks and the other will just be a line that spans the entire graph
#NEITHER of these have labels because we have to use the text string we created above
axis.POSIXct(1,range(graph.dat$DateTime), labels=F)
axis(1, at=c(0,50000000000), labels=F)
#add in that text string below the axis
text(x=(graph.dat$DateTime), par("usr")[3] - 1, srt = 45, adj = 1,labels = axisdate, xpd = TRUE)
#plot the legend
legend(x = 'topright', inset = -0.05, legend=c("PO4","Depth","Sal"),lty=c(1,2,3), xpd = TRUE)
#close the pdf file
dev.off()

I was interested in multiple months, so I just copied and pasted that chunk of code and changed graph.dat at the beginning. (Any advice on looping? I'm sure there's a better way to do this than copy-and-paste, but I don't know it.)

You can spit all the graphs out into the same file; just save that dev.off() line for the very end, after you've made all the graphs. I'll insert a file with a few graphs next, then paste that whole script so you can all see what I mean.

File Attachment:

File Name: BCISCOplots2012.pdf
File Size:7 KB

#Multiple y-axis plots for BC ISCO in 2012
#2015-06-11 kac
#following example from: http://www.r-bloggers.com/multiple-y-axis-in-a-r-plot/
#set working directory
#this is where R will look for files
#it's also where R will save files you create in this session
setwd('C:/Users/kimberly.cressman/Desktop/Modified after emergency backup/Phosphate Publication/data')
#read in the data file
#don't forget R is case sensitive
dat <- read.csv('BC_WQandNUT_2012.csv')
#make the date/time column POSIXct, so it will be treated
#like an actualy date/time and not a factor
dat$DateTime <- as.POSIXct(dat$DateTimeStamp, format = "%m/%d/%Y %H:%M")
#subset into months of interest
Sept.dat <- subset(dat, dat$DateTime >= '2021-09-01' & dat$DateTime <= '2021-09-30')
Oct.dat <- subset(dat, dat$DateTime >= '2021-10-01' & dat$DateTime <= '2021-10-31')
Nov.dat <- subset(dat, dat$DateTime >= '2021-11-01' & dat$DateTime <= '2021-11-30')
#start a pdf file for the graphs to go into
pdf(file='BCISCOplots2012.pdf')
#make wide margins for multiple y axes
#also increase the bottom margin so the x-axis labels can be angled
par(mar=c(7, 12, 4, 4) + 0.1)
###--------------START GRAPH-----------------
###
###
#CHANGE ME
graph.dat <- Sept.dat
###
###
#plot PO4
plot(PO4~DateTime, data=graph.dat, axes=F, ylim=c(0, max(graph.dat$PO4, na.rm = T)), 
     xlab="", ylab="", type="l", col="black", main="")
axis(2, at=c(0,1000),lwd=2, labels = F)
axis(2, ylim=c(0,max(graph.dat$PO4, na.rm = T)), col="black",lwd=2)
mtext(2,text="PO4 conc (mg/L)",line=2)
#plot depth
par(new=T)
plot(Depth~DateTime, data=graph.dat, axes=F, ylim=c(0,max(graph.dat$Depth, na.rm = T)), xlab="", ylab="", 
     type="l",lty=2, main="",lwd=2)
axis(2, at=c(0,1000),line = 3.5, lwd = 2, labels = F)
axis(2, ylim=c(0,max(graph.dat$Depth, na.rm = T)),lwd=2,line=3.5)
mtext(2,text="Sonde Depth (m)",line=5.5)
#plot salinity
par(new=T)
plot(Sal~DateTime, data=graph.dat, axes=F, ylim=c(0,max(graph.dat$Sal, na.rm = T)), xlab="", ylab="", 
     type="l",lty=3, main="",lwd=2)
axis(2, at=c(0,1000), lwd=2, line = 7, labels = F)
axis(2, ylim=c(0,max(graph.dat$Sal, na.rm = T)),lwd=2,line=7)
mtext(2,text="Salinity (psu)",line=9)
#plot x axis
#make a text string with date/times
axisdate <- seq(min(graph.dat$DateTime), max(graph.dat$DateTime), by = '2 hours')
axisdate <- format(axisdate, '%m/%d %H:%M')
#draw the axis
axis.POSIXct(1,range(graph.dat$DateTime), labels=F)
axis(1, at=c(0,50000000000), labels=F)
#add in the text string below the axis
text(x=(graph.dat$DateTime), par("usr")[3] - 1, srt = 45, adj = 1,labels = axisdate, xpd = TRUE)
#plot the legend
legend(x = 'topright', inset = -0.05, legend=c("PO4","Depth","Sal"),lty=c(1,2,3), xpd = TRUE)
###-------------END GRAPH-----------------------
###--------------START GRAPH-----------------
###
###
#CHANGE ME
graph.dat <- Oct.dat
###
###
#plot PO4
plot(PO4~DateTime, data=graph.dat, axes=F, ylim=c(0, max(graph.dat$PO4, na.rm = T)), 
     xlab="", ylab="", type="l", col="black", main="")
axis(2, at=c(0,1000),lwd=2, labels = F)
axis(2, ylim=c(0,max(graph.dat$PO4, na.rm = T)), col="black",lwd=2)
mtext(2,text="PO4 conc (mg/L)",line=2)
#plot depth
par(new=T)
plot(Depth~DateTime, data=graph.dat, axes=F, ylim=c(0,max(graph.dat$Depth, na.rm = T)), xlab="", ylab="", 
     type="l",lty=2, main="",lwd=2)
axis(2, at=c(0,1000),line = 3.5, lwd = 2, labels = F)
axis(2, ylim=c(0,max(graph.dat$Depth, na.rm = T)),lwd=2,line=3.5)
mtext(2,text="Sonde Depth (m)",line=5.5)
#plot salinity
par(new=T)
plot(Sal~DateTime, data=graph.dat, axes=F, ylim=c(0,max(graph.dat$Sal, na.rm = T)), xlab="", ylab="", 
     type="l",lty=3, main="",lwd=2)
axis(2, at=c(0,1000), lwd=2, line = 7, labels = F)
axis(2, ylim=c(0,max(graph.dat$Sal, na.rm = T)),lwd=2,line=7)
mtext(2,text="Salinity (psu)",line=9)
#plot x axis
#make a text string with date/times
axisdate <- seq(min(graph.dat$DateTime), max(graph.dat$DateTime), by = '2 hours')
axisdate <- format(axisdate, '%m/%d %H:%M')
#draw the axis
axis.POSIXct(1,range(graph.dat$DateTime), labels=F)
axis(1, at=c(0,50000000000), labels=F)
#add in the text string below the axis
text(x=(graph.dat$DateTime), par("usr")[3] - 1, srt = 45, adj = 1,labels = axisdate, xpd = TRUE)
#plot the legend
legend(x = 'topright', inset = -0.05, legend=c("PO4","Depth","Sal"),lty=c(1,2,3), xpd = TRUE)
###-------------END GRAPH-----------------------
###--------------START GRAPH-----------------
###
###
#CHANGE ME
graph.dat <- Nov.dat
###
###
#plot PO4
plot(PO4~DateTime, data=graph.dat, axes=F, ylim=c(0, max(graph.dat$PO4, na.rm = T)), 
     xlab="", ylab="", type="l", col="black", main="")
axis(2, at=c(0,1000),lwd=2, labels = F)
axis(2, ylim=c(0,max(graph.dat$PO4, na.rm = T)), col="black",lwd=2)
mtext(2,text="PO4 conc (mg/L)",line=2)
#plot depth
par(new=T)
plot(Depth~DateTime, data=graph.dat, axes=F, ylim=c(0,max(graph.dat$Depth, na.rm = T)), xlab="", ylab="", 
     type="l",lty=2, main="",lwd=2)
axis(2, at=c(0,1000),line = 3.5, lwd = 2, labels = F)
axis(2, ylim=c(0,max(graph.dat$Depth, na.rm = T)),lwd=2,line=3.5)
mtext(2,text="Sonde Depth (m)",line=5.5)
#plot salinity
par(new=T)
plot(Sal~DateTime, data=graph.dat, axes=F, ylim=c(0,max(graph.dat$Sal, na.rm = T)), xlab="", ylab="", 
     type="l",lty=3, main="",lwd=2)
axis(2, at=c(0,1000), lwd=2, line = 7, labels = F)
axis(2, ylim=c(0,max(graph.dat$Sal, na.rm = T)),lwd=2,line=7)
mtext(2,text="Salinity (psu)",line=9)
#plot x axis
#make a text string with date/times
axisdate <- seq(min(graph.dat$DateTime), max(graph.dat$DateTime), by = '2 hours')
axisdate <- format(axisdate, '%m/%d %H:%M')
#draw the axis
axis.POSIXct(1,range(graph.dat$DateTime), labels=F)
axis(1, at=c(0,50000000000), labels=F)
#add in the text string below the axis
text(x=(graph.dat$DateTime), par("usr")[3] - 1, srt = 45, adj = 1,labels = axisdate, xpd = TRUE)
#plot the legend
legend(x = 'topright', inset = -0.05, legend=c("PO4","Depth","Sal"),lty=c(1,2,3), xpd = TRUE)
###-------------END GRAPH-----------------------
dev.off()

Have fun with this!

Attachments:

...lot.pdf (5KB)

Last Edit: 10 months 3 weeks ago by Todd.OBrien. Reason: Standardizing Titles

The administrator has disabled public write access.

Plot of the Month #2: Multiple y-axes 1 year 2 weeks ago #23

Marcus Beck Offline Administrator Posts: 32 Thank you received: 5 Karma: 4	Great plots Kim! I plan on creating a function to include in SWMPr that does something similar. Stay tuned!
	The administrator has disabled public write access.

Plot of the Month #2: Multiple y-axes 1 year 2 weeks ago #24

Kim_Cressman
Offline
Moderator
Posts: 19
Karma: 2

Thanks Marcus! A SWMPr function for this is going to be AWESOME!

I've been playing some with the legend position and had a breakthrough.

In my original script, I had the legend in the top right corner, with x='topright' and inset= -0.05 (which was in the original blog post I was copying from).

Here's what happened when I changed it to x='topleft':

Me:

I looked at the code some more and decided that negative inset was probably pushing it to the left, so I started playing with it. And sure enough, inset=0.0 put it right on the axis:

I ended up liking an inset of 0.01 the best for this. I also shrank the text just a bit by adding cex=0.8 into the script.

Final code:

legend(x = 'topleft', inset = 0.01, legend=c("PO4","Depth","Sal"),lty=c(1,2,3), cex=0.8, xpd = TRUE)

Final plot:

Final me:

Last Edit: 1 year 2 weeks ago by Kim_Cressman.

The administrator has disabled public write access.

Forum

The Plot of the Month

POTM: Multiple Variable Exploration

POTM-02 (June): Multiple y-axes

Time to create page: 0.151 seconds