the Air Vent

Because the world needs another opinion

Gridded Daily Sea Ice – NSIDC

Posted by Jeff Id on October 3, 2009

I’ve been working on improving the noise rejection of the sea ice area and anomaly from the NSIDC gridded data. These plots are created by summation of the gridded satellite sea ice concentration data the following table is up to date as of 2008

Platform and Instrument Time Period
Nimbus-7 SMMR 26 October 1978 through 20 August 1987
DMSP-F8 SSM/I 9 July 1987 through 31 December 1991
DMSP-F11 SSM/I 3 December 1991 through 30 September 1995
DMSP-F13 SSM/I 3 May 1995 through the March 26 2008
DMSP-F17 SSM/I 26 Mar 2008 through present

This table is not the official table presented by the NSIDC as I have added the end point for F13 and the beginning point for NOAA-17. I should note that for a time in 2008 overlapping data is presented for both F13 and F17.

My plots still aren’t perfect matches to the NSIDC or Cryosphere but they are very close. Some reasons for differences include:

  • No access to final processing routines used to blend and filter daily trends.
  • The anomaly for these plots is computed over the whole length of the dataset whereas other scientific institutions limit the anomaly calculation period to different fixed sets of years.
  • Filtering differences
  • Differences in methods for locating and removing bad points from the near real time data

Either way, differences in the results are quite small and the net trends and absolute magnitude are close. Thanks to the NSIDC for making the data available, their openness gives me an improved comfort in the quality and completeness of the work done by this agency.

Antarctic Area

Antarctic Extent

Arctic Area

Arctic Extent

Global Area

Global Extent

Global Ice anomaly offset

The code for creating these graphs is below. It requires you to download a huge pile of sea ice data to your local harddrive which if I recall took about 15 hours per hemisphere so I doubt anyone else will try but for documentation purposes I’ve included it here (reformatting courtesy wordpress).

filenames=list.files(path=”C:/agw/sea ice/north sea ice/nasateam daily/”, pattern = NULL, all.files = TRUE, full.names = FALSE, recursive = TRUE)
trend=array(0,dim=length(filenames)-1)
extent=array(0,dim=length(filenames)-1)
date=array(0,dim=length(filenames)-1)
masktrend=array(0,dim=length(filenames)-1)

for(i in 1:(length(filenames)-1))
{
#open filename
fn=paste(“C:/agw/sea ice/north sea ice/nasateam daily/”,filenames[i],sep=””)
a=file(fn,”rb”)

#extract header info -not used portion
header= readBin(a,n=102,what=integer(),size=1,endian=”little”,signed=FALSE)

#assemble year and day
year=readChar(a,n=6)
print(year)
day=readChar(a,n=6)
print(day)
header=readChar(a,n=300-114)#finish loading 300 char header

#get pixel data – little endian form
data=readBin(a,n=304*448,what=integer(),size=1,endian=”little”,signed=FALSE)
close(a)

#fix date format to be year plus day expressed as decimal of 365 day year – ignore leap year
if(as.integer(year)+1900<=2500)
{
date[i]=1900+as.integer(year)+as.integer(day)/365
}else
{
date[i]=as.integer(year)+as.integer(day)/365
}

#use hole mask for early satellite in entire trend. Insures no steps in trend
if(i==1)
{
holemask= !(data==251)
}

#Limit data/land mask to all info <250 – 100% and >37 – 15%
datamask=data<251# & data>37

#store sum in trend value
trend[i]=sum(data[(datamask*holemask)==1])/250*625
extentmask=data[(datamask*holemask)==1]>37
extent[i]=sum(extentmask)
}
#create mask for NOAA 15 values
#satname=substring(filenames,18,20)
#satmask= satname==”f15″

########################################
#insert data into equal spaced timeseries
areatrend=rep(NA,(2010-1978)*365)
extenttrend=rep(NA,(2010-1978)*365)

for (i in 1:length(data))
{
index=(date[i]-1978)*365+.5
##1799*625 is satellite hole mask
areatrend[index]=trend[i]+1799*625
extenttrend[index]=extent[i]*625+1799*625
}

