Skip to content

Instantly share code, notes, and snippets.

@ilarischeinin
Created November 30, 2017 17:51
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save ilarischeinin/a91d697415403cba78524911d762af09 to your computer and use it in GitHub Desktop.
Save ilarischeinin/a91d697415403cba78524911d762af09 to your computer and use it in GitHub Desktop.
Plot snow depths for Helsinki winters
# snow.R: Plot snow depths for Helsinki winters
# Ilari Scheinin
# firstname.lastname@gmail.com
# MIT License
library(dplyr)
library(fmi)
library(ggplot2)
library(lubridate)
library(R.cache)
cached_fmi <- function(query, fmisid, date, sleep = 5L) {
fmisid <- as.integer(fmisid)
date <- format(date, "%F")
key <- list(query = query, fmisid = fmisid, date = date)
suffix <- paste0("::", query, "::", fmisid, "::", date)
cached <- loadCache(key = key, suffix = suffix, dirs = "fmi")
if (!is.null(cached)) {
message("Loaded cached ", suffix, ".")
return(cached)
}
request <- FMIWFSRequest$new(apiKey = apikey)
request$setParameters(request = "getFeature",
storedquery_id = query,
fmisid = fmisid,
starttime = paste0(date, "T00:00:00"),
endtime = paste0(date, "T23:59:59"))
client <- FMIWFSClient$new(request = request)
suppressMessages({
layers <- client$listLayers()
})
if (length(layers) > 0L) {
suppressMessages({
response <- client$getLayer(layer = layers[1L],
crs = "+proj = longlat +datum = WGS84", swapAxisOrder = TRUE,
parameters = list(splitListFields = TRUE))
})
data <- response@data %>%
tbl_df() %>%
transmute(
fmisid = fmisid,
time = ymd_hms(Time),
variable = ParameterName,
value = as.numeric(ParameterValue)
)
} else {
data <- data_frame(
fmisid = integer(),
time = as.POSIXct(character()),
variable = character(),
value = numeric())
}
saveCache(data, key = key, suffix = suffix, dirs = "fmi", compress = TRUE)
message("Downloaded ", suffix, ".")
Sys.sleep(sleep)
return(data)
}
apikey <- ""
query <- "fmi::observations::weather::daily::simple"
fmisid <- 100971L # Helsinki Kaisaniemi
dates <- seq.Date(from = as.Date("1911-01-01"), to = Sys.Date() - 1L,
by = "day")
depth <- rep(NA_real_, length.out = length(dates))
for (i in rev(seq_along(dates))) {
try({
depth[i] <- cached_fmi(query, fmisid, dates[i]) %>%
filter(variable == "snow") %>%
pull(value)
}, silent = TRUE)
}
snow <- data_frame(date = dates, depth = depth)
breaks <- dates[month(dates) == 7L & mday(dates) == 1L]
if (max(breaks) < Sys.Date()) {
breaks <- c(breaks, max(breaks) + years(1L))
}
labels <- year(breaks)[-length(breaks)]
snow$winter <- cut(dates, breaks = breaks, labels = labels) %>%
as.character() %>%
as.integer()
winters <- snow %>%
filter(!is.na(winter)) %>%
group_by(winter) %>%
summarise(
days = sum(depth > 0, na.rm = TRUE)
)
png("talvet.png", width = 297, height = 210, units = "mm", res = 150)
winters %>%
ggplot(aes(winter, days)) +
geom_line(size = 2L, color = "#18447e") +
annotate("text",
label = paste(format(max(dates), format = "%d.%m.%Y"), "asti"),
x = last(winters$winter), y = last(winters$days),
hjust = 1.1, size = rel(3)
) +
scale_x_continuous(labels = function(x) paste0(x, "-", x + 1L)) +
labs(
title = paste("Helsingin talvien lumisuus vuodesta", first(winters$winter)),
subtitle = paste(
"Talvikohtaiset kokonaismäärät vuorokausista,",
"joiden lumensyvyys Kaisaniemen havaintoasemalla > 0 cm."
),
x = "talvi",
y = "vuorokausien kokonaismäärä",
caption =
"Lähde: Ilmatieteen laitos, Tekijä: Ilari Scheinin, Lisenssi: CC-BY"
) +
theme_bw() +
theme(
plot.caption = element_text(size = rel(0.5))
)
dev.off()
# EOF
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment