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

313

Ää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. No nytkö tuli lähtö Orpolle?

      Pieniä oli Marinin aamupalasilakat joulukaloiksi vrt. Orpon 35 miljoonan euron kähmintä johonkin Vapaavuoren urheiluhall
      Maailman menoa
      242
      2111

    2. Mikä teidän jutussa on ongelmana?

      Missä meni pieleen?
      Ikävä
      188
      1627

    3. Kauhavan häiriköijistä

      Juttua Iltalehdessä. Pakko sanoa että noi nuoret on kyllä ihan pimeitä. Putkin peltoja jupksevat kiusaamaan kun ei tietä
      Kauhava
      46
      1282

    4. Haluan sinut, kuuletko minua.

      Haluan sinut. Toivon, että voisimme olla yhdessä. Mietin pystynkö täyttämään toiveesi, olemaan arvoisesi. Voisitko saad
      Ikävä
      50
      992

    5. Hän on tosi

      hyvännäköinen. Ei edes ryppyi oo. :D
      Ikävä
      41
      816

    6. Auto ajoi päälle?

      Ja pakeni luin iltapäivälehdestä. ! Ken on kuski joka tuollee teki
      Kuusankoski
      14
      790

    7. Miksi Lapset kiusaa yöllä

      Miksi Lapset kiusaa yöllä ihmisiä? Miksi vanhemmat antaa tämän tapahtua? Eikö ne huomaa ettei lapset ole kotona vai eivä
      Kauhava
      30
      771

    8. Sama ransetti taas!

      Keikkui tällä kertaa Honkavaaran tien varressa muutaman sadan metrin päässä Louhenkoskelta.. Otin rekisterin ylös ja ver
      Hyrynsalmi
      21
      750

    9. Viimeinen lankafest

      Käykää viimeisessä lanka festissä. Ensivuonna sitä ei enää ole. Rahat on loppu. Harmi .
      Puolanka
      24
      735

    10. Tehdäänkö tänään toiveista totta?

      Poikkea tänä illasta siinä lähellä ja annetaan silmien puhua ja sen jälkeen puhu sinä lopulta mitä ajattelet..
      Ikävä
      46
      637
    Aihe
    TietosuojaselosteKäyttöehdotEvästeasetuksetSäännöt