fareatrend=array(0,dim=length(areatrend))
fextenttrend=array(0,dim=length(areatrend))

#filter short for removal of bad data points
for(i in 1:(length(areatrend)))
{
sumdat=0
sumdatb=0
count=0
for(j in -6:6)
{
k=i+j
if(k<1)k=1
if(k>length(fareatrend)-1)k=length(fareatrend)-1
if(!is.na(areatrend[k]))
{
sumdat=sumdat+areatrend[k]
sumdatb=sumdatb+extenttrend[k]
count=count+1
}
}
if(count>0)
{
fareatrend[i]=sumdat/count
fextenttrend[i]=sumdatb/count

}else
{
fareatrend[i]=NA
fextenttrend[i]=NA
}
}

tt=fareatrend-areatrend
mask=abs(tt)>250000

mareatrend=areatrend
mareatrend[mask]=NA
mextenttrend=extenttrend
mextenttrend[mask]=NA

#filter long for result
fareatrend=array(0,dim=length(areatrend))
fextenttrend=array(0,dim=length(areatrend))

for(i in 1:(length(areatrend)))
{
sumdat=0
sumdatb=0
count=0
for(j in -6:6)
{
k=i+j
if(k<1)k=1
if(k>length(fareatrend)-1)k=length(fareatrend)-1
if(!is.na(mareatrend[k]))
{
sumdat=sumdat+mareatrend[k]
sumdatb=sumdatb+mextenttrend[k]
count=count+1
}
}
if(count>0)
{
fareatrend[i]=sumdat/count
fextenttrend[i]=sumdatb/count

}else
{
fareatrend[i]=NA
fextenttrend[i]=NA
}
}

NorthIceArea=fareatrend
NorthIceExtent=fextenttrend

##############################################
##############################################
#calc daily anomaly ignore leap years

dim(NorthIceArea)=c(365,32)
an=rowMeans(NorthIceArea,na.rm=TRUE)
NorthIceAreaAnomaly = NorthIceArea-an
dim(NorthIceArea)=365*32
dim(NorthIceAreaAnomaly)=365*32
plot(NorthIceAreaAnomaly)

dim(NorthIceExtent)=c(365,32)
an=rowMeans(NorthIceExtent,na.rm=TRUE)
NorthIceExtentAnomaly = NorthIceExtent-an
dim(NorthIceExtent)=365*32
dim(NorthIceExtentAnomaly)=365*32
plot(NorthIceExtentAnomaly)

NorthIceArea=ts(NorthIceArea,start=1978,deltat=1/365)
NorthIceAreaAnomaly=ts(NorthIceAreaAnomaly,start=1978,deltat=1/365)
NorthIceExtent=ts(NorthIceExtent,start=1978,deltat=1/365)
NorthIceExtentAnomaly=ts(NorthIceExtentAnomaly,start=1978,deltat=1/365)

##########################################
##########################################
#southern hemisphere sea ice
filenames=list.files(path=”C:/agw/sea ice/south sea ice/nasateam daily/”, pattern = NULL, all.files = TRUE, full.names = FALSE, recursive = TRUE)
trend=array(0,dim=length(filenames)-1)
extent=array(0,dim=length(filenames)-1)
date=array(0,dim=length(filenames)-1)
masktrend=array(0,dim=length(filenames)-1)

for(i in 1:(length(filenames)-1))
{
fn=paste(“C:/agw/sea ice/south sea ice/nasateam daily/”,filenames[i],sep=””)

a=file(fn,”rb”)
header= readBin(a,n=102,what=integer(),size=1,endian=”little”,signed=FALSE)
year=readChar(a,n=6)
print(year)
day=readChar(a,n=6)
print(day)
header=readChar(a,n=300-114)
data=readBin(a,n=304*448,what=integer(),size=1,endian=”little”,signed=FALSE)
close(a)

if(as.integer(year)+1900<=2500)
{
date[i]=1900+as.integer(year)+as.integer(day)/365
}else
{
date[i]=as.integer(year)+as.integer(day)/365
}
if(i==1)
{
holemask= !(data==251)
}
#Limit data/land mask to all info <250 – 100% and >37 – 15%
datamask=data<251# & data>37

#store sum in trend value
trend[i]=sum(data[(datamask*holemask)==1])/250*625
extentmask=data[(datamask*holemask)==1]>37
extent[i]=sum(extentmask)
}

