jstults/ice-code.R

## ice-code.R
ssq = function (x) {sum(x ^ 2)}

plotTrend = function (x, st = 1978, en = c(2012, 12), y.pos = NA, x.pos = NA, main.t = "Untitled") {

	### Get trend
	trend = lm(window(x, st, en) ~ I(time(window(x, st, en))))

	### Initialize variables
	N = length(window(x, st, en))
	I = seq(1:N) / frequency(x)
	rmI = ssq(I - mean(I))

	### Calculate sum of squared errors
	SSE = ssq(trend[[2]]) / (N - 2)

	### Calculate OLS standard errors
	seOLS = sqrt(SSE / rmI) * 10

	### Calculate two-tailed 95% CI
	ci95 = seOLS * 1.964

	### Get lag-1 ACF
	resids = residuals(trend)
	r1 = acf(resids, lag.max=1, plot=FALSE)$acf[[2]]

	### Calculate CIs with Quenouille (Santer) adjustment for autocorrelation
	Q = (1 - r1) / (1 + r1)
	ciQ = ci95 / sqrt(Q)

	### Plot data
	plot(window(x, st, en), main=main.t, ylab="Anomaly Km^2")
	grid()

	### Add trendline and x-axis
	abline(h = 0, col = 2)
	abline(trend, col = 4, lwd = 2)

	### Annotate text
	text(x = ifelse(is.na(x.pos), st, x.pos), y = ifelse(is.na(y.pos), max(x, na.rm = TRUE) - 0.2 * max(x, na.rm = TRUE), y.pos), paste("Slope: ", signif(trend$coef[2] * 10, 3), "+/-", signif(ciQ, 3), "Km^2/Decade\n(corrected for AR(1) serial correlation)\nLag-1 value:", signif(r1, 3)), cex = 0.8, col = 4, pos = 4)

	###  Return
	r = list(trend$coef[2] * 10, ciQ)
}


cal.Isig=function(dat=dat, st=0,en=2020)
{
	###  Get trend
	fm=lm(window(dat, st, en)~I(time(window(dat, st, en))))

	###  Initialize variables
	N=length(window(dat, st, en))
	I=seq(1:N)/frequency(dat)
	rmI=ssq(I-mean(I))

	###  Calculate sum of squared errors
	SSE=ssq(fm[[2]])/(N-2)

	###  Calculate OLS standard errors
	seOLS=sqrt(SSE/rmI)

	###  Calculate two-tailed 95% CI
	ci95=seOLS*1.964

	###  Get lag-1 ACF
	resids=residuals(lm(window(dat, st, en)~seq(1:N)))
	r1=acf(resids, lag.max=1, plot=FALSE)$acf[[2]]

	###  Calculate CIs with Quenouille (Santer) adjustment for autocorrelation
	Q=(1-r1)/(1+r1)
	ciQ=ci95/sqrt(Q)

	return(c(fm$coef[2],ciQ))
}

##################################################################################
##################################################################################
##################################################################################
##################################################################################
##################################################################################
##################################################################################
#calculate northern sea ice area from gridded satelite data

Nfilenames=list.files(path="C:/agw/sea ice data/north/", pattern = NULL, all.files = TRUE, full.names = FALSE, recursive = TRUE)
startday=as.integer(unclass(as.POSIXct(strptime("1977-12-31", '%Y-%m-%d' )))/86400)
year= substr(Nfilenames,9,12)
month=substr(Nfilenames,13,14)
day=substr(Nfilenames,15,16)
dat=paste(year,month,day,sep="-")
oday=as.integer(unclass(as.POSIXct(strptime(dat, '%Y-%m-%d' )))/86400)-startday

area=array(NA,dim=c(366,35))
offsetday=array(NA,dim=c(366,35))

#set up hole mask here
getnorthicedata=function(file=file)
{
	a=file(file,"rb")
	header= readBin(a,n=300,what=integer(),size=1,endian="little",signed=FALSE)
	data=readBin(a,n=304*448,what=integer(),size=1,endian="little",signed=FALSE)
	close(a)
	data
}

