Skip to content

Embed URL

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Comparison of Northern Hemisphere Perennial and Seasonal Sea Ice
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.