#satname=substring(filenames,18,20)
#satmask= satname==”f15″

########################################
#insert data into equal spaced timeseries
areatrend=rep(NA,(2010-1978)*365)
extenttrend=rep(NA,(2010-1978)*365)

for (i in 1:length(data))
{
index=(date[i]-1978)*365+.5
##1799*625 is satellite hole mask
areatrend[index]=trend[i]+1799*625
extenttrend[index]=extent[i]*625+1799*625
}

fareatrend=array(0,dim=length(areatrend))
fextenttrend=array(0,dim=length(areatrend))

#filter short for removal of bad data points
for(i in 1:(length(areatrend)))
{
sumdat=0
sumdatb=0
count=0
for(j in -6:6)
{
k=i+j
if(k<1)k=1
if(k>length(fareatrend)-1)k=length(fareatrend)-1
if(!is.na(areatrend[k]))
{
sumdat=sumdat+areatrend[k]
sumdatb=sumdatb+extenttrend[k]
count=count+1
}
}
if(count>0)
{
fareatrend[i]=sumdat/count
fextenttrend[i]=sumdatb/count

}else
{
fareatrend[i]=NA
fextenttrend[i]=NA
}
}

tt=fareatrend-areatrend
mask=abs(tt)>600000

mareatrend=areatrend
mareatrend[mask]=NA
mextenttrend=extenttrend
mextenttrend[mask]=NA

#filter long for result
fareatrend=array(0,dim=length(areatrend))
fextenttrend=array(0,dim=length(areatrend))

for(i in 1:(length(areatrend)))
{
sumdat=0
sumdatb=0
count=0
for(j in -6:6)
{
k=i+j
if(k<1)k=1
if(k>length(fareatrend)-1)k=length(fareatrend)-1
if(!is.na(mareatrend[k]))
{
sumdat=sumdat+mareatrend[k]
sumdatb=sumdatb+mextenttrend[k]
count=count+1
}
}
if(count>0)
{
fareatrend[i]=sumdat/count
fextenttrend[i]=sumdatb/count

}else
{
fareatrend[i]=NA
fextenttrend[i]=NA
}
}

SouthIceArea=fareatrend
SouthIceExtent=fextenttrend

##############################################
##############################################
#calc daily anomaly ignore leap years

dim(SouthIceArea)=c(365,32)
an=rowMeans(SouthIceArea,na.rm=TRUE)
SouthIceAreaAnomaly = SouthIceArea-an
dim(SouthIceArea)=365*32
dim(SouthIceAreaAnomaly)=365*32
plot(SouthIceAreaAnomaly)

dim(SouthIceExtent)=c(365,32)
an=rowMeans(SouthIceExtent,na.rm=TRUE)
SouthIceExtentAnomaly = SouthIceExtent-an
dim(SouthIceExtent)=365*32
dim(SouthIceExtentAnomaly)=365*32
plot(SouthIceExtentAnomaly)

SouthIceArea=ts(SouthIceArea,start=1978,deltat=1/365)
SouthIceAreaAnomaly=ts(SouthIceAreaAnomaly,start=1978,deltat=1/365)
SouthIceExtent=ts(SouthIceExtent,start=1978,deltat=1/365)
SouthIceExtentAnomaly=ts(SouthIceExtentAnomaly,start=1978,deltat=1/365)

################################
################################
################################
#compute global values
GA=ts.union(NorthIceAreaAnomaly,SouthIceAreaAnomaly)
plot(GA,main=”Hemispheric Sea Ice Anomalies\nTop = Arctic, Bottom = Antarctic”,cex.main=1,xlab=”Year” )

