-
Notifications
You must be signed in to change notification settings - Fork 0
/
carpet.Rmd
135 lines (109 loc) · 3.94 KB
/
carpet.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
---
title: "A temporal raster plot (carpet diagram) of annual sunlight based on solarpos data"
editor_options:
chunk_output_type: inline
output: github_document
---
First, run [solarpos](https://github.com/klausbrunner/solarpos) and capture its output in a CSV file. We're getting position data for the entire year 2023 in Oulu, Finland. Depending on your setup, you may have to add the "java" command and an absolute path.
```{zsh}
solarpos 65.013 25.473 2023 --timezone UTC --deltat --format=csv --headers position --step=180 > /tmp/sunpositions.csv
```
Now read that CSV file, dropping what we don't need and regrouping the data a bit for convenient plotting.
```{r}
library(tidyverse)
library(lubridate, warn.conflicts = FALSE)
sunpath <- read_csv(
"/tmp/sunpositions.csv",
show_col_types = FALSE
) |>
mutate(
date = date(dateTime),
time = hour(dateTime) + minute(dateTime) / 60,
azimuth = NULL,
daytime = NULL
)
sunpath
```
With the nicely prepared data, plotting should be straightforward now. Getting the axes and labelling right is a bit of a hassle, though.
```{r}
library(scales, warn.conflicts = FALSE)
plot <- ggplot(sunpath, aes(date, time)) +
geom_raster(aes(fill = zenith)) +
scale_y_continuous(expand = c(0, 0), breaks = seq(0, 24, 3)) +
scale_x_date(
labels = date_format("%m"),
breaks = breaks_width("month"),
expand = c(0, 0)
)
plot
```
It's a start, but the default colour palette clearly isn't a good choice.
```{r}
plot <- plot + scale_fill_continuous(type = "viridis", direction = -1)
plot
```
Aesthetically better, but still hard to make sense of. It really seems we need a custom colour scheme with irregular breaks, and I haven't had much luck defining those in R. To make this easier and explicit rather than relying on ggplot's colour mapping magic, let's change the data from continuous to categorical first, and use a manual colour mapping.
```{r}
sunpath <- sunpath |>
mutate(light = cut(
zenith,
labels = c(
"day",
"golden hour",
"civil twl",
"nautical twl",
"astronomical twl",
"night"
),
breaks = c(0, 84, 90, 96, 102, 108, 180)
))
sunpath
plot <- ggplot(sunpath, aes(date, time)) +
geom_raster(aes(fill = light), alpha = 0.85) +
scale_y_continuous(expand = c(0, 0), breaks = seq(0, 24, 3)) +
labs(title = "2023 Sunlight hours in Oulu, Finland") +
scale_x_date(
labels = date_format("%m"),
breaks = breaks_width("month"),
expand = c(0, 0)
) +
scale_fill_manual(
values = c(
"day" = "lightyellow",
"golden hour" = "gold",
"civil twl" = "lightblue",
"nautical twl" = "blue",
"astronomical twl" = "blue4",
"night" = "black"
)
)
plot
```
Doesn't have the same mysterious vibe as the previous ones, but it's a bit more useful now: we can clearly see those very long Northern summer days and depressingly long winter nights, peaking on 21 June and 21 December, respectively. Of course, don't forget all times are UTC. You can also tell solarpos-cli to use the proper local timezone (Europe/Helsinki).
One of the problems here is that the contours aren't smooth. Which is likely because we're using a raster plot: basically a point matrix. What if we used an actual contour plot after all?
```{r}
plot <- ggplot(sunpath, aes(date, time)) +
geom_contour_filled(aes(z = zenith),
breaks = c(0, 84, 90, 96, 102, 108, 180),
alpha = 0.85) +
scale_fill_manual(
name = "light",
values = c("lightyellow", "gold", "lightblue", "blue", "blue4", "black"),
labels = c(
"daylight",
"golden hour",
"civil twilight",
"nautical twilight",
"astronomical twilight",
"night"
)
) +
scale_y_continuous(expand = c(0, 0), breaks = seq(0, 24, 3)) +
labs(title = "2023 Sunlight hours in Oulu, Finland") +
scale_x_date(
labels = date_format("%m"),
breaks = breaks_width("month"),
expand = c(0, 0)
)
plot
`````