R scripts: How did I make this plot?

Figure 1: The "master plot" for the R10mm index, which tells us how many days in year year had at least 10mm of precipitation.

Figure 1: The “master plot” for the R10mm index, which tells us how many days in each year had at least 10mm of precipitation fall.

Organization and data display is proving to be integral in this climate analysis “chapter.” We have 26 climate indices, 11 of which are monthly (meaning there are up to 12 plots per index), for the 20 data products selected I this analysis!

Luckily, R is being a huge helper in data organization and visual display. For each annual-based index, I have set up a simple matrix to group all the information-packed plots per index into one window.

Below I thought I’d make a post on how I create the “Master Plot” shown in Figure 1.

Here is the code which will re-recreate the plot above, with a brief description of what it means. Note: I am new to R, so for experts, this may not be the most efficient script, but it will work!

pdf(“C:/where.you.want.to.save.this.image/name.of.image.pdf”, width=5, height=4.5)

#This line will automatically create and save a PDF of this image to your computer.

#For beginners, if you have a file on your computer you’d love these to be saved to, but don’t know the “code” to get there, try this: Right click on that file, go to Properties, and copy the Location. The last line in the code below (name.of.image.pdf) is the name you want this file to be, so just add that!

layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow = TRUE), widths=c(1,1,1,1), heights=c(1,0.5,1,1))

#This script allows you to plot multiple images in the same window. I have adjusted the ratios so all 4 images we will plot will fit nicely. The matrix, widths, and heights entered here will create the organization of images you see in Figure 1.

par(mar=c(1, 4, 2, 1))

#The par command changes the plot margins to make more room!

plot(x, y, ylim=c(15,60), col=’black’, xlim=range(longest x in your timeseries), type=’l’, main=’R10mm’,lwd=3, xlab=’Year’, ylab=’R10mm (days)’)

#Now to the plot! This is the simple time series plot. Repeat to add any more “lines” to the first panel plot using the par(new=T) line. The xlim=range(x) allows all the times series to be plotted on a common time frame.


R10mm.mean<-rollapply(R10mm.HAD, 21, mean, by=1)

 #create the 21-year rolling mean using zoo.

plot(R10mm.mean, ylim=c(30,45), col=’black’, xlim=range(longest x), type=’l’, main=NULL,lwd=3,xaxt=’n’, xlab=’Year’, ylab=’21-year mean’)

#now plot the “smoothed” mean, again on a common time range.

axis(side = 3, at = c(1900,1920,1940,1960,1980,2000))

#This axis command creates that neat “shared” time axis for space saving!

par(mar=c(4, 3, .5, 1))

#Now onto the density plots! First, we need to again change the margins to get everything to fit!

R10mm.1 <- ((coredata(R10mm.HAD))[(((index(R10mm.HAD))> “1980”)) & (((index(R10mm.HAD)) < “2011”))])

#Now, to create the probability density functions, we must “break” our time series into our specific climate normal from 1951-1980 and 1981-2010. You can change those dates to whichever time ranges you want to compare!

d.R10mm.1 <- density(R10mm.1,na.rm=T) # returns the density data

plot(d.R10mm.1,col=”black”, ylim=c(0,0.09), lwd=3,xlim=c(15,60), main=NA, xlab=”R10mm (days)”) # plots the results

 #now use the density function to create those nice probability density functions that allow us to compare changes in distribution (mean, tails, highest density value of data, etc).

NORTH<-cbind((R10mm.MD.12),(R10mm.MD.3), (R10mm.VA.5),(R10mm.MD.9), (R10mm.VA.8),(R10mm.MD.1), (R10mm.MD.11),(R10mm.MD.5), (R10mm.MD.8))

#First, I grouped the stations into North and South by combining the individual stations into two sets of 9 stations (ordered by latitude).

#Now, we repeat the density code, but for the North and South grouping GHCN-Daily stations.


#Now end you image; this saves it to your file!

Now, depending on how you label your own files and commands, and which time series you choose to plot, you have the script to create Figure 1.

Was this helpful to you? Would you like to see more R programing posts in the future? What does Figure 1 say to you? Voice your opinion in the comments!

Kari Pohl

About Kari Pohl

I am a post-doctoral researcher at NOAA and the University of Maryland (Center for Environmental Science at Horn Point Laboratory). My work investigates how climate variability and extremes affect the diverse ecosystems in Chesapeake Bay. I received a Ph.D. in oceanography from the University of Rhode Island (2014) and received a B.S. in Environmental Science and a B.A. in Chemistry from Roger Williams University (2009). When I am not busy being a scientist, my hobbies include running, watching (and often yelling at) the Boston Bruins, and taking photos of my cat.
This entry was posted in Data, R script and tagged , , , . Bookmark the permalink.

Leave a Reply

Your email address will not be published. Required fields are marked *