GlobalIceArea=SouthIceArea+NorthIceArea
plot(GlobalIceArea)
GlobalIceAreaAnomaly=SouthIceAreaAnomaly+NorthIceAreaAnomaly
plot(GlobalIceAreaAnomaly)
GlobalIceExtent=SouthIceExtent+NorthIceExtent
plot(GlobalIceExtent)
GlobalIceExtentAnomaly=SouthIceExtentAnomaly+NorthIceExtentAnomaly
plot(GlobalIceExtentAnomaly)

mean(GlobalIceArea,na.rm=TRUE)#[1] 20194725

###########################################
#North plots
#par(oma=c(1,0,0,0), mar=c(2,4,3,1),mfrow=c(2,1), bg=”grey80″)

plot(NorthIceArea,main=”Arctic Ice Area”,xlab=”Year”,ylab=”Km^2″,ylim=c(0,1.9e7))
grid(col=”black”)
abline(lsfit(time(NorthIceArea),NorthIceArea,NULL), col=”red”,lwd=2)
text(1985,1e6,paste(“Slp =”,round( coef(lsfit(time(NorthIceArea),NorthIceArea,NULL))[2],0) ),col=”red”)
#savePlot(paste(“c:/agw/sea ice/North Ice area.jpg”),type=”jpg”)

plot(NorthIceAreaAnomaly,main=”Arctic Ice Area Anomaly”,xlab=”Year”,ylab=”Km^2″)
grid(col=”black”)
abline(lsfit(time(NorthIceAreaAnomaly),NorthIceAreaAnomaly,NULL), col=”red”,lwd=2)
text(1985,-2e6,paste(“Slp =”,round( coef(lsfit(time(NorthIceAreaAnomaly),NorthIceAreaAnomaly,NULL))[2],0) ),col=”red”)
#savePlot(paste(“c:/agw/sea ice/North Ice area anomaly.jpg”),type=”jpg”)

#savePlot(paste(“c:/agw/sea ice/Arctic Area.jpg”),type=”jpg”)

plot(NorthIceExtent,main=”Arctic Ice Etent”,xlab=”Year”,ylab=”Km^2″,ylim=c(0,1.9e7))
grid(col=”black”)
abline(lsfit(time(NorthIceExtent),NorthIceExtent,NULL), col=”red”,lwd=2)
text(1985,1e6,paste(“Slp =”,round( coef(lsfit(time(NorthIceExtent),NorthIceExtent,NULL))[2],0) ),col=”red”)
#savePlot(paste(“c:/agw/sea ice/North Ice Extent.jpg”),type=”jpg”)

plot(NorthIceExtentAnomaly,main=”Arctic Ice Extent Anomaly”,xlab=”Year”,ylab=”Km^2″)
grid(col=”black”)
abline(lsfit(time(NorthIceExtentAnomaly),NorthIceExtentAnomaly,NULL), col=”red”,lwd=2)
text(1985,-2e6,paste(“Slp =”,round( coef(lsfit(time(NorthIceExtentAnomaly),NorthIceExtentAnomaly,NULL))[2],0) ),col=”red”)
#savePlot(paste(“c:/agw/sea ice/North Ice Extent anomaly.jpg”),type=”jpg”)
#savePlot(paste(“c:/agw/sea ice/Arctic Extent.jpg”),type=”jpg”)

###########################################
#South plots

plot(SouthIceArea,main=”Antarctic Ice Area”,xlab=”Year”,ylab=”Km^2″,ylim=c(0,2.2e7))
grid(col=”black”)
abline(lsfit(time(SouthIceArea),SouthIceArea,NULL), col=”red”,lwd=2)
text(1985,1e6,paste(“Slp =”,round( coef(lsfit(time(SouthIceArea),SouthIceArea,NULL))[2],0) ),col=”red”)
#savePlot(paste(“c:/agw/sea ice/South Ice Area.jpg”),type=”jpg”)

