##SEAICE FUNCTIONS ### ### get.seaice.monthly - retrieves daily data from ftp://sidads.colorado.edu/DATASETS/NOAA/G02135 ### data in directories: ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/Jun/N_200806_area.txt ### make.anomaly - function that makes monthly anomaly relative to all monthly measurements ### ### plot.seaice - twopanel and three panel plots ### ### functions to retrieve NOAA and NSIDC that need summarization ############################################################################ ### NOAA.G02135 ############################################################################ numbertochar=function(x) {if(x<10) paste ("0",x,sep="") else paste(x) } get.seaice.monthly=function(method="NOAA.G02135") { url="ftp://sidads.colorado.edu/DATASETS/NOAA/G02135" # month1=as.character(1:12); y=nchar(month1);temp=(y<2);month1[temp]=paste("0",month1[temp],sep="") month0=c("Jan", "Feb", "Mar", "Apr" ,"May" ,"Jun", "Jul", "Aug" ,"Sep", "Oct" ,"Nov", "Dec") hemi=c("N","S") seaice=NULL for(j in 1:12) { #month for (k in 1:2) { #hemi loc=file.path(url, month0[j],paste(hemi[k],"_",month1[j],"_area.txt",sep="") ) #construct location fred=readLines(loc) year=as.numeric(substr(fred,1,4)) temp1=!is.na(year)&(year>=1978) fred=fred[temp1];year=year[temp1] test=data.frame(year, as.numeric(substr(fred,6,7)),substr(fred,27,27), as.numeric(substr(fred,29,34)),as.numeric(substr(fred,36,41)) ) names(test)=c("year","month","N","area","extent") seaice=rbind(seaice,test) } #k } #j # dim(seaice) #724 5 # 352 5 N=nrow(seaice) order1=order(seaice$N,seaice$year,seaice$month) seaice=seaice[order1,] temp=(seaice$N=="N") X=seaice[temp,] names(X)[4:5]=paste(names(X)[4:5],"N",sep="_") X$area_S=seaice$area[!temp] X$extent_S=seaice$extent[!temp] X$time=round(X$year+(X$month-1)/12,2) N=nrow(X) # X[1,] # year month N area_N extent_N area_S extent_S time #589 1978 11 N 12.02 8.95 16.4 11.56 1978.83 time0=ts(seq(X$year[1]+(X$month[1]-1)/12, X$year[N]+(X$month[N]-1)/12,1/12), ,start=c(X$year[1],X$month[1]),end=c(X$year[N],X$month[N]),freq=12) seaice=data.frame(round(time0,2));names(seaice)="time" test=match(seaice$time,X$time) temp=!is.na(test);sum(!temp) #2 seaice$year=floor(seaice$time); seaice$month= 1+round( 12*(time0-seaice$year)) seaice$area_N=NA;seaice$extent_N=NA;seaice$area_S=NA;seaice$extent_S=NA seaice[temp,4:7] = X[test[temp],4:7] #this copes with a couple of missing values in 1987-1988 seaice$area_GLB=seaice$area_N+seaice$area_S; seaice$extent_GLB=seaice$extent_S+seaice$extent_N return(seaice) } ##seaice=get.seaice.monthly() #Function to extract anomaly make.anomaly=function(x) { #function to make anomaly using all data m0=tapply(x,seaice$month,mean,na.rm=TRUE) y=factor(seaice$month);levels(y)=m0 y=as.numeric(as.character(y)) #monthly average make.anomaly=round(x-y,2) make.anomaly } #Function to extract anomaly make.anomaly=function(x,reference=NULL ) { #function to make anomaly using all data N=nrow(seaice) if( is.null(reference) ) temp=rep(TRUE,N) else temp= (seaice$year>=reference[1])&(seaice$year<(reference[2]+1)) # this is per UIUC jpg m0=tapply(x[temp],seaice$month[temp],mean,na.rm=TRUE) y=factor(seaice$month);levels(y)=m0 y=as.numeric(as.character(y)) #monthly average make.anomaly=round(x-y,2) make.anomaly } plot.seaice=function(seaice,case="extent",GDD_option=FALSE,plot_type="twopanel",reference=c(1978,2000) ) { N=nrow(seaice) url="ftp://sidads.colorado.edu/DATASETS/NOAA/G02135" #source of data month0=c("Jan", "Feb", "Mar", "Apr" ,"May" ,"Jun", "Jul", "Aug" ,"Sep", "Oct" ,"Nov", "Dec") monthto=paste(month0[ seaice$month[N]], seaice$year[N]) name0=dimnames(seaice)[[2]];name0 g=function(x) make.anomaly(x,reference) X=apply(seaice[,grep(case,name0)],2,g) X=data.frame(X) dimnames(X)[[2]]=c("anom_NH","anom_SH","anom_GLB") X$time=seaice$time Trend=list() if(plot_type=="twopanel") { if(GDD_option) GDD(file=paste("d:/climate/images/2009/seaice/seaice_",plot_type,".",seaice$year[N],"_",seaice$month[N],".gif",sep=""),type="gif",w=480,h=480) nf=layout(array(1:2,dim=c(2,1)),heights=c(1.1,1.3)) par(mar=c(0,3,2,1)) plot(X$time,X$anom_SH,type="l",ylab="SH",xlab="",axes=FALSE) axis(side=1,labels=FALSE);axis(side=2,las=1);box();abline(h=0,lty=2) title(paste("Sea Ice Anomaly-",case,monthto)) text(1979,1.2,"SH",pos=4,font=2) Trend$SH=fm=lm(anom_SH~time,data=X);summary(fm) x=predict(fm,newdata=X) #lines(c(seaice$time),x,col=2) par(mar=c(3.5,3,0,1)) plot(X$time,X$anom_NH,type="l",ylab="NH",xlab="",axes=FALSE) axis(side=1,labels=T);axis(side=2,las=1);box();abline(h=0,lty=2) text(1979,-1.0,"NH",pos=4,font=2) mtext(side=1,paste("Source: ",url),line=2,cex=.7) Trend$NH=fm=lm(anom_NH~time,data=X);summary(fm) # -0.038527 0.003445 -11.18 <2e-16 *** x=predict(fm,newdata=X) #lines(c(seaice$time),x,col=2) if(GDD_option) dev.off() } #plot_type if(plot_type=="threepanel") { if(GDD_option) GDD(file=paste("d:/climate/images/2009/seaice/seaice_",plot_type,".",seaice$year[N],"_",seaice$month[N],".gif",sep=""),type="gif",w=480,h=640) nf=layout(array(1:3,dim=c(3,1)),heights=c(1.2,1.,1.3)) par(mar=c(0,3,2,1)) plot(X$time,X$anom_GLB,type="l",ylab="GLB",xlab="",axes=FALSE) axis(side=1,labels=FALSE);axis(side=2,las=1);box();abline(h=0,lty=2) title(paste("Sea Ice Anomaly-",case,monthto)) text(1979,1.5,"GLB",pos=4,font=2) Trend$GLB=fm=lm(anom_GLB~time,data=X);summary(fm) # -0.038527 0.003445 -11.18 <2e-16 *** x=predict(fm,newdata=X) lines(c(X$time),x,col=2) par(mar=c(0,3,0,1)) plot(X$time,X$anom_SH,type="l",ylab="SH",xlab="",axes=FALSE) axis(side=1,labels=FALSE);axis(side=2,las=1);box();abline(h=0,lty=2) text(1979,1.2,"SH",pos=4,font=2) Trend$SH=fm=lm(anom_SH~time,data=X);summary(fm) x=predict(fm,newdata=X) lines(c(X$time),x,col=2) par(mar=c(3.5,3,0,1)) plot(X$time,X$anom_NH,type="l",ylab="NH",xlab="",axes=FALSE) axis(side=1,labels=T);axis(side=2,las=1);box();abline(h=0,lty=2) text(1979,-1.0,"NH",pos=4,font=2) mtext(side=1,paste("Source: ",url),line=2,cex=.7) Trend$NH=fm=lm(anom_NH~time,data=X);summary(fm) # -0.038527 0.003445 -11.18 <2e-16 *** x=predict(fm,newdata=X) lines(c(X$time),x,col=2) if(GDD_option) dev.off() } #plot_type return(Trend) } #function #plot.seaice(seaice,case="extent",GDD_option=FALSE,plot_type="threepanel") ################################### ### JAXA ######################################### #AMSR Ehttp://www.ijis.iarc.uaf.edu/en/home/seaice_extent.htm #The sea-ice extent is calculated as the areal sum of sea ice covering the ocean where sea-ice concentration (SIC) exceeds 15%. SIC data of JAXA’s AMSR-E standard products are used for this purpose (http://sharaku.eorc.jaxa.jp/AMSR/products/pdf/alg_des.pdf). #The area of sea-ice cover is often defined in two ways, i.e., sea-ice “extent” and sea-ice “area.” #These multiple definitions of sea-ice cover may sometimes confuse data users. The former is defined #as the areal sum of sea ice covering the ocean (sea ice + open ocean), whereas #the latter “area” definition counts only sea ice covering a fraction of the ocean (sea ice only). #Thus, the sea-ice extent is always larger than the sea-ice area. Because of the possible errors in SIC mentioned above, # satellite-derived sea-ice concentration can be underestimated, particularly in summer. #In such a case, the sea-ice area is more susceptible to errors than the sea-ice extent. Thus, we adopt the #definition of sea-ice extent to monitor the variation of the Arctic sea ice on this site. get.jaxa=function(method="extent") { url="http://www.ijis.iarc.uaf.edu/seaice/extent/plot.csv" #AMSR-E sea ice concentration algorithm developed by Dr. Comiso in NASA/GSFC. daily=read.csv(url,header=FALSE) #dim(daily) #[1] 2406 4 names(daily)=c("month","day","year","ice") daily$ice[daily$ice== -9999]=NA daily$date<-as.Date(paste(daily$year,daily$month,daily$day,sep="-")) daily$julian<-julian(daily$date,start="1970-01-01") daily$mm=100*daily$year+daily$month daily$dd=factor(daily$year) m0=tapply(daily$julian,daily$year,min) m0[1]=julian(as.Date("2002-01-01"),start="1970-01-01") levels(daily$dd)=m0 daily$dd=as.numeric(as.character(daily$dd)) daily$dd=daily$julian-daily$dd+1 daily$ice=daily$ice/1E6 daily$diff=c(NA,diff(daily$ice)) f=approxfun(daily$julian,daily$ice) daily$ice=f(daily$julian) daily$diff=c(NA,diff(daily$ice)) #month day year ice date julian mm dd diff #2564 6 7 2009 11.028906 2009-06-07 14402 200906 158 -0.05781300 #2565 6 8 2009 11.004688 2009-06-08 14403 200906 159 -0.02421800 daily=daily[!is.na(daily$ice),] return(daily) } #daily=get.jaxa() plot.jaxa=function(daily,GDD_option=FALSE) { url="http://www.ijis.iarc.uaf.edu/seaice/extent/plot.csv" layout(1) N=nrow(daily) if(GDD_option) GDD(file=paste("d:/climate/images/2009/seaice/seaice_jaxa.",daily$date[N],".gif",sep=""),type="gif",w=480,h=360) col0=c(8:3,1,2) temp= daily$year==2008 par(mar=c(3.5,4,2,1)) plot(daily$dd[temp],daily$ice[temp],type="n",ylim=c(4,15),xlab="Julian Day",ylab="million sq km") for ( year in 2002:2009) { temp= daily$year==year lines(daily$dd[temp],daily$ice[temp],col=col0[year-2001] ) } temp= daily$year==2009 lines(daily$dd[temp],daily$ice[temp],col=2,lwd=3) mtext(side=1,paste("Source: ",url),line=2,cex=.7) title("Arctic Seaice Extent (JAXA)") if(GDD_option) dev.off() } # plot.jaxa(daily,GDD_option=FALSE) ################################ ## NSIDC FUNCTIONS ##################################### #the hole is coded as 251; #the hole function counts the number of cells with code 251 #the area function counts cells at 250 and under, thus EXCLUDING the hole; #requires a matrix Cells separately collated for cell areas #multiplies Cell area by percentage cover where pct exceeds 15 #discontinuity when "hole" changes in 1986-87 #the extent function counts cells with code 251 and under, with more than 15% #thus INCLUDES hole at 100% ##COMMENT the extent in 1987 is nearly identical to the extent in 1986 with changed hole #but the area increased by 5139-4584= 555 which is approx amount of discontinuity perhaps # tapply(Arctic$area,Arctic$year,min,na.rm=TRUE) # 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 #4.530 4.644 4.193 4.294 4.491 3.887 4.148 4.584 5.139 4.999 4.698 4.492 4.331 4.906 4.314 4.665 4.277 5.123 4.793 4.124 4.074 # 2000 2001 2002 2003 2004 2005 2006 2007 2008 #4.030 4.405 3.888 3.982 4.146 3.920 3.907 2.793 4.243 # tapply(Arctic$extent,Arctic$year,min,na.rm=TRUE) # 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 #7.844 8.396 7.816 8.009 8.168 7.238 7.457 8.070 8.079 8.193 7.920 7.053 7.300 8.244 7.322 8.057 6.961 8.170 7.432 7.173 6.689 # 2000 2001 2002 2003 2004 2005 2006 2007 2008 #6.934 7.618 6.498 6.875 6.537 6.151 6.774 4.969 7.149 # # tapply(Arctic$hole,Arctic$year,min,na.rm=TRUE) # 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 #1.1926 1.1926 1.1926 1.1926 1.1926 1.1926 1.1926 1.1926 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 # 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 #0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 0.3108 ##SATELLITE INFO make.stat.info=function() { sat.info=array(NA,dim=c(5,3));sat.info=data.frame(sat.info);names(sat.info)=c("sat","start","end") sat.info$sat=c("n07","f08","f11","f13","f15") sat.info$start=c("1979-01-01","1987-07-09","1991-12-03","1995-05-03","2008-01-01") sat.info$end=c("1987-08-20","1991-12-31","1995-09-30","2009-12-31","2009-02-12") sat.info$jstart=julian(as.Date(sat.info$start) ) #1970-01-01 start sat.info$jend=julian(as.Date(sat.info$end) ) #1970-01-01 start return(sat.info) } # sat start end jstart jend #1 n07 1979-01-01 1987-08-20 3287 6440 #2 f08 1987-07-09 1991-12-31 6398 8034 #3 f11 1991-12-03 1995-09-30 8006 9403 #4 f13 1995-05-03 2008-12-31 9253 14244 #5 f15 2008-01-01 2008-12-31 13879 14244 ##CELLS ##POLAR STEORTGTAPH http://nsidc.org/data/grids/ps_grid.html #The polar stereographic projection specifies a projection plane or grid tangent to the Earth at 70 degrees. The planar grid is designed #so that the grid cells at 70 degrees latitude are 25 km by 25 km. #url="ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/tools/pss25area_v2.dat issomething different update.nsidc=function(working=Arctic,end0=Sys.Date() ) { N0=1+max( (1:nrow(working))[!is.na(working$hole)]);N0 N1=(1:nrow(working))[working$julian==end0 ];N1 for (k in N0:N1) { Data=scrape.arctic(day=working$dd[k],month=working$mm[k],year=working$year[k],sat=working$sat[k]) dest.file=paste(working$year[k],working$mm[k],working$dd[k],working$sat[k],sep="_") save(Data,file=file.path("d:/climate/data/seaice/nsidc",paste(dest.file,"tab",sep=".") )) if ( !class(Data)=="logical") { working$hole[k]=hole(Data) working$extent[k]=extent_seaice(Data) working$area[k]= area_seaice(Data) } } #k update.nsidc=working update.nsidc } ## this function needs to be fixed up annual.update=function(year=substr(SysDate(),1,4) ) { N=nrow(Arctic) Arctic[N,] #May 2 2009 # id date year mm dd sat jday area extent extent15 hole julian #10866 2008_09_30_f15.tab 2008-09-30 2008 09 30 f15 274 NA NA NA NA 14152 end0=julian(as.Date(Arctic$date[N])) #14152 dummy=array(NA,dim=c(31+30+31,ncol(Arctic))) dummy=data.frame(dummy); names(dummy)=names(Arctic) dummy[1,] dummy$julian=end0+(1:nrow(dummy)) dummy$date=as.Date(dummy$julian,origin= attributes(end0)$origin) dummy$year=as.numeric(substr(dummy$date,1,4)) dummy$mm=as.numeric(substr(dummy$date,6,7)) dummy$dd=as.numeric(substr(dummy$date,9,10)) dummy$sat="f15" dummy$id= paste( gsub("-","_",dummy$date),"_",dummy$sat,".tab",sep="") dummy$jday=dummy$julian-julian(as.Date("2007-12-31")) Arctic=rbind(Arctic,dummy) } #Arctic=update.nsidc() ##information at http://nsidc.org/data/polar_stereo/ps_grids.html make.Cells=function(hemi="north"){ if(hemi=="north") { url="ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/tools/psn25area_v2.dat" download.file(url,"temp.bin",mode="wb") handle=file("temp.bin","rb") X=readBin(handle, "int", 448*304, size=4, signed = FALSE) close(handle) make.Cells=array(X,dim=c(304,448))/1000 make.Cells=make.Cells[,448:1] } if(hemi=="south") { url="ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/tools/pss25area_v2.dat" download.file(url,"temp.bin",mode="wb") handle=file("temp.bin","rb") X=readBin(handle, "int", 316*332, size=4, signed = FALSE) close(handle) make.Cells=array(X,dim=c(316,332))/1000 make.Cells=make.Cells[,332:1] } make.Cells } Cells=make.Cells() #304 448 #save(Cells,file="d:/climate/data/seaice/Cells.tab") #dim 304 448 #par(mar=c(3,3,2,1)) #image.plot (1:304,1:448,Cells) #Distortion in the grid increases as the latitude decreases, because more of the Earth's surface falls #into any given grid cell, which can be quite significant at the edge of the northern SSM/I grid #where distortion reaches 31 percent. extent_seaice = function(X,method="N") { temp=(X>251);X[temp]=NA temp2=(X>=.15)&!is.na(X) #this includes hole at 100% extent_seaice=sum( Cells[temp2])/1E6 extent_seaice } # extent_seaice(Data) # 10.0122 hole = function(X,method="N") { temp=(X==251); hole=sum( Cells[temp])/1E6 hole } # hole(Data) # ] 0.3107784 area_seaice = function(X,method="N") { temp=(X>250);X[temp]=NA X[!temp]=X[!temp]/250 temp2=(X>=.15)&!is.na(X) area_seaice=sum( (Cells*X)[temp2])/1E6 area_seaice } # area_seaice(Data) # 6.93482 scrape.arctic=function(day,month,year=2009,sat="override") { if(sat=="override") sat=sat_select(DD) if (year>=2008) { loc= "ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/near-real-time/north" ; orig.file=paste("nt_",year,pastedate(month),pastedate(day),"_",sat,"_nrt_n.bin",sep="") url=file.path(loc,orig.file) } else { loc=file.path("ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/final-gsfc/north/daily",year) orig.file=paste("nt_",year,pastedate(month),pastedate(day),"_",sat,"_v01_n.bin",sep="") url=file.path( loc,orig.file)} dest.file=paste(year,month,day,sat,sep="_") test =try( url(url,"rb")) if( length(class(test))==1) scrape.arctic=NA else { download.file(url,"d:/temp/temp.bin",mode="wb") unlink(url); handle=file("d:/temp/temp.bin","rb") header= readBin(handle, "int", 1*300, size=1, signed = FALSE, endian = "big") X = readBin(handle, "int", 448*304, size=1, signed = FALSE, endian = "big") unlink(handle);close(handle) X=array(X,dim=c(304,448)) scrape.arctic=X[,448:1] } scrape.arctic } pastedate=function(x) {if (nchar(x)==1) pastedate=paste("0",x,sep="") else pastedate=paste(x);pastedate} #day=Arctic$dd[k];month=Arctic$mm[k];year=Arctic$year[k];sat=Arctic$sat[k] ####BREMEN # http://iup.physik.uni-bremen.de:8084/amsredata/asi_daygrid_swath/l1a/n6250/README.TXT #monthly: ftp://ftp-projects.zmaw.de/seaice/AMSR-E_ASI_IceConc/area-extent/area_ASI-AMSR-E_v5.5i_Ant.txt