Jag har i år tänkt börja använda statistikprogrammet R lite mer än tidigare. Förut har jag vanligtvis i huvudsak använt Stata för analyser, men jag har insett att R har andra möjligheter. För att komma igång har jag bland annat startat den här bloggen, som är skapad i R med hjälp av paketet Blogdown. När man ska lära sig ett nytt statistikprogram är det bäst att ge sig i kast med någon uppgift och då jag länge har tänkt göra en svensk version av Edward Tuftes klassiska visualisering av vädret i New Yor kändes det som en lämplig uppgift. För mer bakgrund, se min vanliga blogg.
Det visade sig dock att Brad Boehmke redan tidigare hade gjort något liknande i R. Hans graf skiljer sig något från NY Times graf, framför allt genom att han inte hade tillgång till högsta och lägsta temperatur för respektive dygn, utan endast dygnsmedeltemperaturen. I hans grafen blir därmed det aktuella årets temperatur över dagar illustrerat som en linje. Dessutom hade han inte tillgång till nederbörd, så den delen av den ursprungliga grafen finns inte med.
å jag har utgått från Boehmkes kod, men jag har ändrat den en del. Bland annat fick jag inte hans numrering av dagar att fungera och jag förstod inte varför han ritade samma graf för varje historiskt år. Dessutom valde Boehmke att låta intervallet avseende normala år representeras av konfidensintervallet för det historiska medelvärdet. Detta då han som sagt inte hade tillgång till varje dygns högsta och lägsta temperatur. Jag tycker inte att det är helt lyckat val då det visar osäkerheten i en skattning av ett historiskt medelvärde snarare än ett intervall som innefattar normala värden.
Man skulle i stället kunna låta intervallet motsvara två standardavvikelser i datamaterialet (dvs. inte dela med roten ur n) i förhållande till medelvärdet. Men jag föredrar att låta intervallet “normala temperaturer” avse de mittersta 50 procenten av tidigare års medeltemperaturer. Med andra ord: om man sorterar alla tidigare års medeltemperaturer för en viss dag från den lägsta till den högsta temperaturen så blir de procent av dagarna som ligger i mitten de dagar som betecknas “normala”. Staplarna i grafen blir därmed en variant av ett lådagram (eller boxplot på engelska), där den mörkare färgen blir själva lådan.
Samtidigt är jag inte någon rutinerad användare av R, så det finns ytterligare utrymme för förbättring. I vilket fall, jag har utgått från hans kod för att skapa en graf för Stockholm 2016. Vi börjar med att läsa in data från SMHI över medeltemperatur i Stockholm över tid. Jag hittad bara data uppdelat i två filer: en för historiska data (från 2916 och framåt) och en för de senaste 4 månaderna. Jag misstänker att det kan finnas något api någonstans som jag inte har sett, så den här delen av koden är lite stökig i nuvarande form. Men vi läser in, lägger ihop och fixar till de två filerna.
library(dplyr)
library(tidyr)
library(tidyverse)
library(lubridate)
library(ggplot2)
# Läs in data från SMHI, två semikolon-separerade filer - en med historiska data, en med mer aktuella
# Ingen snygg lösning, men funkar
file_hist <- "C:/Dropbox/data/väder/smhi-opendata_2_98210_20170118_165135.csv"
sthm_hist <- read_csv2(file_hist,skip=9) #ges felmeddelande för vissa kolumner, men det är ok.
file_recent <- "C:/Dropbox/data/väder/smhi-opendata_2_98210_20170122_172654.csv"
sthm_recent <- read_csv2(file_recent,skip=9) #ges felmeddelande för vissa kolumner, men det är ok.
# Byt namn på kolumner
names(sthm_hist) <- c("date_full", "date2", "date", "temp", "quality", "noise", "noise2")
names(sthm_recent) <- c("date_full", "date2", "date", "temp", "quality", "noise", "noise2")
sthm_recent2 <- sthm_recent %>% filter(date > "2016-09-30") # Ta bort de datum som redan finns i historiska data
sthm <- bind_rows(sthm_hist, sthm_recent2) # Lägg ihop dataseten
# Fixa till och rensa bort en del variabler
weather <- sthm %>%
mutate(year = year(date),
month = month(date),
day = day(date),
temp = as.numeric(temp)) %>%
select(temp, year, month, day, date)Sedan skapar vi fyra nya dataset: ett med genomsnittliga temperaturer historiskt sett, ett med temperaturer 2016, ett rekordhöga temperaturer 2016 och ett med rekordlåga temperaturer 2016. Här avviker min kod ganska mycket från Boehmkes.
# skapa ett dataset med aggregerade data för åren 1961-2015
past <- weather %>%
mutate(leap=(month==2 & day==29))%>%
filter(leap==FALSE) %>% # Ta bort extra dag vid skottår
arrange(year,month, day) %>%
group_by(year) %>%
mutate(day_no = seq_len(n())) %>%
ungroup() %>%
filter(year < 2016) %>% # ta bort data för 2016 och senare
group_by(day_no) %>%
summarise(upper = max(temp), # Högsta temperatur
lower = min(temp), # Lägsta temperatur
q1 = quantile(temp, probs = 0.25), # Kvartil 1
q3 = quantile(temp, probs = 0.75)) # Kvatil 3
# Skapa ett dataset med temperaturer 2016
present <- weather %>%
mutate(leap=(month==2 & day==29)) %>%
filter(leap==FALSE) %>%
arrange(year,month, day) %>%
group_by(year) %>%
mutate(day_no = seq_len(n())) %>%
ungroup() %>%
filter(year == 2016) # Selektera 2016
# Skapa ett dataset med de dagar 2016 som hade lägre temperaturer än tidigare år
present_lows <- present %>%
left_join(past, by = "day_no") %>% # matcha historiska data med 2016 års data
mutate(record = ifelse(temp<lower, "Y", "N")) %>% # identifiera dagar m rekordlåg temp
filter(record == "Y") # selektera dagar m rekordlåg temp
# Skapa ett dataset med de dagar 2016 som hade lägre temperaturer än tidigare år
present_highs <- present %>%
left_join(past, by = "day_no") %>% # matcha historiska data med 2016 års data
mutate(record = ifelse(temp>upper, "Y", "N")) %>% # identifiera dagar m rekordlåg temp
filter(record == "Y") # selektera dagar m rekordlåg tempSedan behöver vi fixa till lite med lablarna på y-axeln och skapa ett litet dataset med slumpmässig information för legenden.
# funktion för att lägga till grad-symbol till lablar på y-axeln
dgr_fmt <- function(x, ...) {
parse(text = paste(x, "*degree", sep = ""))
}
# skapa skalsteg på y-axeln
a <- dgr_fmt(seq(-25,35, by=5))
# skapa ett litet dataset för legenden avseende 2016 år temperatur
legend_data <- data.frame(x=seq(173,182),y=rnorm(10,-10,1))Då är det bara att börja bygga ihop grafen. I ett första steg lägger vi in våra tre sorters värden: högsta/lägsta historiskt, normala värden historiskt (de 50 procent av åren som ligger i mitten av fördelningen, dvs. Q1-Q3) och dygnsmedeltemperaturen 2016.
p <- ggplot(past, aes(day_no, upper)) +
theme(plot.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank()) +
geom_linerange(past, mapping=aes(x=day_no, ymin=lower, ymax=upper), colour = "wheat2")
p <- p +
geom_linerange(past, mapping=aes(x=day_no, ymin=q1, ymax=q3), colour = "wheat4")
p <- p +
geom_line(present, mapping=aes(x=day_no, y=temp, group=1)) +
geom_vline(xintercept = 0, colour = "wheat4", linetype=1, size=1)Sedan fixar vi till det lite som Tufte gillar, med vita horisontella linjer och prickade vertikala linjer, samt lägger till månadernas namn längs x-axeln.
p <- p +
geom_hline(yintercept = -25, colour = "white", linetype=1) +
geom_hline(yintercept = -20, colour = "white", linetype=1) +
geom_hline(yintercept = -15, colour = "white", linetype=1) +
geom_hline(yintercept = -10, colour = "white", linetype=1) +
geom_hline(yintercept = -5, colour = "white", linetype=1) +
geom_hline(yintercept = 0, colour = "white", linetype=1) +
geom_hline(yintercept = 5, colour = "white", linetype=1) +
geom_hline(yintercept = 10, colour = "white", linetype=1) +
geom_hline(yintercept = 15, colour = "white", linetype=1) +
geom_hline(yintercept = 20, colour = "white", linetype=1) +
geom_hline(yintercept = 25, colour = "white", linetype=1) +
geom_hline(yintercept = 30, colour = "white", linetype=1) +
geom_hline(yintercept = 35, colour = "white", linetype=1)
p <- p +
geom_vline(xintercept = 31, colour = "wheat4", linetype=3, size=.5) +
geom_vline(xintercept = 59, colour = "wheat4", linetype=3, size=.5) +
geom_vline(xintercept = 90, colour = "wheat4", linetype=3, size=.5) +
geom_vline(xintercept = 120, colour = "wheat4", linetype=3, size=.5) +
geom_vline(xintercept = 151, colour = "wheat4", linetype=3, size=.5) +
geom_vline(xintercept = 181, colour = "wheat4", linetype=3, size=.5) +
geom_vline(xintercept = 212, colour = "wheat4", linetype=3, size=.5) +
geom_vline(xintercept = 243, colour = "wheat4", linetype=3, size=.5) +
geom_vline(xintercept = 273, colour = "wheat4", linetype=3, size=.5) +
geom_vline(xintercept = 304, colour = "wheat4", linetype=3, size=.5) +
geom_vline(xintercept = 334, colour = "wheat4", linetype=3, size=.5) +
geom_vline(xintercept = 365, colour = "wheat4", linetype=3, size=.5)
p <- p +
coord_cartesian(ylim = c(-25,35)) +
scale_y_continuous(breaks = seq(-25,35, by=5), labels = a) +
scale_x_continuous(expand = c(0, 5),
breaks = c(15,45,75,105,135,165,195,228,258,288,320,350),
labels = c("Januari", "Februari", "Mars", "April",
"Maj", "Juni", "Juli", "Augusti", "September",
"Oktober", "November", "December"))Sedan kan vi markera de dagar som hade rekordhög eller rekordlåg temperatur.
p <- p +
geom_point(data=present_lows, aes(x=day_no, y=temp), colour="blue3") +
geom_point(data=present_highs, aes(x=day_no, y=temp), colour="firebrick3")Och slutligen lägger vi till rubrik, förklarande texter och legend. Den här delen av koden kan utvecklas så att det inte är så mycket manuell anpassning. Det vore att föredra om man vill köra samma kod för någon annan stad. Men det får ligga på listan över saker att göra.
p <- p +
ggtitle("Vädret i Stockholm 2016") +
theme(plot.title=element_text(face="bold",hjust=.012,vjust=.8,colour="#3C3C3C",size=20)) +
annotate("text", x = 5, y = 37.5, label = "Medeltemperatur", size=4, fontface="bold", hjust = 0)
#print(p)
p <- p +
annotate("text", x = 5, y = 35,
label = "Data avser dygnsmedeltemperatur för Stockholm år 2016 i förhållande till",
size=3, colour="gray30", hjust = 0) +
annotate("text", x = 5, y = 33,
label = "hur det har sett ut 1961-2015. Genomsnittstemperaturen för hela 2016 var 8,2°,",
size=3, colour="gray30", hjust = 0) +
annotate("text", x = 5, y = 31,
label = "vilket är 1° över genomsnittet för 1961-2015. Data kommer från SMHI.",
size=3, colour="gray30", hjust = 0) +
annotate("text", x = 5, y = 29, label = "Grafen är gjord av Richard Öhrvall 2017, @richardohrvall",
size=3, colour="gray30", hjust = 0)
p <- p +
annotate("segment", x = 311, xend = 308, y = -2.8, yend = -7, colour = "blue3") +
annotate("text", x = 306, y = -7, label = "Vi hade 4 dagar som var de", size=3, colour="blue3", hjust = 1) +
annotate("text", x = 306, y = -9, label = "kallaste sedan 1961", size=3, colour="blue3", hjust = 1) +
annotate("segment", x = 258, xend = 263, y = 18.2, yend = 25, colour = "firebrick3") +
annotate("text", x = 265, y = 28, label = "Vi hade 7 dagar som var de", size=3, colour="firebrick3", hjust = 0) +
annotate("text", x = 265, y = 26, label = "varmaste sedan 1961", size=3, colour="firebrick3", hjust = 0)
p <- p +
annotate("segment", x = 181, xend = 181, y = -15, yend = -5, colour = "wheat2", size=3) +
annotate("segment", x = 181, xend = 181, y = -12, yend = -8, colour = "wheat4", size=3) +
geom_line(data=legend_data, aes(x=x,y=y)) +
annotate("segment", x = 183, xend = 185, y = -11.9, yend = -11.9, colour = "wheat4", size=.5) +
annotate("segment", x = 183, xend = 185, y = -7.9, yend = -7.9, colour = "wheat4", size=.5) +
annotate("segment", x = 185, xend = 185, y = -7.9, yend = -11.9, colour = "wheat4", size=.5) +
annotate("text", x = 186, y = -9, label = "Normal medeltemperatur", size=2.5, colour="gray30", hjust = 0) +
annotate("text", x = 186, y = -11, label = "(innefattar hälften av åren)", size=2.5, colour="gray30", hjust = 0) +
annotate("text", x = 172, y = -10, label = "2016 års medeltemperatur", size=2.5, colour="gray30", hjust = 1) +
annotate("text", x = 184, y = -5.2, label = "Högsta medeltemperatur", size=2.5, colour="gray30", hjust = 0) +
annotate("text", x = 184, y = -14.3, label = "Lägsta medeltemperatur", size=2.5, colour="gray30", hjust = 0)Sedan är det bara att ta fram grafen. Visst ser den helt ok ut?
print(p)