plot(SouthIceAreaAnomaly,main=”Antarctic Ice Area Anomaly”,xlab=”Year”,ylab=”Km^2″)
grid(col=”black”)
abline(lsfit(time(SouthIceAreaAnomaly),SouthIceAreaAnomaly,NULL), col=”red”,lwd=2)
text(1985,1.5e6,paste(“Slp =”,round( coef(lsfit(time(SouthIceAreaAnomaly),SouthIceAreaAnomaly,NULL))[2],0) ),col=”red”)
#savePlot(paste(“c:/agw/sea ice/South Ice Area anomaly.jpg”),type=”jpg”)
#savePlot(paste(“c:/agw/sea ice/Antarctic Area.jpg”),type=”jpg”)

plot(SouthIceExtent,main=”Antarctic Ice Etent”,xlab=”Year”,ylab=”Km^2″,ylim=c(0,2.2e7))
grid(col=”black”)
abline(lsfit(time(SouthIceExtent),SouthIceExtent,NULL), col=”red”,lwd=2)
text(1985,1e6,paste(“Slp =”,round( coef(lsfit(time(SouthIceExtent),SouthIceExtent,NULL))[2],0) ),col=”red”)
#savePlot(paste(“c:/agw/sea ice/South Ice Extent.jpg”),type=”jpg”)

plot(SouthIceExtentAnomaly,main=”Antarctic Ice Extent Anomaly”,xlab=”Year”,ylab=”Km^2″)
grid(col=”black”)
abline(lsfit(time(SouthIceExtentAnomaly),SouthIceExtentAnomaly,NULL), col=”red”,lwd=2)
text(1985,1.5e6,paste(“Slp =”,round( coef(lsfit(time(SouthIceExtentAnomaly),SouthIceExtentAnomaly,NULL))[2],0) ),col=”red”)
#savePlot(paste(“c:/agw/sea ice/North Ice extent anomaly.jpg”),type=”jpg”)
#savePlot(paste(“c:/agw/sea ice/Antarctic Extent.jpg”),type=”jpg”)

###########################################
#Global plots

plot(GlobalIceArea,main=”Global Sea Ice Area “,xlab=”Year”,ylab=” Km^2″,ylim=c(0,2.4e7))
grid(col=”black”)
abline(lsfit(time(GlobalIceArea),GlobalIceArea,NULL), col=”red”,lwd=2)
text(1985,1000000,paste(“Slp =”,round( coef(lsfit(time(GlobalIceArea),GlobalIceArea,NULL))[2],0) ),col=”red”)
#savePlot(paste(“c:/agw/sea ice/Global Ice Area.jpg”),type=”jpg”)

plot(GlobalIceAreaAnomaly,main=”Global Sea Ice Area Anomaly”,xlab=”Year”,ylab=”Anomaly Km^2″)
grid(col=”black”)
abline(lsfit(time(GlobalIceAreaAnomaly),GlobalIceAreaAnomaly,NULL), col=”red”,lwd=2)
text(1985,-2e6,paste(“Slp =”,round( coef(lsfit(time(GlobalIceAreaAnomaly),GlobalIceAreaAnomaly,NULL))[2],0) ),col=”red”)
#savePlot(paste(“c:/agw/sea ice/Global Ice Area anomaly.jpg”),type=”jpg”)
#savePlot(paste(“c:/agw/sea ice/Global Area.jpg”),type=”jpg”)

plot(GlobalIceExtent,main=”Global Sea Ice Extent “,xlab=”Year”,ylab=” Km^2″,ylim=c(0,3e7))
grid(col=”black”)
abline(lsfit(time(GlobalIceExtent),GlobalIceExtent,NULL), col=”red”,lwd=2)
text(1985,1e6,paste(“Slp =”,round( coef(lsfit(time(GlobalIceExtent),GlobalIceExtent,NULL))[2],0) ),col=”red”)
#savePlot(paste(“c:/agw/sea ice/Global Ice Extent.jpg”),type=”jpg”)

