Code for global sea ice calcs

The code used for generating the global ice trend plots is below. Instructions are in the code. This should be easily reproducible in excel with some time spent breaking the data into annual basis for conversion into anomaly.

WordPress italicized my quotes again.  If you use R for replication you need to fix all the quotes — sorry.

I also apologize for the heavy C accent in my code.  I gave up trying to work out the R syntax after several hours.  Hints for improvement would be appreciated.

library(gplots)
source(“http://www.climateaudit.org/scripts/utilities.txt”)
nb=read.csv(“c:/agw/sea ice/global sea ice compiled.csv”, header=FALSE)

#I used microsoft excel to represent years as decimals of days. I.e. 1978.23 This was done in the CSV file before running.
#the data is in 3 columns, Year, global nasateam area summation, global bootstrap summation
#missing data may require cells to be “cleared” not just whitespaced
#data can be downloaded to excel from the follwing websites
#http://nsidc.org/data/smmr_ssmi_ancillary/area_extent.html — main page with data description
#ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/trends-climatologies/ice-extent/nasateam/ — nasa team data add .n and .s together for global numbers
#ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/trends-climatologies/ice-extent/bootstrap/ — bootstrap data add .n and .s together for global numbers

nb[,2]=nb[,2]/1000000
nb[,3]=nb[,3]/1000000

ll=length(nb[,1])

annual=array(0,dim=365)
anncount=array(0,dim=365)
anomaly=array(0,dim=ll)
anomalyb=array(0,dim=ll)
nbts=ts(nb[,2],start=nb[1,1])

##########nasateam plots
#break data into annual basis
for(j in 1:ll)
{
ind=as.integer((nb[j,1]-as.integer(nb[j,1]))*365+.5)+1
if(!is.na(nb[j,2]))
{
annual[ind]=annual[ind]+nb[j,2]
anncount[ind]=anncount[ind]+1
}
}
#calc annualied average
for(j in 1:365)
{
if(anncount[j]!=0)
{
annual[j]=(annual[j]/anncount[j])
}
}

#use average to calc anomaly
for(j in 1:ll)
{
ind=as.integer((nb[j,1]-as.integer(nb[j,1]))*365+.5)+1

anomaly[j]=nb[j,2]-annual[ind]
}

#plot results
plot(annual,type=”l”,ylab=”Ice Area 1e6 Km^2″, xlab=”Day of Year”,main=”Global 30 yr Average Sea Ice Area\nNasaTeam Algorithm”)
grid()
savePlot(paste(“c:/agw/sea ice/Global 30 year average NasaTeam Algorithm AREA.jpg”),type=”jpg”)

meanwind=nb[,1]<2000
anomalyb=anomaly-mean(anomaly[meanwind],na.rm=TRUE)
plot(nb[,1],anomaly,type=”l”,ylab=”Ice Area Anomaly 1e6 Km^2″, xlab=”Year”,main=”Global Sea Ice Area Anomaly\nNasaTeam Algorithm”)

grid()
abline(lsfit(nb[,1], anomaly,NULL), col=”red”)
coef(lsfit(nb[,1], anomaly,NULL))
savePlot(paste(“c:/agw/sea ice/Global Sea Ice Area Anomaly NasaTeam Algorithm.jpg”),type=”jpg”)

plot(nb[,1], nb[,2], type = “l”, ylab=”Ice Area 1e6 Km^2″,xlab = “Year”,main=”Global Sea Ice Area\nNasaTeam Algorithm”,ylim=c(0,24))
grid()
abline(lsfit(nb[,1], nb[,2],NULL), col=”red”)
coef(lsfit(nb[,1], nb[,2],NULL))
mean(nb[,2],na.rm=TRUE)
savePlot(paste(“c:/agw/sea ice/Global Sea Ice NasaTeam Algorithm AREA.jpg”),type=”jpg”)
###########bootstrap plots
#break data into annual basis

annual[1:365]=0
anncount[1:365]=0

for(j in 1:ll)
{
ind=as.integer((nb[j,1]-as.integer(nb[j,1]))*365+.5)+1
if(!is.na(nb[j,3]))
{
annual[ind]=annual[ind]+nb[j,3]
anncount[ind]=anncount[ind]+1
}
}
#calc annual average
for(j in 1:365)
{
if(anncount[j]!=0)
{
annual[j]=annual[j]/anncount[j]
}
}
#calc anomaly
for(j in 1:ll)
{
ind=as.integer((nb[j,1]-as.integer(nb[j,1]))*365+.5)+1

anomalyb[j]=(nb[j,3]-annual[ind])
}

plot(annual,type=”l”,ylab=”Ice Area 1e6 Km^2″, xlab=”Day of Year”,main=”Global 30 yr Average Sea Ice Area\nBootstrap Algorithm”)
grid()
savePlot(paste(“c:/agw/sea ice/Global 30 year average bootstrap Algorithm AREA.jpg”),type=”jpg”)

meanwind=nb[,1]<2000
anomalyb=anomalyb-mean(anomalyb[meanwind],na.rm=TRUE)
plot(nb[,1],anomalyb,type=”l”,ylab=”Ice Area Anomaly 1e6 Km^2″, xlab=”Year”,main=”Global Sea Ice Area Anomaly\nBootstrap Algorithm”,ylim=c(-3,3))
grid()
abline(lsfit(nb[,1], anomalyb,NULL), col=”red”)
coef(lsfit(nb[,1], anomalyb,NULL))

savePlot(paste(“c:/agw/sea ice/Global Sea Ice Area Anomaly bootstrap Algorithm.jpg”),type=”jpg”)

plot(nb[,1], nb[,3], type = “l”, ylab=”Ice Area 1e6 Km^2″,xlab = “Year”,main=”Global Sea Ice Area\nBootstrap Algorithm”,ylim=c(0,24))
grid()
abline(lsfit(nb[,1], nb[,3],NULL), col=”red”)
b=nb[,1]<2008
coef(lsfit(nb[,1], nb[,3],NULL))
coef(lsfit(nb[,1][b], nb[,3][b],NULL))
abline(lsfit(nb[,1][b], nb[,3][b],NULL), col=”blue”)
mean(nb[,3],na.rm=TRUE)
savePlot(paste(“c:/agw/sea ice/Global Sea Ice Area bootstrap Algorithm.jpg”),type=”jpg”)