starthmask=!(getnorthicedata(paste("C:/agw/sea ice data/north/",Nfilenames[10],sep=""))==251)
endhmask=!(getnorthicedata(paste("C:/agw/sea ice data/north/",Nfilenames[9000],sep=""))==251)
diffmask=!(starthmask==endhmask)
holemask=starthmask #insures no step

lat=readBin("C:/agw/sea ice data/sea ice coordinates/psn25lats_v3.dat", integer(), endian = "little",n=304*448)/100000
lon=readBin("C:/agw/sea ice data/sea ice coordinates/psn25lons_v3.dat", integer(), endian = "little",n=304*448)/100000
latmask=lat>66.5  #arctic circle


#################################
#load area values into 2d matrix 366days each year x 35 days
for(i in 1:(length(Nfilenames)-1))
{
	#open filename
	fn=paste("C:/agw/sea ice data/north/",Nfilenames[i],sep="")
	data=getnorthicedata(fn)
	print(dat[i])
	#Limit data/land mask to all info <250 - 100% and >25 - 10%
	datamask=data<251# & data>25

	#store sum in trend value
	#area[ydayindex[i],yearindex[i]]=sum(data[(datamask&holemask&latmask)])/250*625
	area[ydayindex[i],yearindex[i]]=sum(data[(datamask&holemask)])/250*625
	offsetday[ydayindex[i],yearindex[i]]=oday[i]
}

areatimeseries=area
dim(areatimeseries)=c(366*35,1)
plot(areatimeseries)

fareatimeseries=rep(NA,length(areatimeseries))
######filter timeseries for anomaly calculation only
######Use more recent years to avoid missing values from early record
fwid=10
for (i in 4000:length(areatimeseries))
{
	if(!is.na(areatimeseries[i]))
	{
		miave=i-10
		maave=i+10

		fareatimeseries[i]=mean(areatimeseries[miave:maave],na.rm=TRUE)
	}
}

#resequence all matrix data into time series to complete correction for leap years
dim(fareatimeseries)=c(366,35)
anomval=rowMeans(fareatimeseries,na.rm=TRUE)
NIAnom=area-anomval
NIArea=area
dim(NIAnom)=c(366*35,1)
dim(NIArea)=c(366*35,1)
dim(offsetday)=c(366*35,1)
NorthIceArea=rep(NA,366*35)
NorthIceAreaAnomaly=rep(NA,366*35)

for (i in 1:length(NorthIceArea))
{
	if(!is.na(offsetday[i]))
	{
		NorthIceArea[offsetday[i]]=NIArea[i]
		NorthIceAreaAnomaly[offsetday[i]]=NIAnom[i]

	}
}

NorthIceArea=ts(NorthIceArea,start=1978,deltat=1/365.25)  #set time series to get time values only
decdate=time(NorthIceArea)						#copy time values
NorthIceArea=approx(decdate,NorthIceArea,xout=decdate)	#interpolate missing values
NorthIceArea=ts(NorthIceArea$y,start=1978,deltat=1/365.25)	#reassemble series
NorthIceArea=NorthIceArea + sum(!holemask)*625
NorthIceAreaAnomaly=approx(decdate,NorthIceAreaAnomaly,xout=decdate)		#interpolate missing values
NorthIceAreaAnomaly=ts(NorthIceAreaAnomaly$y,start=1978,deltat=1/365.25)	#reassemble series
NorthIceAreaAnomaly=NorthIceAreaAnomaly-mean(NorthIceAreaAnomaly,na.rm=TRUE)	#center series


########## plotting calls may be a little messy due to continuous changes
plotTrend(NorthIceAreaAnomaly,y.pos=-1.4e6,main.t = "Arctic Ice Area Anomaly")
#arcticcircleanomaly=NorthIceAreaAnomaly
#arcticcirclearea=NorthIceArea
#totalanomaly=NorthIceAreaAnomaly
#totalarea=NorthIceArea

plotTrend(NorthIceAreaAnomaly,y.pos=-1.3e6,main.t = "Sea Ice Area North")