plot(GlobalIceExtentAnomaly,main=”Global Sea Ice Extent Anomaly”,xlab=”Year”,ylab=”Anomaly Km^2″)
grid(col=”black”)
abline(lsfit(time(GlobalIceExtentAnomaly),GlobalIceExtentAnomaly,NULL), col=”red”,lwd=2)
text(1985,-2e6,paste(“Slp =”,round( coef(lsfit(time(GlobalIceAreaAnomaly),GlobalIceAreaAnomaly,NULL))[2],0) ),col=”red”)
#savePlot(paste(“c:/agw/sea ice/Global Ice Extent anomaly.jpg”),type=”jpg”)
#savePlot(paste(“c:/agw/sea ice/Global Extent.jpg”),type=”jpg”)

GSA=GlobalIceAreaAnomaly+mean(GlobalIceArea,na.rm=TRUE)
plot(GSA,main=”Global Sea Ice Area Anomaly\nAverage Area Added Back In”,xlab=”Year”,ylab=”Anomaly Km^2″,ylim=c(0,3e7))
grid(col=”black”)
abline(lsfit(time(GSA),GSA,NULL), col=”red”,lwd=2)
text(1985,2e6,paste(“Slp =”,round( coef(lsfit(time(GSA),GSA,NULL))[2],0) ),col=”red”)
#savePlot(paste(“c:/agw/sea ice/Global Ice Area anomaly offset.jpg”),type=”jpg”)

GSE=GlobalIceExtentAnomaly+mean(GlobalIceExtent,na.rm=TRUE)
plot(GSE,main=”Global Sea Ice Extent Anomaly\nAverage Exent Added Back In”,xlab=”Year”,ylab=”Anomaly Km^2″,ylim=c(0,3e7))
grid(col=”black”)
abline(lsfit(time(GSE),GSE,NULL), col=”red”,lwd=2)
text(1985,2e6,paste(“Slp =”,round( coef(lsfit(time(GSE),GSE,NULL))[2],0) ),col=”red”)
#savePlot(paste(“c:/agw/sea ice/Global Ice Extent anomaly offset.jpg”),type=”jpg”)
#savePlot(paste(“c:/agw/sea ice/Global Ice anomaly offset.jpg”),type=”jpg”)


6 Responses to “Gridded Daily Sea Ice – NSIDC”

  1. DeWitt Payne said

    Some sort of moving average seems to be used by most reporting agencies. JAXA, for example, uses a two day average. My preference would be a three or five day centered average on a pixel by pixel basis. No data pixels would be ignored when constructing the average. Using three days should mean the probability of no data for a pixel on all three days would be very small.

  2. Jeff, the work in this is utterly astounding.

    Now the requests. What are the y-axis units? What is slp? and finally, What is the meaning and significance of “Average area/extent added back in”?

  3. Jeff Id said

    Thanks,

    Kilometers squared, Slope and the anomaly is offset from zero to show how fast we’re running out of sea ice.

  4. Thanks. ah, I see now, e+o7 means 10 to the power of 7 ie 10 million. So if you labelled y-axis as “millions of square km” the markers would be 10,20,30… And the last graphs are just to show the sea ice changes are a “storm in a teacup”, with fluctuations well in excess of trend.

    Some kind of mean interannual fluctuation ought to be included alongside trend, IMHO, to get the feeling for polar regions’ huge variability. But of course, why just take one average… when what would be far more significant is to show the presence of cycles that can be correlated with solar cycles.

  5. Charlie said

    You probably know all of the info below, but just in case it slipped by you:

    As for the processing, there is are glitches in the data on June 1 and October 15 as some parameter is changed to account for the reflection of melt water on top of sea ice.

    At http://www.ijis.iarc.uaf.edu/en/home/seaice_extent.htm IARC-JAXA says “The algorithm for calculating SIC was developed and provided by Dr. Comiso of NASA GSFC through a cooperative relationship between NASA and JAXA.”

    Considering the tight relationship between NSIDC and NASA I’d bet that Josefino Comiso’s algorithms are used for the NSIDC data also.

    His contact info is at http://eospso.gsfc.nasa.gov/directory/get_person.php?eosDirID=eos0000301

  6. DJA said

    The Catlin expedition results have emerged see http://www.theaustralian.news.com.au/story/0,25197,26213406-601,00.html

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

 
%d bloggers like this: