R-ohjelmointi.org

Tilastotieteellistä ohjelmointia R-kielellä

Viimaindeksi ja Ilmatieteenlaitoksen avoin data

Viimaindeksi pyrkii kuvaamaan sitä, miten lämpötila ja tuulennopeus vaikuttavat ihmisen aistimaan kylmyyteen. Selkeän paleltumisvaaran rajana voidaan pitää indeksin, siis aistitun kylmyyden arvoa -35 Celsius-astetta. Ilmatieteen laitoksen sivulta löytyy interaktiivinen visualisointi viimaindeksistä, mutta omaan järkeeni istuu paremmin kaavio, jossa lämpötila on pystyakselilla, ja tuulen nopeus vaaka-akselilla. Tällainen kaavio voidaan laatia R:ssä helposti:

# Viimaindeksi laskettuna arvoille väliltä 0...-52 Celsius-astetta ja 1...31 m/s.
t<- seq(-52, 0)
v<- seq(1, 31)
# Muodostaa kaikki yhdistelmät lämpötilasta ja tuulennopeudesta
m<-expand.grid(t,v)
colnames(m)<-c("t", "v")
# Viimaindeksin kaava Wikipediasta
wcn<-13.12+0.6125*m$t-13.956*m$v^0.16+0.4867*m$t*m$v^0.16
pm<-matrix(ncol=length(v), nrow=length(t), data=wcn)
colnames(pm)<-paste("v", v, sep="")
rownames(pm)<-paste("t", t, sep="")
# Kierretään matriisi oikeaan muotoon
pm2<-pm[53:1,]
tpm2<-t(pm2)
tpm3<-tpm2[,53:1]
 
# Kuva
# Sinisen sävyt
require(RColorBrewer)
cols2<-brewer.pal(9, "Blues")
par(mar=c(4,4,4,5))
# Piirretään käyrät
contour(x=1:nrow(tpm3), y=-(ncol(tpm2):1), z=tpm3, col="grey75", levels=seq(0,-86,by=-5), xlab="Tuulen nopeus (m/s)", ylab="Lämpötila (Celsius-astetta)", main="Viimaindeksi", asp=30/52, bty="c")
# Väritetään käyrät paleltumisvaaran mukaan
# contourLines laskee vain käyrien koordinaatit, muttei piirrä mitään, mutta polygon tekee värityksen
a<-contourLines(x=1:nrow(tpm3), y=-(ncol(tpm2):1), z=tpm3, levels=c(max(tpm3),-25))
polygon(x=c(a[[1]]$x, 31, rev(a[[2]]$x)), y=c(a[[1]]$y, -1, rev(a[[2]]$y)), col=cols2[1], border=cols2[1])
a<-contourLines(x=1:nrow(tpm3), y=-(ncol(tpm2):1), z=tpm3, levels=c(-25,-35))
polygon(x=c(a[[1]]$x, rev(a[[2]]$x)), y=c(a[[1]]$y, rev(a[[2]]$y)), col=cols2[2], border=cols2[2])
a<-contourLines(x=1:nrow(tpm3), y=-(ncol(tpm2):1), z=tpm3, levels=c(-35,-60))
polygon(x=c(a[[1]]$x, rev(a[[2]]$x)), y=c(a[[1]]$y, rev(a[[2]]$y)), col=cols2[4], border=cols2[4])
a<-contourLines(x=1:nrow(tpm3), y=-(ncol(tpm2):1), z=tpm3, levels=c(-60,-86.7))
polygon(x=c(a[[1]]$x, rev(a[[2]]$x)), y=c(a[[1]]$y, rev(a[[2]]$y)), col=cols2[6], border=cols2[6])
polygon(x=c(1, 1, 1.44, 1), y=c(-32.08079, -53, -53, -32.08079), col=cols2[4], border=cols2[4])
# Lisätään vielä uudelleen käyrät
contour(x=1:nrow(tpm3), y=-(ncol(tpm2):1), z=tpm3, col="grey50", levels=seq(0,-86,by=-5), xlab="Tuulen nopeus (m/s)", ylab="Lämpötila (Celsius-astetta)", main="Viimaindeksi", asp=30/52, bty="c", add=T)
contour(x=1:nrow(tpm3), y=-(ncol(tpm2):1), z=tpm3, col="black", levels=c(-25,-35,-60), xlab="Tuulen nopeus (m/s)", ylab="Lämpötila (Celsius-astetta)", main="Viimaindeksi", asp=30/52, add=T)
mtext(side=4, at=-5, text="Kylmää", las=1)
mtext(side=4, at=-13, text="Erittäin\nkylmää", las=1)
mtext(side=4, at=-25, text="Paleltumis-\nvaara", las=1)
mtext(side=4, at=-45, text="Suuri\npaleltumis-\nvaara", las=1)

