Tee oma lämpötila-analyysisi

R-skriptejä

Postaan tähän viestiketjuun R-kielisen lähdekoodin, jolla loin esim. täällä esitetyt kuvat:

http://keskustelu.suomi24.fi/node/8443809

Seuraavissa kolmessa postauksessa tulee 4 skriptiä, jotka kaikki tarvitaan: funktiot.R, laskut.R, kuva1.R ja kuva2.R

Laitan viestin otsikoksi tiedoston nimen ja sisällöksi skriptin. Eli kopioi viestin sisältö otsikon nimiseen tiedostoon (esim. viestin funktiot.R sisällön kopioit tiedostoon nimeltä funktiot.R).

Kun olet kopioinut skriptit koneellesi, voit luoda kuvat ajamalla skrpitit kuva1.R ja kuva2.R. Kirjoita R-komentoriville:

source("kuva1.R")
source("kuva2.R")

Sen jälkeen sinulla on kaksi kuvaa: kuva1.png ja kuva2.png

Jos haluat kokeilla eri pituisia trendejä, voit muuttaa skriptin laskut.R alussa olevaa muuttujaa num.months eli kuukausien lukumäärä.

11

311

Äänestä

    Vastaukset

    Anonyymi (Kirjaudu / Rekisteröidy)
    5000
    • aloittaja

      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

    • aloittaja

      ##############################################################################
      # 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)

    • aloittaja

      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()

    • aloittaja

      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()

    • aloittaja

      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ä.

    • rtynrtynbr

      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ä.

      • trshsrsrhsr

        Eli kaikki maailman ilmastontutkimuslaitokset ovat katalassa salaliitossa teitä hörhöjä vastaan, vai?

      • trshsrsrhrs
        trshsrsrhsr kirjoitti:

        Eli kaikki maailman ilmastontutkimuslaitokset ovat katalassa salaliitossa teitä hörhöjä vastaan, vai?

        Siltä se vaikuttaa, koska jopa mitattuja lämpötiloja pitää pantata. Todella outoa.

      • tyjdtyjdtyj
      • tyjdtyjdtjy
      • allinko
        tyjdtyjdtjy kirjoitti:

        Edelleenkään sitä muokkaamatonta dataa ei ole saatavilla.

        joten kyseessä on katala juonittelu. Muuten paljastuisi, että olemme menossa kohti seuraavaa jääkautta ja äkkiä.

    Ketjusta on poistettu 0 sääntöjenvastaista viestiä.

    Takaisin ylös

    Luetuimmat keskustelut

    1. Jussi Halla-aho huolissaan Sofia Virrasta

      Jussihan on vanha vihreä. Onko tässä kyse alkukesän kiimasta, kun aidan toisella puolella oleva vihreä alkaa kiinnostama
      Maailman menoa
      91
      7600

    2. Sofia Virta kadonnut....onko juomassa?

      Virran poissaolo eduskunnasta on herättänyt huomiota. Esimerkiksi Ilta-Sanomat kertoi aiemmin, että Virta on ollut tällä
      Maailman menoa
      123
      6324

    3. Julkista rahaa ei tule antaa senttiäkään yksityisille yrityksille

      Julkinen raha on meidän yhteistä rahaa, ja se raha on tarkoitettu yhteiseen käyttöön, kuten esimerkiksi tuottamaan palve
      Maailman menoa
      162
      4452

    4. Tytti Tuppurainen: Suomen pakolaiskiintiö pitäisi nostaa 10 000 vuodessa

      asia on faktaa, noin Tytti sanoi aiemmin. Kun taas Orpon hallitusohjelman mukaisesti Suomen pakolaiskiintiö on pudotettu
      Maailman menoa
      235
      3061

    5. Halla-aho sivaltaa edustajantyöstään lintsaavaa Sofia Virtaa

      https://www.iltalehti.fi/politiikka/a/937c74d7-f905-4466-b9b4-abd017fe5b63 Kansanedustajan on ilmoitettava poissaolosta
      Maailman menoa
      90
      2638

    6. Ruotsissa uusi monikulttuurisuusongelma: Mummonraiskuut

      Ilmiö räjähti käsiin ja nyt painetaan paniikkinappulaa. Moni vanhustenhoivayhtiö on joutunut jopa lopettamaan, koska keh
      Maailman menoa
      35
      1935

    7. Tunnustusten lauantai

      Mitä haluat sanoa kaivatullesi?
      Ikävä
      129
      1551

    8. Yhteydenotto

      Tiedätkö tai ymmärrätkö syyn, miksi kaivattusi ei ota sinuun yhteyttä? Mikä se syy on?
      Ikävä
      197
      1534

    9. Sofia Virta SUURI POLIITIKKO

      Osallistumalla Erikoisjoukkoihin nostaa Vihreät kauaksi ohi perussuomalaisista, joka on muutenkin hajoamassa omaan mahdo
      Maailman menoa
      41
      823

    10. Minkä arvosanan 4-10 annat Susanna Laineelle Farmi-juontohommista?

      Minkä arvosanan 4-10 annat Suskille? Tätä ei tv:ssä: Susanna Laine paljastaa, mikä yksi asia hermostuttaa Farmi-kuvauk
      Tv-sarjat
      21
      782
    Aihe
    TietosuojaselosteKäyttöehdotEvästeasetuksetSäännöt