plot(arcticcirclearea,y.pos=-.7e6,main = "Arctic Circle Sea Ice Area",ylim=c(0,1.3e7),ylab="Ice Area km^2")
grid()
plotTrend(arcticcircleanomaly,y.pos=-.9e6,main.t = "Arctic Circle Sea Ice Area Anomaly")
plotTrend(NorthIceArea-above60area,y.pos=-.7e6,main.t = "Sea Ice Area Below Arctic Circle")
plotTrend(NorthIceAreaAnomaly-above60anomaly,y.pos=-.7e6,main.t = "Arctic Sea Ice Area Anomaly Below Arctic Circle")
lines(totalanomaly,col="red")

plotTrend(NorthIceAreaAnomaly,y.pos=-1.2e6,main.t = "Northern Hemisphere Sea Ice Area Anomaly Total vs Non-Arctic")
lines(NorthIceAreaAnomaly-above60anomaly,col="red")

plotTrend(NorthIceAreaAnomaly,y.pos=-1.2e6,main.t = "Northern Hemisphere Sea Ice Area Anomaly")

nonarctic=NorthIceArea-arcticcirclearea

seasonalratio = nonarctic/arcticcirclearea*100
plot(seasonalratio, main="NH Ratio of Non-Arctic Ice Area to Arctic Ice Area",ylab="Percent")
grid()

##########################################
##########################################
#southern hemisphere sea ice
Sfilenames=list.files(path="C:/agw/sea ice data/south/", pattern = NULL, all.files = TRUE, full.names = FALSE, recursive = TRUE)
startday=as.integer(unclass(as.POSIXct(strptime("1977-12-31", '%Y-%m-%d' )))/86400)
year= substr(Sfilenames,9,12)
month=substr(Sfilenames,13,14)
day=substr(Sfilenames,15,16)
dat=paste(year,month,day,sep="-")
oday=as.integer(unclass(as.POSIXct(strptime(dat, '%Y-%m-%d' )))/86400)-startday

area=array(NA,dim=c(366,35))
offsetday=array(NA,dim=c(366,35))

getsouthicedata=function(file=file)
{
	a=file(file,"rb")
	header= readBin(a,n=300,what=integer(),size=1,endian="little",signed=FALSE)
	data=readBin(a,n=332*316,what=integer(),size=1,endian="little",signed=FALSE)
	close(a)
	data
}

#################################
#load area values into 2d matrix 366days each year x 35 days
for(i in 1:(length(Sfilenames)-1))
{
	#open filename
	fn=paste("C:/agw/sea ice data/south/",Sfilenames[i],sep="")
	data=getsouthicedata(fn)
	print(dat[i])
	#Limit data/land mask to all info <250 - 100% and >25 - 10%
	datamask=data<251# & data>25

	#store sum in trend value
	area[ydayindex[i],yearindex[i]]=sum(data[datamask])/250*625
	offsetday[ydayindex[i],yearindex[i]]=oday[i]
}

areatimeseries=area
dim(areatimeseries)=c(366*35,1)
fareatimeseries=rep(NA,length(areatimeseries))

######filter timeseries for anomaly calculation only
######Use more recent years to avoid missing values from early record
fwid=10
for (i in 4000:length(areatimeseries))
{
	if(!is.na(areatimeseries[i]))
	{
		miave=i-10
		maave=i+10

		fareatimeseries[i]=mean(areatimeseries[miave:maave],na.rm=TRUE)
	}
}

#resequence all matrix data into time series to complete correction for leap years
dim(fareatimeseries)=c(366,35)
anomval=rowMeans(fareatimeseries,na.rm=TRUE)
SIAnom=area-anomval
SIArea=area
dim(SIAnom)=c(366*35,1)
dim(SIArea)=c(366*35,1)
dim(offsetday)=c(366*35,1)
SouthIceArea=rep(NA,366*35)
SouthIceAreaAnomaly=rep(NA,366*35)