Tulosmatriisi piti ensin pyörittää sopivaan muotoon, jotta komento contour() osaa piirtää siitä oikeanlaisen kuvan. Tuloksena syntyvä kuva on esitetty alla:

Ilmatieteen laitoksen avoimen datan API mahdollistaa eri paikkakuntien lämpötila- ja tuulennopeus-havaintojen hakemisen ilmaiseksi. Palvelun käyttö vaatii rekisteröitymisen avaimen saamiseksi. Haetaan seuraavaksi Helsingin seudun säätiedot viimeisen viikon ajalta, ja lisätään ne käyräksi kuvioon, jotta nähdään, miten pakkasen purevuus on viime viikolla muuttunut. Manuaalista löytyy esimerkkejä, miten hakuja voidaan tehdä. Otetaan havainnot klo 06:00 aamulla, koska silloin aletaan heräillä töihin lähtöä varten:

# Haetaan ensin tuuli- ja lämpötilatiedot Ilmatieteenlaitoksen avointa API:a käyttäen
# Helsinki
require(RCurl)
dates<-seq.Date(from=as.Date("2013-11-01"), to=as.Date("2014-01-19"), by="day")
m<-data.frame(matrix(nrow=length(dates)-1, ncol=3, data=NA))
rownames(m)<-dates[1:(length(dates)-1)]
colnames(m)<-c("time", "wind", "temp")
for(i in 1:(length(dates)-1)) { 
   dlurl<-paste("http://data.fmi.fi/fmi-apikey/lisää_oma_avaimesi_tähän/wfs?request=getFeature&storedquery_id=fmi::observations::weather::timevaluepair&place=kumpula,Helsinki&parameters=windspeedms,temperature&starttime=", dates[i], "T06:00:00Z&endtime=", dates[i+1], "T05:00:00Z&timestep=60&", sep="")
   test<-getURL(dlurl)
   write(test, "test.txt")
   t1<-readLines("test.txt")
   da<-data.frame( time=t1[grep("<wml2:time>", t1)],
                  value=t1[grep("<wml2:value>", t1)])
   m[i,]<-c(c(substr(as.character(da$time), 50, 59)[1], gsub("</wml2:value> ", "", gsub("\t\t\t\t      <wml2:value>", "", as.character(da$value)))[c(1,25)]))
}
helsinki<-m

Tuloksena saadut tiedot voidaan lisätä kuvaan seuraavalla koodilla:

lines(x=as.numeric(helsinki[73:79,]$wind), y=as.numeric(helsinki[73:79,]$temp), type="b", pch=20)
text(x=as.numeric(helsinki[73:79,]$wind), y=as.numeric(helsinki[73:79,]$temp), labels=paste(substr(helsinki[73:79,]$time, 9, 10), ".", sep=""), cex=0.75, pos=4)

Tuloksena on seuraava kuva:

Tuloksena näyttää olevan, että viime viikon on ollut vain kylmää, vaikka onkin ehkä tuntunut siltä, että ollaan paleltumisen partaalla.


Vastaa

Sähköpostiosoitettasi ei julkaista. Pakolliset kentät on merkitty *