-
Notifications
You must be signed in to change notification settings - Fork 0
/
rclimate-forecast.r
79 lines (67 loc) · 4.21 KB
/
rclimate-forecast.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
######### RClimate Script: Arctic Sea Ice Extent Trend & Forecast for Sept, 2011 ##########
## ## ##
## Original Source data from NSIDC organized into 12 monthly files ##
## Sept data file Link: ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/Sep/N_09_area.txt ##
## I have built and posted a consoldidated monthly csv file available at this link ##
## http://processtrends.com/files/RClimate_NSIDC_sea_ice_extent.csv ##
## R script & forecast developed by D Kelly O'Day ##
## http:chartsgraphs.wordpress.com 2/22/11 ##
## I got the basic idea of quadratic regression from Tamino, Open Mind Post: ##
## http://tamino.wordpress.com/2010/09/14/death-spiral/ ##
## Updated 2/22/11: Reader Edi noticed that I did not define my predict vector a. ##
## Fred added code to generate a png graph (out-1.png)
############################################################################################
link_n <- "http://processtrends.com/files/RClimate_NSIDC_sea_ice_extent.csv"
## Original source data @
link_orig <- "ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/Sep/N_09_area.txt"
sie <- read.table(link_n, skip = 1, sep = ",", header = F,
colClasses = rep("numeric", 5),
comment.char = "#", na.strings = "NA",
col.names = c("yr_frac", "yr", "mo", "ext", "area") )
sie <- sie[,c(1,2,3,4)]
sie <- subset(sie, yr_frac >= 1980) # get 1st full year of data
mn <- 9
sdf <- subset(sie, sie$mo==mn) ## get Sept data
### Quadratic regression: y = a + bx + cx^2
## Center x values # See here for explaination on centering for polynomial regression:
## http://rtutorialseries.blogspot.com/2010/02/r-tutorial-series-basic-polynomial.html
sdf$yr_c <- sdf$yr - mean(sdf$yr,na.rm=T)
sdf$yr_c2 <- sdf$yr_c^2
q_m <- lm(ext ~ yr_c + yr_c2, data=sdf)
a <- predict(q_m)
# Predict 2011 value
future <- data.frame(2011,2011-1994.4, (2011-1994.5)^2)
names(future) <- c("yr", "yr_c", "yr_c2")
f <- predict(q_m, future,interval="confidence")
f_range <- range(f[,2]:f[,3])
###########################################################################
# send to png - original size 1200x1000 is too big
png(file = "out-1.png", width = 800, height = 700)
# Plot annotation strings
y_lab <- expression(paste("Arctic SIE - millions k",m^2))
pl_title <- "September Arctic Sea Ice Extent\n 1980-2010 Trend & Forecast for 2011 - on www.arctic-graphs.net"
proj_2011 <- paste("Forecast Sept, 2011 SIE @ ",signif(f[,1],3), " mil km2 ", "(CI: "
, signif(f[,2],3), " - ", signif(f[,3],3),")",sep="")
## Basic plot
par(las=1); par(ps=10); par(oma=c(3.5,2,1,1)); par(mar=c(2,5,2,1))
par(bg="white"); par(mfrow=c(1,1))
plot(sdf$yr, sdf$ext, type="o",col="grey",pch=16,xlim=c(1980, 2012), ylim=c(3.7,8),
xlab="", ylab = y_lab, main=pl_title, cex=0.75, las=1, pty="m")
# Add quadratic trend line
points(sdf$yr, a, type="l", pch=16,col="blue",lwd=2)
# Add forecast points & line
points(2011, f[,1], type="p", col ="red", pch=16)
points(2011, f[,2], type="p", col="red", pch=16, cex=0.75)
points(2011, f[,3], type="p", col="red", pch=16, cex=0.75)
lines(x=c(2011,2011), y=c(f[,2], f[,3]), type="l", col="red", lwd=1.5)
lines(2011,f[,2], 2011, f[,3], col="black", type="l")
source_note <- paste("Original Data Source: ", link_orig, sep="")
mtext(source_note, 1,0, outer=T, adj=0.5, cex=0.75)
legend(1979,4.5, c("September SIE - mil km2:1980-2010", "Quadratic Trend Line (1980-2010)", proj_2011), col = c("grey", "blue","red"),
text.col = "black",lty = c(1,1,1),lwd=c(1,2,2),pch=c(16,0,16),pt.cex=c(1,0,1),
merge = F, bg = "white", bty="o", box.col = "white", cex = .7)
# Generate and add bottom footer KOD & System data notes
mtext("D Kelly O'Day - http://chartsgraphs.wordpress.com", 1,1, adj = 0, cex = 0.8, outer=T)
mtext(format(Sys.time(), "%m/%d/ %Y"), 1, 1, adj = 1, cex = 0.8, outer=T)
# file's end
dev.off()