Created
February 4, 2012 16:26
-
-
Save jstults/1738795 to your computer and use it in GitHub Desktop.
Comparison of Northern Hemisphere Perennial and Seasonal Sea Ice
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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