forskning.no
function (uah = TRUE, rss = TRUE, giss = TRUE, ncdc = TRUE, hadcrut3 = TRUE, 
    lwd = 1, ylim = c(-0.4, 1), xlim = c(1995, 2012), base.period = 1981:2010, 
    type = "s") 
{
    anomaly <- function(X, yr, mo, base.period) {
        anom <- X + NA
        iiref <- is.element(yr, base.period)
        for (im in 1:12) {
            iim <- is.element(mo, im)
            anom[iim] <- X[iim] - mean(X[iim & iiref], na.rm = TRUE)
        }
        anom
    }
    annual <- function(X, date) {
        yr <- trunc(date)
        yrs <- as.numeric(rownames(table(yr)))
        anm <- yrs + NA
        for (i in 1:length(anm)) {
            ii <- is.element(yr, yrs[i])
            if (sum(ii) == 12) 
                anm[i] <- mean(X[ii])
        }
        anm
    }
    L <- function(yrs) {
        l = as.numeric(rownames(table(trunc(yrs))))
        l
    }
    strleg <- c("UAH", "RSS", "GISS", "NCDC", "HadCRUT3")
    colleg <- c("blue", "black", "red", "violet", "darkgreen")
    ileg <- c(uah, rss, giss, ncdc, hadcrut3)
    uahmsuurl <- "http://vortex.nsstc.uah.edu/data/msu/t2lt/uahncdc.lt"
    rssmsuurl <- "http://www.remss.com/data/msu/monthly_time_series/RSS_Monthly_MSU_AMSU_Channel_TLT_Anomalies_Land_and_Ocean_v03_3.txt"
    ncdcurl <- "ftp://ftp.ncdc.noaa.gov/pub/data/anomalies/monthly.land_ocean.90S.90N.df_1901-2000mean.dat"
    hadcrut3url <- "http://www.metoffice.gov.uk/hadobs/hadcrut3/diagnostics/global/nh+sh/monthly"
    gissurl <- "http://data.giss.nasa.gov/gistemp/tabledata_v3/GLB.Ts+dSST.txt"
    ref <- paste(min(base.period), max(base.period), sep = " - ")
    dev.new()
    par(bty = "n", las = 2)
    plot(x = xlim, y = ylim, ylim = ylim, xlim = xlim, type = "n", 
        main = "Replication of Solheim et al in forsking.no (March 2012)", 
        sub = "forskning.no (replicationDemos)", xlab = "time", 
        ylab = paste("anomaly wrt", ref, "(deg C)"))
    grid()
    par(las = 0)
    mtext("http://www.forskning.no/artikler/2012/mars/316178", 
        side = 4, col = "grey")
    legend(xlim[1], ylim[2], strleg[ileg], col = colleg[ileg], 
        bty = "n", lty = 1, lwd = lwd, cex = 0.75)
    if (uah) {
        print("UAH")
        test <- readLines(uahmsuurl)
        lwy <- grep("Year", test)
        nrows <- lwy[2] - 2
        uahmsu <- read.table(uahmsuurl, header = TRUE, nrows = nrows)
        uahtemp <- anomaly(uahmsu$Globe, uahmsu$Year, uahmsu$Mo, 
            base.period)
        lines(uahmsu$Year + (uahmsu$Mo - 0.5)/12, uahtemp, lwd = lwd, 
            col = "blue", type = type)
    }
    if (rss) {
        print("RSS")
        rsscn <- c("year", "mon", "S70.0toN82.5", "S20.0toN20.0", 
            "N20.0toN82.5", "S70.0toS20.0", "N60.0toN82.5", "S70.0toS60", 
            "Cont.USA", "EqtoN82.5", "S70.0toEq")
        rssmsu <- read.table(rssmsuurl, skip = 3, col.names = rsscn)
        rsstemp <- anomaly(rssmsu$S70.0toN82.5, rssmsu$year, 
            rssmsu$mon, base.period)
        lines(rssmsu$year + (rssmsu$mon - 0.5)/12, rsstemp, lwd = lwd, 
            type = type)
    }
    if (giss) {
        print("GISS")
        test <- readLines(gissurl)
        cyrs <- substr(test, 1, 4)
        vdat <- is.element(cyrs, as.character(1880:2050))
        writeLines(con = "giss.test.txt", text = test[vdat])
        widths = c(rep(5, 13), 7, rep(5, 6))
        gisscn <- c("year", "Jan", "Feb", "Mar", "Apr", "May", 
            "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec", 
            "Jan.Dec", "Dec.Nov", "DJF", "MAM", "JJA", "SON", 
            "Year")
        GISS <- read.fwf("giss.test.txt", widths = widths, col.names = gisscn, 
            na.strings = "*****")
        gisstemp <- c(t(GISS[, 2:13])) * 0.01
        gissdate <- sort(rep(GISS[, 1], 12)) + (rep(1:12, length(GISS[, 
            1])) - 0.5)/12
        gistemp <- anomaly(gisstemp, sort(rep(GISS[, 1], 12)), 
            rep(1:12, length(GISS[, 1])), base.period)
        lines(gissdate, gistemp, lwd = lwd, col = "red", type = type)
    }
    if (ncdc) {
        print("NCDC")
        NCDC <- read.table(ncdcurl, col.names = c("year", "mon", 
            "global"))
        NCDC$global[NCDC$global < -99] <- NA
        ncdctemp <- anomaly(NCDC$global, NCDC$year, NCDC$mon, 
            base.period)
        lines(NCDC$year + (NCDC$mon - 0.5)/12, ncdctemp, col = "violet", 
            lwd = lwd, type = type)
    }
    if (hadcrut3) {
        print("HadCRUT3")
        HadCRUT3 <- read.table(hadcrut3url)
        crudate <- as.numeric(substr(HadCRUT3$V1, 1, 4)) + (as.numeric(substr(HadCRUT3$V1, 
            6, 7)) - 0.5)/12
        crutemp <- anomaly(HadCRUT3$V2, as.numeric(substr(HadCRUT3$V1, 
            1, 4)), as.numeric(substr(HadCRUT3$V1, 6, 7)), base.period)
        lines(crudate, crutemp, col = "darkgreen", type = type)
    }
    dev.new()
    par(bty = "n", las = 2)
    plot(x = xlim, y = ylim, ylim = ylim, xlim = xlim, type = "n", 
        main = "Solheim et al in forsking.no: annual means", 
        sub = "forskning.no (replicationDemos)", xlab = "time", 
        ylab = paste("anomaly wrt", ref, "(deg C)"))
    grid()
    par(las = 0)
    mtext("http://www.forskning.no/artikler/2012/mars/316178", 
        side = 4, col = "grey")
    trendleg <- rep(" ", length(strleg))
    if (uah) {
        UAH <- data.frame(x = L(uahmsu$Year), y = annual(uahtemp, 
            uahmsu$Year))
        lines(UAH, col = "blue", type = "b", pch = 19, cex = 1.2)
        uahtrend <- lm(y ~ x, data = UAH)
        uahstat <- summary(uahtrend)
        lines(UAH$x[is.finite(UAH$y)], predict(uahtrend), col = "blue", 
            lty = 2)
        trendleg[1] <- paste(round(uahstat$coefficients[2] * 
            10, 2), " [", round(2 * uahstat$coefficients[4] * 
            10, 2), "]", sep = "")
    }
    if (rss) {
        RSS <- data.frame(x = L(rssmsu$year), y = annual(rsstemp, 
            rssmsu$year))
        lines(RSS, type = "b", pch = 19, cex = 1.2)
        rsstrend <- lm(y ~ x, data = RSS)
        rssstat <- summary(rsstrend)
        lines(RSS$x[is.finite(RSS$y)], predict(rsstrend), col = "black", 
            lty = 2)
        trendleg[2] <- paste(round(rssstat$coefficients[2] * 
            10, 2), " [", round(2 * rssstat$coefficients[4] * 
            10, 2), "]", sep = "")
    }
    if (giss) {
        x <- L(gissdate)
        y <- annual(gistemp, gissdate)
        GIS <- data.frame(x = x[(x >= xlim[1]) & (x <= xlim[2])], 
            y = y[(x >= xlim[1]) & (x <= xlim[2])])
        lines(GIS, col = "red", type = "b", pch = 19, cex = 1.2)
        gisstrend <- lm(y ~ x, data = GIS)
        gissstat <- summary(gisstrend)
        lines(GIS$x[is.finite(GIS$y)], predict(gisstrend), col = "red", 
            lty = 2)
        trendleg[3] <- paste(round(gissstat$coefficients[2] * 
            10, 2), " [", round(2 * gissstat$coefficients[4] * 
            10, 2), "]", sep = "")
    }
    if (ncdc) {
        x <- L(NCDC$year)
        y <- annual(ncdctemp, NCDC$year)
        NCD <- data.frame(x = x[(x >= xlim[1]) & (x <= xlim[2])], 
            y = y[(x >= xlim[1]) & (x <= xlim[2])])
        lines(NCD, col = "violet", type = "b", pch = 19, cex = 1.2)
        ncdctrend <- lm(y ~ x, data = NCD)
        ncdcstat <- summary(ncdctrend)
        lines(NCD$x[is.finite(NCD$y)], predict(ncdctrend), col = "violet", 
            lty = 2)
        trendleg[4] <- paste(round(ncdcstat$coefficients[2] * 
            10, 2), " [", round(2 * ncdcstat$coefficients[4] * 
            10, 2), "]", sep = "")
    }
    if (hadcrut3) {
        x <- L(crudate)
        y <- annual(crutemp, crudate)
        CRU <- data.frame(x = x[(x >= xlim[1]) & (x <= xlim[2])], 
            y = y[(x >= xlim[1]) & (x <= xlim[2])])
        lines(CRU, col = "darkgreen", type = "b", pch = 19, cex = 1.2)
        crutrend <- lm(y ~ x, data = CRU)
        crustat <- summary(crutrend)
        lines(CRU$x[is.finite(CRU$y)], predict(crutrend), col = "darkgreen", 
            lty = 2)
        trendleg[5] <- paste(round(crustat$coefficients[2] * 
            10, 2), " [", round(2 * crustat$coefficients[4] * 
            10, 2), "]", sep = "")
    }
    legend(xlim[1], ylim[2], paste(strleg[ileg], trendleg[ileg]), 
        col = colleg[ileg], bty = "n", lty = 1, lwd = lwd, cex = 0.75)
}