R for Statistical Analysis

Home

Extent

Volume

Flux

R for Statistical Analysis

Data

Arctic Sea Ice
by Tracy Rogers
4-21-14
UAF Phys 212

#R for Statistical Analysis
#Download R at the R Project for Statistical Computing
#This page contains all the code (plus some) used to make the figures and calculations
#You can copy this text into R, and, with the data (in data) you can copy sections of this code
#directly into R
#I use Notepad++ to write code, and then copy that code into R.
#This allows me to compile several commands and then run it all in one go.
#This is particularly useful when running multiple trials of something, but with different data.
# the # symbol is used for intertextual notes when writing in R. This line won't do anything.

#For any command, type > help(command)
#This will bring up information on that command
dec.msie<-c(15.6,15.3,15.6,15.4)
#<- is the same as =, and is used for storing data in a variable.
#In this case, the variable name is dec.msie
#This is storing the values of March decadal sea ice extent
dec.ssie<-c(7.1,6.3,6.4,5)
dec.asiv<-c(32,30,27,24.5)
dec.ssiv<-c(16.5,14,11,5)
#These are sample datasets to practice working in R, and cover 1980, 1990, 2000, 2010 SIE and SIV

dec.year<-c(1980,1990,2000,2010)
plot(dec.year,dec.msie)
#This is a simple plot function, but we want to make this a little fancier
plot(dec.year,dec.msie,ylim=c(12,16),xlab='Year',ylab='March Decadal SIE in million sq km',ty='l')
#With this dataset, we clearly don't have significant a trend, so it would be pointless to add a regression.

plot(dec.year,dec.ssie,ylim=c(3,8),xlab='Year',ylab='Sept Decadal SIE in million sq km',ty='l')
#However, for this, we can see something of a trend.

#To look at the correlation between sea ice data, we can use:
summary(lm(dec.ssie~dec.year))
#lm is the regression itself, but summary prints out more information about it,
#including the intercept, slope, and statistical significance

#If we want to plot that out, we can use the lines command, and then $fitted.values after
#lm to plot the fitted values of that regression line
#In R, adding $ after a variable calls specific elements of that variable
#In this case, we are interested in the fitted values of the regression line

lines(dec.year,lm(dec.ssie~dec.year)$fitted.values,col='blue')
#The lines command adds a line on top of the current plot,
#while points will add points on top of your current plot.
#In this case, we are using the regression line command, lm, to put fitted values
#in the form of a line on top of our graph. There's a simpler way to do this though:
plot(dec.year,dec.ssie,ylim=c(3,8),xlab='Year',ylab='Sept Decadal SIE in million sq km',ty='l')
abline(lm(dec.ssie~dec.year),col='blue')

plot(dec.year,dec.asiv,ylim=c(20,35),xlab='Year',ylab='April Decadal SIV in 1000 cubic km',ty='l')
abline(lm(dec.asiv~dec.year),col='blue')
plot(dec.year,dec.ssiv,ylim=c(4,20),xlab='Year',ylab='Sept Decadal SIV in 1000 cubic km',ty='l')
abline(lm(dec.ssiv~dec.year),col='blue')
#These were to show the basic decadal information without cluttering the graph too much.
#I didn't use these graphs, but this shows most of the commands I'll use for the graphs lower down

#Now we'll look at the full SIV and SIE datasets
asiv<-scan()
#scan() is a useful command that lets you enter a dataset
#There are several other ways to do this, such as read.csv(), but scan works well for just 4 datasets
#With scan(), we can take a column of excel data, and paste it directly into R
ssiv<-scan()
msie<-scan()
ssie<-scan()
#I entered the datasets in x.x 10⁶ km² for SIE, and
#xx.x 10³ km³ for SIV, which makes for better looking tables and figures

#We can also look at the mean, median, and standard deviation of our data:
mean(asiv)
median(asiv)
sd(asiv)

mean(ssiv)
median(ssiv)
sd(ssiv)

mean(ssie)
median(ssie)
sd(ssie)

mean(msie)
median(msie)
sd(msie)

year<-1979:2013
#This is binding the years 1979 to 2013 to the variable year.
#The : tells R to include those years and all years in between as a string of numbers.

plot(year,msie,ylim=c(12,18),xlab='Year',ylab='March SIE in million sq km',ty='l')
plot(year,ssie,ylim=c(3,8),xlab='Year',ylab='Sept SIE in million sq km',ty='l')
abline(lm(ssie~year),col='blue')

#Another way to analyze this data is to look at it as a percent of the 1979 data,
#which allows comparing sea ice between different months
plot(year,ssie/ssie[1],lwd=2,ylim=c(.4,1.2),xlab='Year',ylab='Sept and March SIE as percent',ty='l')
#lwd increases the width of the line
lines(year,msie/msie[1],lwd=2,col='blue')
abline(lm(msie/msie[1]~year),col='blue')
abline(lm(ssie/ssie[1]~year))
#This block of plot/lines gives SSIE and MSIE as thick black and blue lines,
#while the regression lines are thin black and blue lines

plot(year,asiv,ylim=c(20,35),xlab='Year',ylab='April SIV in 1000 cubic km',ty='l')
abline(lm(asiv~year),col='blue')
plot(year,ssiv,ylim=c(4,20),xlab='Year',ylab='Sept SIV in 1000 cubic km',ty='l')
abline(lm(ssiv~year),col='blue')

plot(year,ssiv/ssiv[1],lwd=2,ylim=c(0,1.3),xlab='Year',ylab='Sept and April SIV as percent',ty='l')
lines(year,asiv/asiv[1],lwd=2,col='blue')
abline(lm(asiv/asiv[1]~year),col='blue')
abline(lm(ssiv/ssiv[1]~year))

#Next, we are interested in the heat flux of sea ice, which is best represented
#by the difference in volume between the maximum (April) and minimum (September) months.
siv.diff<-(asiv-ssiv)*1000

#First we need to determine the mass of sea ice
#Ice density varies based on age (older ice is more dense), but we can approximate it with
#.91 Mg (megagrams)/m^3 (Timco and Frederking, 1996), which means we need to convert to m^3
siv.meter<-siv.diff*1000*1000*1000
siv.mass<-siv.meter*.91
siv.mass<-siv.mass*1000000 #megagrams to grams
siv.flux<-siv.mass*334 #mass * 334 J/g

#There are a few ways to explore the data. First, to check our work, PIOMAS lists that to melt
#16,400 km^3 of ice, it requires approximately 5*10^21 Joules. To check our values:
range(siv.flux)
# 3.95e+21 6.23e+21
#This means we are within an acceptable range.

mean(siv.flux)
median(siv.flux)
sd(siv.flux)
# 5.11e+21, 5.02e+21, 5.22e+20

plot(year,siv.flux,xlab='Year',ylab='Sea Ice Heat Flux', ty='l',lwd=2)
abline(lm(siv.flux~year))

Home