#remove two extreme values which are obviously errors
mask= SIAnom< -5000000  | SIAnom>3000000
SIAnom[mask]=NA
SIArea[mask]=NA

for (i in 1:length(SouthIceArea))
{
	if(!is.na(offsetday[i]))
	{
		SouthIceArea[offsetday[i]]=SIArea[i]
		SouthIceAreaAnomaly[offsetday[i]]=SIAnom[i]

	}
}


SouthIceArea=ts(SouthIceArea,start=1978,deltat=1/365.25)  #set time series to get time values only
decdate=time(SouthIceArea)						#copy time values
SouthIceArea=approx(decdate,SouthIceArea,xout=decdate)	#interpolate missing values
SouthIceArea=ts(SouthIceArea$y,start=1978,deltat=1/365.25)	#reassemble series
SouthIceAreaAnomaly=approx(decdate,SouthIceAreaAnomaly,xout=decdate)		#interpolate missing values
SouthIceAreaAnomaly=ts(SouthIceAreaAnomaly$y,start=1978,deltat=1/365.25)	#reassemble series
SouthIceAreaAnomaly=SouthIceAreaAnomaly-mean(SouthIceAreaAnomaly,na.rm=TRUE)	#center series

plotTrend(SouthIceAreaAnomaly,y.pos=-2e6,main.t = "Antarctic Ice Area Anomaly")


GlobalIceArea=NorthIceArea+SouthIceArea
mean(GlobalIceArea,na.rm=TRUE)#19020751
GlobalIceAreaAnomaly= NorthIceAreaAnomaly+SouthIceAreaAnomaly
plotTrend(GlobalIceAreaAnomaly,y.pos=-2e6,main.t = "Global Ice Area Anomaly")

