Tee oma lämpötila-analyysisi

fetchGISS = function(file) {
url = "http://data.giss.nasa.gov/gistemp/tabledata/GLB.Ts dSST.txt"
if (file.exists("giss.txt") == FALSE) download.file(url, "giss.txt")
lines = readLines(file)
lines = lines[9:(length(lines)-11)]
time = double()
temp = double()
k = 1
for (i in 1:length(lines)) {
if (substr(lines[i], 1, 4) != "Year" && nchar(lines[i]) > 65) {
year = as.integer(substr(lines[i], 1, 4))
for (j in 1:12) {
token = substr(lines[i], j*5 1, (j 1)*5)
if (substr(token, 1, 5) == "*****") break
time[k] = year (j-1)/12.0
temp[k] = as.double(token)/100
k = k 1
}
}
}
return(data.frame(cbind(time, temp)))
}
fetchNCDC = function(file) {
url = "ftp://ftp.ncdc.noaa.gov/pub/data/anomalies/monthly.land_ocean.90S.90N.df_1901-2000mean.dat"
if (file.exists("ncdc.txt") == FALSE) download.file(url, "ncdc.txt")
lines = readLines(file)
time = double()
temp = double()
for (i in 1:length(lines)) {
if (substr(lines[i], 9, 18) == "-999.0000") break
year = as.integer(substr(lines[i], 1, 4))
month = as.integer(substr(lines[i], 5, 7))
time[i] = year (month-1)/12.0
temp[i] = as.double(substr(lines[i], 8, 18))
}
return(data.frame(cbind(time, temp)))
}
fetchHadC = function(file) {
url = "http://www.cru.uea.ac.uk/cru/data/temperature/hadcrut3gl.txt"
if (file.exists("hadc.txt") == FALSE) download.file(url, "hadc.txt")
lines = readLines(file)
data = lines[seq(1, length(lines), by=2)]
freq = lines[seq(2, length(lines), by=2)]
time = double()
temp = double()
k = 1
for (i in 1:(length(lines)/2)) {
year = as.integer(substr(data[i], 1, 5))
for (j in 0:11) {
count = as.integer(substr(freq[i], j*7 6, (j 1)*7 5))
if (count == 0) break
time[k] = year j/12.0
temp[k] = as.double(substr(data[i], j*7 6, (j 1)*7 5))
k = k 1
}
}
return(data.frame(cbind(time, temp)))
}
fetchRSS = function(file) {
url = "ftp://ftp.ssmi.com/msu/monthly_time_series/rss_monthly_msu_amsu_channel_tlt_anomalies_land_and_ocean_v03_2.txt"
if (file.exists("rss.txt") == FALSE) download.file(url, "rss.txt")
lines = readLines(file)
lines = lines[4:length(lines)]
time = double()
temp = double()
for (i in 1:length(lines)) {
year = as.integer(substr(lines[i], 1, 5))
month = as.integer(substr(lines[i], 6, 8))
time[i] = year (month-1)/12.0
temp[i] = as.double(substr(lines[i], 9, 17))
}
return(data.frame(cbind(time, temp)))
}
fetchUAH = function(file) {
url = "http://vortex.nsstc.uah.edu/public/msu/t2lt/uahncdc.lt"
if (file.exists("uah.txt") == FALSE) download.file(url, "uah.txt")
lines = readLines(file)
lines = lines[2:(length(lines)-3)]
time = double()
temp = double()
for (i in 1:length(lines)) {
year = as.integer(substr(lines[i], 1, 5))
month = as.integer(substr(lines[i], 6, 8))
time[i] = year (month-1)/12.0
temp[i] = as.double(substr(lines[i], 9, 14))
}
return(data.frame(cbind(time, temp)))
}
truncateData = function(data, max.time) {
time = data$time
temp = data$temp
index = time

##############################################################################
# Kertoo kuinka monen kuukauden trendeja lasketaan.

num.months = 120

##############################################################################
# Lataa tarvittavat kirjastot ja funktiot.

library(nlme)

source("funktiot.R")

##############################################################################
# Hakee aikasarjat netista tai levylta.

giss = fetchGISS("giss.txt")
ncdc = fetchNCDC("ncdc.txt")
hadc = fetchHadC("hadc.txt")
rss = fetchRSS("rss.txt")
uah = fetchUAH("uah.txt")

##############################################################################
# Varmistaa aikasarjojen loppuvan viimeiseen yhteiseen ajanhetkeen.

last.giss = giss$time[nrow(giss)]
last.ncdc = ncdc$time[nrow(ncdc)]
last.hadc = hadc$time[nrow(hadc)]
last.rss = rss$time[nrow(rss)]
last.uah = uah$time[nrow(uah)]

max.time = min(last.giss, last.ncdc, last.hadc, last.rss, last.uah)

giss = truncateData(giss, max.time)
ncdc = truncateData(ncdc, max.time)
hadc = truncateData(hadc, max.time)
rss = truncateData(rss, max.time)
uah = truncateData(uah, max.time)

##############################################################################
# Lyhentaa aikasarjat.

giss = tail(giss, num.months)
ncdc = tail(ncdc, num.months)
hadc = tail(hadc, num.months)
rss = tail(rss, num.months)
uah = tail(uah, num.months)

##############################################################################
# Keskittaa aikasarjat.

giss$temp = giss$temp-mean(giss$temp)
ncdc$temp = ncdc$temp-mean(ncdc$temp)
hadc$temp = hadc$temp-mean(hadc$temp)
rss$temp = rss$temp-mean(rss$temp)
uah$temp = uah$temp-mean(uah$temp)

##############################################################################
# Laskee Lowess-trendit.

smooth1 = lowess(giss, f=24/num.months)
smooth2 = lowess(ncdc, f=24/num.months)
smooth3 = lowess(hadc, f=24/num.months)
smooth4 = lowess(rss, f=24/num.months)
smooth5 = lowess(uah, f=24/num.months)

##############################################################################
# Laskee lineaariset trendit huomioiden autokorrelaation.

model.giss = gls(temp~time, data=giss, corr=corAR1(), method="ML")
model.ncdc = gls(temp~time, data=ncdc, corr=corAR1(), method="ML")
model.hadc = gls(temp~time, data=hadc, corr=corAR1(), method="ML")
model.rss = gls(temp~time, data=rss, corr=corAR1(), method="ML")
model.uah = gls(temp~time, data=uah, corr=corAR1(), method="ML")

##############################################################################
# Laskee 95 % luottamusvalit.

conf.giss = confint(model.giss, level=0.95)
conf.ncdc = confint(model.ncdc, level=0.95)
conf.hadc = confint(model.hadc, level=0.95)
conf.rss = confint(model.rss, level=0.95)
conf.uah = confint(model.uah, level=0.95)

##############################################################################
# Tulostaa lineaariset trendit asteina per vuosikymmen. Suluissa 95 %
# luottamusvalit.

coeff.giss = model.giss$coefficients[2]*10
coeff.ncdc = model.ncdc$coefficients[2]*10
coeff.hadc = model.hadc$coefficients[2]*10
coeff.rss = model.rss$coefficients[2]*10
coeff.uah = model.uah$coefficients[2]*10

lo.giss = conf.giss[2, 1]*10
hi.giss = conf.giss[2, 2]*10

lo.ncdc = conf.ncdc[2, 1]*10
hi.ncdc = conf.ncdc[2, 2]*10

lo.hadc = conf.hadc[2, 1]*10
hi.hadc = conf.hadc[2, 2]*10

lo.rss = conf.rss[2, 1]*10
hi.rss = conf.rss[2, 2]*10

lo.uah = conf.uah[2, 1]*10
hi.uah = conf.uah[2, 2]*10

old.opts = options(digits=2)

cat("Trendit asteina per vuosikymmen (suluissa 95 % luottamusvalit)\n")

cat("GISS = ", coeff.giss, " (", lo.giss, "...", hi.giss, ")\n", sep="")
cat("NCDC = ", coeff.ncdc, " (", lo.ncdc, "...", hi.ncdc, ")\n", sep="")
cat("HADC = ", coeff.hadc, " (", lo.hadc, "...", hi.hadc, ")\n", sep="")
cat("RSS = ", coeff.rss, " (", lo.rss, "...", hi.rss, ")\n", sep="")
cat("UAH = ", coeff.uah, " (", lo.uah, "...", hi.uah, ")\n", sep="")

options(old.opts)

source("laskut.R")

vert.offset = 0.35

month.names = c("tammi", "helmi", "maalis", "huhti", "touko", "kesa",
"heina", "elo", "syys", "loka", "marras", "joulu")

col.giss = "red"
col.ncdc = "#009900"
col.hadc = "blue"
col.rss = "magenta"
col.uah = "#009999"

min.time = giss$time[1]
max.time = giss$time[nrow(giss)]
delta.time = (max.time-min.time)/3
month.index = round((max.time-floor(max.time))*12) 1
num.years = num.months/12

if (num.months == 12) {
plot.title = paste("Viimeisen vuoden")
} else if (num.years == floor(num.years)) {
plot.title = paste("Viimeisten", num.years, "vuoden")
} else {
plot.title = paste("Viimeisten", num.months, "kuukauden")
}
plot.title = paste(plot.title, "havainnot ja trendit")
plot.title = paste(plot.title, "\n", month.names[month.index], "kuuhun", sep="")
plot.title = paste(plot.title, floor(max.time), "asti")

min.temp = min(giss$temp)
max.temp = max(uah$temp 4*vert.offset)

png("kuva1.png", width=600, height=600, pointsize=14)

plot(giss$time, giss$temp, col=col.giss, type="l",
xlim=c(min.time, max.time), ylim=c(min.temp, max.temp),
main=plot.title, xlab="Vuosi", ylab="Anomalia (astetta)")
lines(giss$time, smooth1$y, col=col.giss, lwd=3)
lines(giss$time, model.giss$fitted, col=col.giss, lwd=3, lty="63")

lines(ncdc$time, ncdc$temp vert.offset, col=col.ncdc)
lines(ncdc$time, smooth2$y vert.offset, col=col.ncdc, lwd=3)
lines(ncdc$time, model.ncdc$fitted vert.offset, col=col.ncdc, lwd=3, lty="63")

lines(hadc$time, hadc$temp 2*vert.offset, col=col.hadc)
lines(hadc$time, smooth3$y 2*vert.offset, col=col.hadc, lwd=3)
lines(hadc$time, model.hadc$fitted 2*vert.offset, col=col.hadc, lwd=3, lty="63")

lines(rss$time, rss$temp 3*vert.offset, col=col.rss)
lines(rss$time, smooth4$y 3*vert.offset, col=col.rss, lwd=3)
lines(rss$time, model.rss$fitted 3*vert.offset, col=col.rss, lwd=3, lty="63")

lines(uah$time, uah$temp 4*vert.offset, col=col.uah)
lines(uah$time, smooth5$y 4*vert.offset, col=col.uah, lwd=3)
lines(uah$time, model.uah$fitted 4*vert.offset, col=col.uah, lwd=3, lty="63")

lims = par("usr")

legend(min.time, lims[4], lty=1, lwd=1, legend="Havainnot", bty="n", yjust=0.75)
legend(min.time 0.85*delta.time, lims[4], lty=63, lwd=3,
legend="Lineaariset trendit", bty="n", yjust=0.75)
legend(min.time 2*delta.time, lims[4], lty=1, lwd=3,
legend="Lowess-trendit", bty="n", yjust=0.75)

left.time = (lims[1] min.time)/2
rite.time = (lims[2] max.time)/2

text(left.time, model.giss$fitted[1], "GISTEMP", col=col.giss, adj=c(0.5, 0.5), srt=90)
text(rite.time, model.ncdc$fitted[nrow(ncdc)] vert.offset, "NCDC", col=col.ncdc, adj=c(0.5, 0.5), srt=90)
text(left.time, model.hadc$fitted[1] 2*vert.offset, "HadCRUT3", col=col.hadc, adj=c(0.5, 0.5), srt=90)
text(rite.time, model.rss$fitted[nrow(rss)] 3*vert.offset, "RSS", col=col.rss, adj=c(0.5, 0.5), srt=90)
text(left.time, model.uah$fitted[1] 4*vert.offset, "UAH", col=col.uah, adj=c(0.5, 0.5), srt=90)

dev.off()

source("laskut.R")

library(gplots)

png("kuva2.png", width=600, height=600, pointsize=14)

month.names = c("tammi", "helmi", "maalis", "huhti", "touko", "kesa",
"heina", "elo", "syys", "loka", "marras", "joulu")

min.time = giss$time[1]
max.time = giss$time[nrow(giss)]
delta.time = (max.time-min.time)/3
month.index = round((max.time-floor(max.time))*12)+1
num.years = num.months/12

col.giss = "red"
col.ncdc = "#009900"
col.hadc = "blue"
col.rss = "magenta"
col.uah = "#009999"

if (num.months == 12) {
plot.title = paste("Viimeisen vuoden")
} else if (num.years == floor(num.years)) {
plot.title = paste("Viimeisten", num.years, "vuoden")
} else {
plot.title = paste("Viimeisten", num.months, "kuukauden")
}
plot.title = paste(plot.title, "trendit 95 % luottamusvaleineen")
plot.title = paste(plot.title, "\n", month.names[month.index], "kuuhun", sep="")
plot.title = paste(plot.title, floor(max.time), "asti")

max.trend = max(hi.giss, hi.ncdc, hi.hadc, hi.rss, hi.uah)
min.trend = min(lo.giss, lo.ncdc, lo.hadc, lo.rss, lo.uah)

plotCI(1, coeff.giss, li=lo.giss, ui=hi.giss,
gap=0.2, col=col.giss, lwd=3,xaxt="n",
xlim=c(1, 5), ylim=c(min.trend, max.trend),
main=plot.title,
xlab="", ylab="Trendi (astetta/vuosikymmen)")
plotCI(2, coeff.ncdc, li=lo.ncdc, ui=hi.ncdc,
gap=0.2, add=TRUE, col=col.ncdc, lwd=3)
plotCI(3, coeff.hadc, li=lo.hadc, ui=hi.hadc,
gap=0.2, add=TRUE, col=col.hadc, lwd=3)
plotCI(4, coeff.rss, li=lo.rss, ui=hi.rss,
gap=0.2, add=T, col=col.rss, lwd=3)
plotCI(5, coeff.uah, li=lo.uah, ui=hi.uah,
gap=0.2, add=T, col=col.uah, lwd=3)

abline(h=0, lty=1)

axis(1, 1:5, c("GISTEMP","NCDC","HadCRUT3","RSS","UAH"))

dev.off()

Jotta voisit ajaa skriptit, tarvitset seuraavat:
(1) R-ohjelmiston
(2) NLME-kirjaston R:lle
(3) gplots-kirjaston R:lle

R löytyy Linuxille, Windowsille ja Macille esim. täältä:

http://mirrors.dotsrc.org/cran/

Kun olet asentanut R:n, asenna NLME- ja gplots-kirjastot. Käynnistä R ja kirjoita komentoriville:

install.packages("nlme", dependencies=TRUE)
install.packages("gplots", dependencies=TRUE)

Ja sitten vain ajelemaan skriptejä.

Mainittakoon vielä, että funktiot.R-skriptissä haetaan netistä (WWW ja FTP) lämpötiladatat ja ne tallennetaan kovalevylle. Seuraavilla ajokerroilla käytetään kovalevylle aiemmin haettuja datoja. Näin siksi, että ei aina turhaan haeta samoja datoja netistä. Jos haluat esim. parin kuukauden päästä kokeilla uusimmilla datoilla, tuhoa giss.txt, ncdc.txt, hadc.txt, rss.txt ja uah.txt nimiset tiedostot kovalevyltäsi, jolloin seuraavalla ajokerralla datat haetaan uudelleen netistä.

Meinaatko, että tuosta valmiiksi muokatusta (adjust) lähdeaineistosta saa mitään järkevää irti?
Raakadata olisi hyvä lähde, mutta nuo on jo valmiiksi väännelty. Jostain syystä muokkaamatonta dataa on melkein mahdotonta saada mistään.
Melko turhaa touhua väännellä tuollaisella mitään käppyröitä.