############################################################# ##Part VI - Examining Given Years and Given Sea level rises## ############################################################# #NOTE: requires the running of "Tracing sources of variability.R" #we can use the results of part V to ask two related questions: #first: what happens in a given year. For example, what does #Olympia sea level rise look like in the year 2075? yr = 2075 yr.index=yr-1999 #we can ask this question about the median model through med.scens #the gcm model through gcm.scens #and the other models with a little more work #the .info files give 5% and 95% quantiles, but if we want a density we have to do the following for the year gcm.scens.big = list(s.1.smooth,s.2.smooth,s.3.smooth,s.4.smooth) d1.full = density(s.1.full[yr.index,1,]) d2.full = density(s.2.full[yr.index,1,]) d3.full = density(s.3.full[yr.index,1,]) d4.full = density(s.4.full[yr.index,1,]) densities=list(d1.full,d2.full,d3.full,d4.full) bi1.full = density(bi.1.full[yr.index,]) bi2.full = density(bi.2.full[yr.index,]) bi3.full = density(bi.3.full[yr.index,]) bi4.full = density(bi.4.full[yr.index,]) bi.densities = list(bi1.full,bi2.full,bi3.full,bi4.full) #to make this picture, note for different years the #scales on this image will need to be reworked png("2075Estimates.png",height=700,width=700) par(mfrow=c(2,2),pty='s') for(i in 1:4) { hist(model.scens[[i]][yr.index,],xlim=c(15,75),prob=T,main=titles[i],xaxt='n',xlab="",ylim=c(0,.14)) axis(1,at=2.54*5*(2:5),labels=seq(from=10,to=25,by=5),col.axis='red',col.ticks='red',col='white') axis(1,at=seq(from=15,to=75,by=15)) mtext("Anomaly",side=1,line=2.35,cex=.9) mtext("(cm) ",side=1,line=3.45,cex=.9) mtext(" (in)",side=1,line=3.45,col='red',cex=.9) lines(density(gcm.scens.big[[i]][yr.index,]),col='red',lwd=2) lines(densities[[i]],col='blue',lwd=2) lines(bi.densities[[i]],col='green',lwd=2) #polygon(bi.densities75shift[[i]],col=rgb(0,1,0,.5)) abline(v=med.scens[[i]][yr.index],lwd=3) } dev.off() #what about years of various rises? #one option is to ask about just 5%, median, and 95% levels #we can do that using the following code: #say we are interested in .5ft,1ft, and 2 ft lvs.ft = c(.5,1,2) lvs = lvs.ft*12*2.54 n.lvs=length(lvs) rise.tables = list(scen2.6=NA,scen4.5=NA,scen6.0=NA,scen8.5=NA) for(j in 1:4) { scen = bi.scens[[j]] med.scen = med.scens[[j]] min.by.year = scen[1,] med.by.year = scen[2,] max.by.year = scen[3,] year = 2000:21000 year.range = matrix(NA,nrow=3,ncol=n.lvs) for(i in 1:n.lvs) { year.range[i,1] <- year[which(max.by.year >= lvs[i])[1]] year.range[i,2] = year[which(med.by.year >= lvs[i])[1]] year.range[i,3] <- year[which(min.by.year >= lvs[i])[1]] } final.table = cbind(level=lvs.ft,year.range) colnames(final.table)[2:4] = c("P0.05","P0.5","P0.95") rise.tables[[j]]=final.table } #we can also make a richer figure for a specific rise, say 1 ft #get.counts will get the number of paths above a certain level get.counts = function(x,level=lv) { return(sum(x>=level)) } lv.ft = 1 lv = lv.ft*12*2.54 med.years = scenario.years = gcm.years = seattle.years = bi.years = matrix(NA,nrow=101,ncol=4) for(j in 1:4) { med.years[,j] = sapply(med.scens[[j]],get.counts) scenario.years[,j] = apply(model.scens[[j]],1,get.counts) gcm.years[,j] = apply(gcm.scens.big[[j]],1,get.counts) } seattle.years[,1] = apply(s.1.full,c(1,2),get.counts) seattle.years[,2] = apply(s.2.full,c(1,2),get.counts) seattle.years[,3] = apply(s.3.full,c(1,2),get.counts) seattle.years[,4] = apply(s.4.full,c(1,2),get.counts) bi.years[,1] = apply(bi.1.full,1,get.counts,level=(lv)) bi.years[,2] = apply(bi.2.full,1,get.counts,level=(lv)) bi.years[,3] = apply(bi.3.full,1,get.counts,level=(lv)) bi.years[,4] = apply(bi.4.full,1,get.counts,level=(lv)) #Now we can plot the proportion of paths above 1 ft #This is not quite a CDF, but it looks like one png("1ftRise.png",height=700,width=700) par(mfrow=c(2,2),pty='s') for(i in 1:4) { plot(yrs.to.plot,med.years[,i],type='l',col='black',lwd=2,main=titles[i],xlab='Year',ylab='Proportion of simlations above 20 cm') lines(yrs.to.plot,scenario.years[,i]/18,type='l',col='grey',lwd=2) lines(yrs.to.plot,gcm.years[,i]/nsamples,col='red',lwd=2) lines(yrs.to.plot,seattle.years[,i]/nsamples2,col='blue',lwd=2) lines(yrs.to.plot,bi.years[,i]/nsamples2,col='green',lwd=2) } dev.off()