plot(GlobalIceAreaAnomaly+19000000,ylim=c(0,22000000),main="Global Sea Ice Area Anomaly\nOffset by Average Global Ice 19M Km^2", ylab="Km^2 Sea Ice")
grid()
	ssq = function (x) {sum(x ^ 2)}

	plotTrend = function (x, st = 1978, en = c(2012, 12), y.pos = NA, x.pos = NA, main.t = "Untitled") {

	### Get trend
	trend = lm(window(x, st, en) ~ I(time(window(x, st, en))))

	### Initialize variables
	N = length(window(x, st, en))
	I = seq(1:N) / frequency(x)
	rmI = ssq(I - mean(I))

	### Calculate sum of squared errors
	SSE = ssq(trend[[2]]) / (N - 2)

	### Calculate OLS standard errors
	seOLS = sqrt(SSE / rmI) * 10

	### Calculate two-tailed 95% CI
	ci95 = seOLS * 1.964

	### Get lag-1 ACF
	resids = residuals(trend)
	r1 = acf(resids, lag.max=1, plot=FALSE)$acf[[2]]

	### Calculate CIs with Quenouille (Santer) adjustment for autocorrelation
	Q = (1 - r1) / (1 + r1)
	ciQ = ci95 / sqrt(Q)

	### Plot data
	plot(window(x, st, en), main=main.t, ylab="Anomaly Km^2")
	grid()

	### Add trendline and x-axis
	abline(h = 0, col = 2)
	abline(trend, col = 4, lwd = 2)

	### Annotate text
	text(x = ifelse(is.na(x.pos), st, x.pos), y = ifelse(is.na(y.pos), max(x, na.rm = TRUE) - 0.2 * max(x, na.rm = TRUE), y.pos), paste("Slope: ", signif(trend$coef[2] * 10, 3), "+/-", signif(ciQ, 3), "Km^2/Decade\n(corrected for AR(1) serial correlation)\nLag-1 value:", signif(r1, 3)), cex = 0.8, col = 4, pos = 4)

	### Return
	r = list(trend$coef[2] * 10, ciQ)
	}


	cal.Isig=function(dat=dat, st=0,en=2020)
	{
	### Get trend
	fm=lm(window(dat, st, en)~I(time(window(dat, st, en))))

	### Initialize variables
	N=length(window(dat, st, en))
	I=seq(1:N)/frequency(dat)
	rmI=ssq(I-mean(I))

	### Calculate sum of squared errors
	SSE=ssq(fm[[2]])/(N-2)

	### Calculate OLS standard errors
	seOLS=sqrt(SSE/rmI)

	### Calculate two-tailed 95% CI
	ci95=seOLS*1.964

	### Get lag-1 ACF
	resids=residuals(lm(window(dat, st, en)~seq(1:N)))
	r1=acf(resids, lag.max=1, plot=FALSE)$acf[[2]]

	### Calculate CIs with Quenouille (Santer) adjustment for autocorrelation
	Q=(1-r1)/(1+r1)
	ciQ=ci95/sqrt(Q)

	return(c(fm$coef[2],ciQ))
	}

	##################################################################################
	##################################################################################
	##################################################################################
	##################################################################################
	##################################################################################
	##################################################################################
	#calculate northern sea ice area from gridded satelite data

	Nfilenames=list.files(path="C:/agw/sea ice data/north/", pattern = NULL, all.files = TRUE, full.names = FALSE, recursive = TRUE)
	startday=as.integer(unclass(as.POSIXct(strptime("1977-12-31", '%Y-%m-%d' )))/86400)
	year= substr(Nfilenames,9,12)
	month=substr(Nfilenames,13,14)
	day=substr(Nfilenames,15,16)
	dat=paste(year,month,day,sep="-")
	oday=as.integer(unclass(as.POSIXct(strptime(dat, '%Y-%m-%d' )))/86400)-startday

	area=array(NA,dim=c(366,35))
	offsetday=array(NA,dim=c(366,35))

	#set up hole mask here
	getnorthicedata=function(file=file)
	{
	a=file(file,"rb")
	header= readBin(a,n=300,what=integer(),size=1,endian="little",signed=FALSE)
	data=readBin(a,n=304*448,what=integer(),size=1,endian="little",signed=FALSE)
	close(a)
	data
	}

	starthmask=!(getnorthicedata(paste("C:/agw/sea ice data/north/",Nfilenames[10],sep=""))==251)
	endhmask=!(getnorthicedata(paste("C:/agw/sea ice data/north/",Nfilenames[9000],sep=""))==251)
	diffmask=!(starthmask==endhmask)
	holemask=starthmask #insures no step

	lat=readBin("C:/agw/sea ice data/sea ice coordinates/psn25lats_v3.dat", integer(), endian = "little",n=304*448)/100000
	lon=readBin("C:/agw/sea ice data/sea ice coordinates/psn25lons_v3.dat", integer(), endian = "little",n=304*448)/100000
	latmask=lat>66.5 #arctic circle


	#################################
	#load area values into 2d matrix 366days each year x 35 days
	for(i in 1:(length(Nfilenames)-1))
	{
	#open filename
	fn=paste("C:/agw/sea ice data/north/",Nfilenames[i],sep="")
	data=getnorthicedata(fn)
	print(dat[i])
	#Limit data/land mask to all info <250 - 100% and >25 - 10%
	datamask=data<251# & data>25

	#store sum in trend value
	#area[ydayindex[i],yearindex[i]]=sum(data[(datamask&holemask&latmask)])/250*625
	area[ydayindex[i],yearindex[i]]=sum(data[(datamask&holemask)])/250*625
	offsetday[ydayindex[i],yearindex[i]]=oday[i]
	}

	areatimeseries=area
	dim(areatimeseries)=c(366*35,1)
	plot(areatimeseries)

	fareatimeseries=rep(NA,length(areatimeseries))
	######filter timeseries for anomaly calculation only
	######Use more recent years to avoid missing values from early record
	fwid=10
	for (i in 4000:length(areatimeseries))
	{
	if(!is.na(areatimeseries[i]))
	{
	miave=i-10
	maave=i+10

	fareatimeseries[i]=mean(areatimeseries[miave:maave],na.rm=TRUE)
	}
	}

	#resequence all matrix data into time series to complete correction for leap years
	dim(fareatimeseries)=c(366,35)
	anomval=rowMeans(fareatimeseries,na.rm=TRUE)
	NIAnom=area-anomval
	NIArea=area
	dim(NIAnom)=c(366*35,1)
	dim(NIArea)=c(366*35,1)
	dim(offsetday)=c(366*35,1)
	NorthIceArea=rep(NA,366*35)
	NorthIceAreaAnomaly=rep(NA,366*35)

	for (i in 1:length(NorthIceArea))
	{
	if(!is.na(offsetday[i]))
	{
	NorthIceArea[offsetday[i]]=NIArea[i]
	NorthIceAreaAnomaly[offsetday[i]]=NIAnom[i]

	}
	}

	NorthIceArea=ts(NorthIceArea,start=1978,deltat=1/365.25) #set time series to get time values only
	decdate=time(NorthIceArea) #copy time values
	NorthIceArea=approx(decdate,NorthIceArea,xout=decdate) #interpolate missing values
	NorthIceArea=ts(NorthIceArea$y,start=1978,deltat=1/365.25) #reassemble series
	NorthIceArea=NorthIceArea + sum(!holemask)*625
	NorthIceAreaAnomaly=approx(decdate,NorthIceAreaAnomaly,xout=decdate) #interpolate missing values
	NorthIceAreaAnomaly=ts(NorthIceAreaAnomaly$y,start=1978,deltat=1/365.25) #reassemble series
	NorthIceAreaAnomaly=NorthIceAreaAnomaly-mean(NorthIceAreaAnomaly,na.rm=TRUE) #center series



	########## plotting calls may be a little messy due to continuous changes
	plotTrend(NorthIceAreaAnomaly,y.pos=-1.4e6,main.t = "Arctic Ice Area Anomaly")
	#arcticcircleanomaly=NorthIceAreaAnomaly
	#arcticcirclearea=NorthIceArea
	#totalanomaly=NorthIceAreaAnomaly
	#totalarea=NorthIceArea

	plotTrend(NorthIceAreaAnomaly,y.pos=-1.3e6,main.t = "Sea Ice Area North")

	plot(arcticcirclearea,y.pos=-.7e6,main = "Arctic Circle Sea Ice Area",ylim=c(0,1.3e7),ylab="Ice Area km^2")
	grid()
	plotTrend(arcticcircleanomaly,y.pos=-.9e6,main.t = "Arctic Circle Sea Ice Area Anomaly")
	plotTrend(NorthIceArea-above60area,y.pos=-.7e6,main.t = "Sea Ice Area Below Arctic Circle")
	plotTrend(NorthIceAreaAnomaly-above60anomaly,y.pos=-.7e6,main.t = "Arctic Sea Ice Area Anomaly Below Arctic Circle")
	lines(totalanomaly,col="red")

	plotTrend(NorthIceAreaAnomaly,y.pos=-1.2e6,main.t = "Northern Hemisphere Sea Ice Area Anomaly Total vs Non-Arctic")
	lines(NorthIceAreaAnomaly-above60anomaly,col="red")

	plotTrend(NorthIceAreaAnomaly,y.pos=-1.2e6,main.t = "Northern Hemisphere Sea Ice Area Anomaly")

	nonarctic=NorthIceArea-arcticcirclearea

	seasonalratio = nonarctic/arcticcirclearea*100
	plot(seasonalratio, main="NH Ratio of Non-Arctic Ice Area to Arctic Ice Area",ylab="Percent")
	grid()

	##########################################
	##########################################
	#southern hemisphere sea ice
	Sfilenames=list.files(path="C:/agw/sea ice data/south/", pattern = NULL, all.files = TRUE, full.names = FALSE, recursive = TRUE)
	startday=as.integer(unclass(as.POSIXct(strptime("1977-12-31", '%Y-%m-%d' )))/86400)
	year= substr(Sfilenames,9,12)
	month=substr(Sfilenames,13,14)
	day=substr(Sfilenames,15,16)
	dat=paste(year,month,day,sep="-")
	oday=as.integer(unclass(as.POSIXct(strptime(dat, '%Y-%m-%d' )))/86400)-startday

	area=array(NA,dim=c(366,35))
	offsetday=array(NA,dim=c(366,35))

	getsouthicedata=function(file=file)
	{
	a=file(file,"rb")
	header= readBin(a,n=300,what=integer(),size=1,endian="little",signed=FALSE)
	data=readBin(a,n=332*316,what=integer(),size=1,endian="little",signed=FALSE)
	close(a)
	data
	}

	#################################
	#load area values into 2d matrix 366days each year x 35 days
	for(i in 1:(length(Sfilenames)-1))
	{
	#open filename
	fn=paste("C:/agw/sea ice data/south/",Sfilenames[i],sep="")
	data=getsouthicedata(fn)
	print(dat[i])
	#Limit data/land mask to all info <250 - 100% and >25 - 10%
	datamask=data<251# & data>25

	#store sum in trend value
	area[ydayindex[i],yearindex[i]]=sum(data[datamask])/250*625
	offsetday[ydayindex[i],yearindex[i]]=oday[i]
	}

	areatimeseries=area
	dim(areatimeseries)=c(366*35,1)
	fareatimeseries=rep(NA,length(areatimeseries))

	######filter timeseries for anomaly calculation only
	######Use more recent years to avoid missing values from early record
	fwid=10
	for (i in 4000:length(areatimeseries))
	{
	if(!is.na(areatimeseries[i]))
	{
	miave=i-10
	maave=i+10

	fareatimeseries[i]=mean(areatimeseries[miave:maave],na.rm=TRUE)
	}
	}

	#resequence all matrix data into time series to complete correction for leap years
	dim(fareatimeseries)=c(366,35)
	anomval=rowMeans(fareatimeseries,na.rm=TRUE)
	SIAnom=area-anomval
	SIArea=area
	dim(SIAnom)=c(366*35,1)
	dim(SIArea)=c(366*35,1)
	dim(offsetday)=c(366*35,1)
	SouthIceArea=rep(NA,366*35)
	SouthIceAreaAnomaly=rep(NA,366*35)

	#remove two extreme values which are obviously errors
	mask= SIAnom< -5000000 \| SIAnom>3000000
	SIAnom[mask]=NA
	SIArea[mask]=NA

	for (i in 1:length(SouthIceArea))
	{
	if(!is.na(offsetday[i]))
	{
	SouthIceArea[offsetday[i]]=SIArea[i]
	SouthIceAreaAnomaly[offsetday[i]]=SIAnom[i]

	}
	}



	SouthIceArea=ts(SouthIceArea,start=1978,deltat=1/365.25) #set time series to get time values only
	decdate=time(SouthIceArea) #copy time values
	SouthIceArea=approx(decdate,SouthIceArea,xout=decdate) #interpolate missing values
	SouthIceArea=ts(SouthIceArea$y,start=1978,deltat=1/365.25) #reassemble series
	SouthIceAreaAnomaly=approx(decdate,SouthIceAreaAnomaly,xout=decdate) #interpolate missing values
	SouthIceAreaAnomaly=ts(SouthIceAreaAnomaly$y,start=1978,deltat=1/365.25) #reassemble series
	SouthIceAreaAnomaly=SouthIceAreaAnomaly-mean(SouthIceAreaAnomaly,na.rm=TRUE) #center series

	plotTrend(SouthIceAreaAnomaly,y.pos=-2e6,main.t = "Antarctic Ice Area Anomaly")



	GlobalIceArea=NorthIceArea+SouthIceArea
	mean(GlobalIceArea,na.rm=TRUE)#19020751
	GlobalIceAreaAnomaly= NorthIceAreaAnomaly+SouthIceAreaAnomaly
	plotTrend(GlobalIceAreaAnomaly,y.pos=-2e6,main.t = "Global Ice Area Anomaly")

	plot(GlobalIceAreaAnomaly+19000000,ylim=c(0,22000000),main="Global Sea Ice Area Anomaly\nOffset by Average Global Ice 19M Km^2", ylab="Km^2 Sea Ice")
	grid()