f = "http://www.psmsl.org/data/obtaining/rlr.annual.data/10.rlrdata"
pmsl = read.table(f,sep=";") # rlr for SF

pmsl[,2] = pmsl[,2]/25.4 # guess Willis doesn't like metric

w = 33/nrow(pmsl) # frac for 33 year window
l = lowess(pmsl[,1],pmsl[,2],f=w)

y1 = l$y[length(l$y)] # good for last point ...
x1 = l$x[length(l$x)] # ... but point was from 2000

y1 = l$y[146] # good for last point ...
x1 = l$x[146] # ... but point was from 2000

# Sea-Level Rise for the Coasts of California, Oregon, and Washington:
# Past, Present, and Future
# data from Table 5.3
# http://www.nap.edu/catalog.php?record_id=13389
x2 = 2030
# outside range
y2a = (y1+4.3/2.54)
y2b = (y1+30/2.54)
# main projection
yy2 = y1 + 14.4/2.54
yy2a = y1 + 9.4/2.54
yy2b = y1 + 19.4/2.54

png("NAS-SF-SLR-2030.png",w=600,h=600)
	plot(pmsl[,2]~pmsl[,1],pch=20,col="lightblue",
		ylab="Sea Level (inches,arbitrary zero point)",
		xlab="",
		xaxt="n",
		yaxt="n",
		xlim=c(1850,2030),
		ylim=c(266,296))
	axis(1, at=seq(1850,2030,20), las=2)
	axis(2, at=seq(266,296,5), las=2)
	title(main="By 2030, Sea Level Off California to Rise\nBetween 2 Inches and 1 Foot\n(6 inch median projection)")
	lines(pmsl[,2]~pmsl[,1],col="lightblue",lwd=2)
	lines(l,col="blue",lwd=2)
	mtext(paste("http://rhinohide.wordpress.com",format(Sys.time(), "%Y%m%d")), 4, side=1, adj=1, cex=0.7, col="dark gray")

	# projection
	polygon(c(x1,x2,x2),c(y1,y2a,y2b),col="pink",lwd=3,border=F)
	polygon(c(x1,x2,x2),c(y1,yy2a,yy2b),col="red",lwd=3,border=F)
	lines(c(x1,2030),c(y1,yy2),col="darkred",lwd=3)

dev.off()