#http://hot-topic.co.nz/five-years-threnody-for-arctic-sea-ice/ # -- get data from PIOMAS model f1 = "data/PIOMAS.vol.daily.1979.2012.Current.v2.dat" piodata = read.table(f1,skip=1,col.names=c("year","day","volume")) piosept = piodata[piodata[,2]>=244&piodata[,2]<=273,] r1 = min(piosept[,1]) r2 = max(piosept[,1]) vol = rep(0,r2-r1+1) for (i in r1:r2) { vol[i-r1+1] = ave(piosept[piosept[,1]==i,3])[1] } # -- get NSIDC September extent/area # "ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/Sep/N_09_area.txt" # reformatted f2 = "data/N_09_area.txt" si = read.table(f2,header=TRUE) # thickness = volume / [extent|area] thicke = vol/si[,5] thicka = vol/si[,6] # -- Get 21 yr Trends # 21 year sets (1990-2010) sie1990 = si[12:32,5] sia1990 = si[12:32,6] vol1990 = vol[12:32] thicke1990 = thicke[12:32] thicka1990 = thicka[12:32] # 21 year trends sie1990 = lm(sie1990~c(1990:2010)) sia1990 = lm(sia1990~c(1990:2010)) lmvol1990 = lm(vol1990~c(1990:2010)) lmthicke1990 = lm(thicke1990~c(1990:2010)) lmthicka1990 = lm(thicka1990~c(1990:2010)) # -- Get 10 yr Trends # ten year sets (2001-2010) sie10 = si[23:32,5] sia10 = si[23:32,6] vol10 = vol[23:32] thicke10 = thicke[23:32] thicka10 = thicka[23:32] # 10 trends lmsie10 = lm(sie10~c(2001:2010)) lmsie10 = lm(sia10~c(2001:2010)) lmvol10 = lm(vol10~c(2001:2010)) lmthicke10 = lm(thicke10~c(2001:2010)) lmthicka10 = lm(thicka10~c(2001:2010)) # -- Get 21 year projection via Gareth's method for EXTENT # -- this pegs projection against 2010 data point #thicke21_2020=lmthicke1990[[1]][2]*c(2011:2020)+lmthicke1990[[1]][1] thicke21_2020=lmthicke1990[[1]][2]*c(1:10)+thicke[32] #thicka21_2020=lmthicka1990[[1]][2]*c(2011:2020)+lmthicka1990[[1]][1] #thicka21_2020=lmthicka1990[[1]][2]*c(1:10)+thicka[32] #vol21_2020=lmvol1990[[1]][2]*c(2011:2020)+lmvol1990[[1]][1] vol21_2020=lmvol1990[[1]][2]*c(1:10)+vol[32] sie21_2020 = vol21_2020/thicke21_2020 #sia21_2020 = vol21_2020/thicka21_2020 # - plot graphic = TRUE if (graphic) png("img/threnody-gareth-21.png",w=600,h=400) plot(c(vol[1:32],vol21_2020)~c(1979:2020), type="o", col="blue", pch=19, ylim=c(0,13), xlim=c(2001,2020), xaxs="i",yaxs="i",xaxt="n", xlab="", ylab="") axis(1, at=c(1979:2020)) abline(v="2010",col="dark green",lwd=2) lines(c(si[1:32,5],sie21_2020)~c(1979:2020), type="o", col="green", pch=19) lines(c(thicke[1:32],thicke21_2020)~c(1979:2020), type="o", col="red", pch=19) points(vol[33:34]~c(2011:2012),col="blue",pch=13,cex=2) points(si[33:34,5]~c(2011:2012),col="green",pch=13,cex=2) title(main="Reconstructed Threnody for Sea Ice Extent\n21 year trend projected from 2010") text(2002,12, "Volume (1000 km^3)",adj=c(0,0)) text(2002,7, "Extent (million km^2)",adj=c(0,0)) text(2002,2.5, "Thickness (m)",adj=c(0,0)) mtext(paste("http://rhinohide.wordpress.com",format(Sys.time(), "%Y%m%d")), 4, side=1, adj=1, cex=0.7, col="dark gray") if (graphic) dev.off() # -- Get 10 year projection via Gareth's method for EXTENT # -- this pegs projection against 2010 data point # 10 year projections of 10 year trend #thicke10_2020=lmthicke10[[1]][2]*c(2011:2020)+lmthicke10[[1]][1] thicke10_2020=lmthicke10[[1]][2]*c(1:10)+thicke[32] #thicka10_2020=lmthicka10[[1]][2]*c(2011:2020)+lmthicka10[[1]][1] #thicka10_2020=lmthicka10[[1]][2]*c(1:10)+thicka[32] #vol10_2020=lmvol10[[1]][2]*c(2011:2020)+lmvol10[[1]][1] vol10_2020=lmvol10[[1]][2]*c(1:10)+vol[32] sie10_2020 = vol10_2020/thicke10_2020 #sia10_2020 = vol10_2020/thicka10_2020 # - plot if (graphic) png("img/threnody-gareth-10.png",w=600,h=400) plot(c(vol[1:32],vol10_2020)~c(1979:2020), type="o", col="blue", pch=19, ylim=c(0,13), xlim=c(2001,2020), xaxs="i",yaxs="i",xaxt="n", xlab="", ylab="") axis(1, at=c(1979:2020)) abline(v="2010",col="dark green",lwd=2) lines(c(si[1:32,5],sie10_2020)~c(1979:2020), type="o", col="green", pch=19) lines(c(thicke[1:32],thicke10_2020)~c(1979:2020), type="o", col="red", pch=19) points(vol[33:34]~c(2011:2012),col="blue",pch=13,cex=2) points(si[33:34,5]~c(2011:2012),col="green",pch=13,cex=2) title(main="Reconstructed Threnody for Sea Ice Extent\n10 year trend projected from 2010") text(2002,12, "Volume (1000 km^3)",adj=c(0,0)) text(2002,7, "Extent (million km^2)",adj=c(0,0)) text(2002,2.5, "Thickness (m)",adj=c(0,0)) mtext(paste("http://rhinohide.wordpress.com",format(Sys.time(), "%Y%m%d")), 4, side=1, adj=1, cex=0.7, col="dark gray") if (graphic) dev.off() # -- Get 21 year projection via STANDARD trend for AREA # -- this pegs projection against the previous trend #thicke21_2020=lmthicke1990[[1]][2]*c(2011:2020)+lmthicke1990[[1]][1] #thicke21_2020=lmthicke1990[[1]][2]*c(1:10)+thicke[32] thicka21_2020=lmthicka1990[[1]][2]*c(2011:2020)+lmthicka1990[[1]][1] #thicka21_2020=lmthicka1990[[1]][2]*c(1:10)+thicka[32] vol21_2020=lmvol1990[[1]][2]*c(2011:2020)+lmvol1990[[1]][1] #vol21_2020=lmvol1990[[1]][2]*c(1:10)+vol[32] #sie21_2020 = vol21_2020/thicke21_2020 sia21_2020 = vol21_2020/thicka21_2020 #- plot if (graphic) png("img/threnody-standard-21.png",w=600,h=400) plot(c(vol[1:32],vol21_2020)~c(1979:2020), type="o", col="blue", pch=19, ylim=c(0,13), xlim=c(2001,2020), xaxs="i",yaxs="i",xaxt="n", xlab="", ylab="") axis(1, at=c(1979:2020)) abline(v="2010",col="dark green",lwd=2) lines(c(si[1:32,6],sia21_2020)~c(1979:2020), type="o", col="green", pch=19) lines(c(thicka[1:32],thicka21_2020)~c(1979:2020), type="o", col="red", pch=19) points(vol[33:34]~c(2011:2012),col="blue",pch=13,cex=2) points(si[33:34,6]~c(2011:2012),col="green",pch=13,cex=2) title(main="Reconstructed Threnody for Sea Ice Area\nprojected 21 year trend") text(2002,12, "Volume (1000 km^3)",adj=c(0,0)) text(2002,5, "Area (million km^2)",adj=c(0,0)) text(2002,3, "Thickness (m)",adj=c(0,0)) mtext(paste("http://rhinohide.wordpress.com",format(Sys.time(), "%Y%m%d")), 4, side=1, adj=1, cex=0.7, col="dark gray") if (graphic) dev.off() # -- Get 10 year projection via standard trend # -- this pegs projection against the previous trend # 10 year projections of 10 year trend #thicke10_2020=lmthicke10[[1]][2]*c(2011:2020)+lmthicke10[[1]][1] #thicke10_2020=lmthicke10[[1]][2]*c(1:10)+thicke[32] thicka10_2020=lmthicka10[[1]][2]*c(2011:2020)+lmthicka10[[1]][1] #thicka10_2020=lmthicka10[[1]][2]*c(1:10)+thicka[32] vol10_2020=lmvol10[[1]][2]*c(2011:2020)+lmvol10[[1]][1] #vol10_2020=lmvol10[[1]][2]*c(1:10)+vol[32] #sie10_2020 = vol10_2020/thicke10_2020 sia10_2020 = vol10_2020/thicka10_2020 if (graphic) png("img/threnody-standard-10.png",w=600,h=400) plot(c(vol[1:32],vol10_2020)~c(1979:2020), type="o", col="blue", pch=19, ylim=c(0,13), xlim=c(2001,2020), xaxs="i",yaxs="i",xaxt="n", xlab="", ylab="") axis(1, at=c(1979:2020)) abline(v="2010",col="dark green",lwd=2) lines(c(si[1:32,6],sia10_2020)~c(1979:2020), type="o", col="green", pch=19) lines(c(thicka[1:32],thicka10_2020)~c(1979:2020), type="o", col="red", pch=19) points(vol[33:34]~c(2011:2012),col="blue",pch=13,cex=2) points(si[33:34,6]~c(2011:2012),col="green",pch=13,cex=2) title(main="Reconstructed Threnody for Sea Ice Area\nprojected 10 year trend") text(2002,12, "Volume (1000 km^3)",adj=c(0,0)) text(2002,5, "Area (million km^2)",adj=c(0,0)) text(2002,3, "Thickness (m)",adj=c(0,0)) mtext(paste("http://rhinohide.wordpress.com",format(Sys.time(), "%Y%m%d")), 4, side=1, adj=1, cex=0.7, col="dark gray") if (graphic) dev.off() # -- plot the whole time series if (graphic) png("img/threnody-extended-10.png",w=600,h=400) plot(c(vol[1:32],vol10_2020)~c(1979:2020), type="o", col="blue", pch=19, ylim=c(0,20), xlim=c(1979,2020), xaxs="i",yaxs="i",xaxt="n", xlab="", ylab="") axis(1, at=c(1979:2020)) abline(v="2010",col="dark green",lwd=2) lines(c(si[1:32,6],sia10_2020)~c(1979:2020), type="o", col="green", pch=19) lines(c(thicka[1:32],thicka10_2020)~c(1979:2020), type="o", col="red", pch=19) title(main="Reconstructed Threnody for Sea Ice Area\nprojected 10 year trend") text(1980,17.5, "Volume (1000 km^3)",adj=c(0,0)) text(1980,6, "Area (million km^2)",adj=c(0,0)) text(1980,1.75, "Thickness (m)",adj=c(0,0)) mtext(paste("http://rhinohide.wordpress.com",format(Sys.time(), "%Y%m%d")), 4, side=1, adj=1, cex=0.7, col="dark gray") if (graphic) dev.off()