A little cleanup of code for gridded sea ice.
Download a free software called FileZilla and copy sea ice files to a known folder.
The following code is a mathematically clean download and summation of gridded microwave satellited data. For the north the pole-hole (area not observed by the satellites) is infilled to the maximum value for its entire timeframe. For the south, the pole hole is not relevant. Two points of the 12000 plus day record are removed in the south and all points are kept in the north. Results are visually comparable to the UIUC cryosphere page. There isn’t much to improve about this for accuracy except for pole hole infilling or an unknown error.
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
}
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
}
##################################################################################
##################################################################################
##################################################################################
##################################################################################
##################################################################################
##################################################################################
#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
yday=as.POSIXlt(dat)$yday #day of the year. Used for anomaly calculations.
sat=substr(Nfilenames,18,20)
satlist=levels(factor(sat))
farea=rep(NA,length(Nfilenames))
#set up hole mask here
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 year
for(i in 1:(length(Nfilenames)))
{
#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
farea[i]=sum(data[(datamask&holemask)])/250*625+sum(!holemask)*625
}
st=oday[1]
en=oday[length(oday)]
Ifarea = approx(oday,farea,st:en)#interpolated series for anom
plot(Ifarea,type="l")
yd=as.POSIXlt((st:en)*86400, origin="1977-12-31")$yday+1 #day of the year. Used for anomaly calculations.
yr=as.POSIXlt((st:en)*86400, origin="1977-12-31")$year-77 #day of the year. Used for anomaly calculations.
#calculate anomaly
yrs=max(yr)
anomarry=array(NA,dim=c(yrs,366))
for(i in 1:length(yr))
{
anomarry[yr[i],yd[i]]=Ifarea$y[i]
}
an=colMeans(anomarry,na.rm=TRUE)
plot(an)
## North sea ice area, no infilling all satellites
NIA=farea
tm=1978+oday/365.25
## North sea ice area anomaly, no infilling all satellites
NIAA=NIA-an[yday+1]
plot(tm,NIAA,type="l",main="Arctic Ice Area Anomaly\nFrom Gridded data - no infilling",ylab="Area Km^2")
grid()
## Assemble all satellites into timeseries
allNSAT=ts(array(NA,dim=c(max(oday),length(satlist))),start=1978,deltat=1/365.25)
## Assemble continuous area series
mask=diff(tm)==0
maska=c(TRUE,!mask)
NIAclean=NIA[maska]
tmclean=tm[maska]
plot(tmclean,NIAclean,type="l")
##Assemble continuous area anomaly series
NIAAclean=NIAA[maska]
tmclean=tm[maska]
plot(tmclean,NIAAclean,type="l")
##Create timeseries from above
NIAts=ts(approx(tmclean,NIAclean,1978+(st:en)/365.25)$y,start=tmclean[1],deltat=1/365.25)
NIAAts=ts(approx(tmclean,NIAAclean,1978+(st:en)/365.25)$y,start=tmclean[1],deltat=1/365.25)
############################################################################################
##################################################################################
##################################################################################
##################################################################################
##################################################################################
##################################################################################
##################################################################################
#calculate southern sea ice area from gridded satelite data
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
yday=as.POSIXlt(dat)$yday #day of the year. Used for anomaly calculations.
sat=substr(Sfilenames,18,20)
satlist=levels(factor(sat))
farea=rep(NA,length(Sfilenames))
#################################
#load area values into 2d matrix 366days each year x 35 years
for(i in 1:(length(Sfilenames)))
{
#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
farea[i]=sum(data[datamask])/250*625
}
st=oday[1]
en=oday[length(oday)]
##########################################################
#calculate all first pass - two bad points must be removed
Ifarea = approx(oday,farea,st:en)#interpolated series for anom
plot(Ifarea,type="l")
yd=as.POSIXlt((st:en)*86400, origin="1977-12-31")$yday+1 #day of the year. Used for anomaly calculations.
yr=as.POSIXlt((st:en)*86400, origin="1977-12-31")$year-77 #day of the year. Used for anomaly calculations.
yrs=max(yr)
anomarry=array(NA,dim=c(yrs,366))
for(i in 1:length(yr))
{
anomarry[yr[i],yd[i]]=Ifarea$y[i]
}
an=colMeans(anomarry,na.rm=TRUE)
plot(an)
## South sea ice area, no infilling all satellites
SIA=farea
## North sea ice area anomaly, no infilling all satellites
SIAA=SIA-an[yday+1]
mask=!(SIAA>5e6 | SIAA< -5e6)
###########################################################
#calculate all second pass - two bad points must be removed
Ifarea = approx(oday[mask],farea[mask],st:en)#interpolated series for anom
plot(Ifarea,type="l")
yd=as.POSIXlt((st:en)*86400, origin="1977-12-31")$yday+1 #day of the year. Used for anomaly calculations.
yr=as.POSIXlt((st:en)*86400, origin="1977-12-31")$year-77 #day of the year. Used for anomaly calculations.
yrs=max(yr)
anomarry=array(NA,dim=c(yrs,366))
for(i in 1:length(yr))
{
anomarry[yr[i],yd[i]]=Ifarea$y[i]
}
an=colMeans(anomarry,na.rm=TRUE)
plot(an)
## South sea ice area, no infilling all satellites
SIA=farea[mask]
tm=1978+oday[mask]/365.25
## North sea ice area anomaly, no infilling all satellites
SIAA=(farea-an[yday+1])[mask]
##########################################################
plot(tm,SIAA,type="l",main="Arctic Ice Area Anomaly\nFrom Gridded data - no infilling",ylab="Area Km^2")
grid()
## Assemble all satellites into timeseries
allSSAT=ts(array(NA,dim=c(max(oday),length(satlist))),start=1978,deltat=1/365.25)
## Assemble continuous area series
mask=diff(tm)==0
maska=c(TRUE,!mask)
SIAclean=SIA[maska]
tmclean=tm[maska]
plot(tmclean,SIAclean,type="l")
##Assemble continuous area anomaly series
SIAAclean=SIAA[maska]
tmclean=tm[maska]
plot(tmclean,SIAAclean,type="l")
##Create timeseries from above
SIAts=ts(approx(tmclean,SIAclean,1978+(st:en)/365.25)$y,start=tmclean[1],deltat=1/365.25)
SIAAts=ts(approx(tmclean,SIAAclean,1978+(st:en)/365.25)$y,start=tmclean[1],deltat=1/365.25)
############################################################################################
the Air Vent 
Jeff,
You don’t need a FTP program like Filezilla to download the data in Windows.
1. Copy the link
2. Open Windows Explorer
3. Paste the link into the address bar and hit enter
4. Copy the files to your hard drive
Greg,
There are thousands of files in different folders and the GB of download takes hours. Filezilla is much easier than the copy paste method.
Jeff,
I don’t see your point. If I open the FTP site you posted I see 3 folders and a readme file. All I have to do is select the 3 folders and the readme, copy, paste. How could it be any easier? The download won’t take any longer then it would with a FTP client so the size is not relevant. I do use Filezilla for FTP servers that require authentication to manage my accounts.
Greg,
You are right of course. The files took me two days to download and I find Filezilla an easy way to clean up any partials from start/stopping.
Doug, Go away. You are not qualified to comment at tAV.
Jeff,
Why do area? Cryosphere Today data seems to be just fine for area. What’s missing is extent. I’ve been getting extent from MASIE since the detector in the AQUA satellite failed. But MASIE hasn’t updated in two days now, and I’m a little worried. MASIE seems to be something of a shoestring organization inside NOAA and they may have run out of money.
The data gives area fairly cleanly but I think you may be right. It introduces noise which is not well quantified.
I have a guest post in the queue which discusses sea ice trends. I’m not sure whether to run it because it is a overly-skeptical about the quality of the sat data. I think it is worthwhile for discussion though. ESMR in particular, has some ugly information that when taken in context of the multi-channel sensors, gives a strong impression of the difficulties faced by these instruments.