sunpath.Rmd

---
title: "Plotting a sun-path diagram in R using data from solarpos"
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 Salzburg, Austria. Depending on your setup, you may have to add the "java" command and an absolute path.

```{zsh}
solarpos 47.795 13.047 2023 --timezone UTC --deltat --format=csv --headers position --step=600 > /tmp/sunpositions.csv
```

Now read that CSV and pick data for one day of each month. While we're at it, remove all observations where the sun isn't visible anyway (i.e. zenith angle is greater than 90°) and add a convenient "month" column. This will come in handy for the diagram.

```{r}
library(tidyverse)
library(lubridate, warn.conflicts = FALSE)

sunpath <- read_csv("/tmp/sunpositions.csv", 
                    show_col_types = FALSE) |> 
  filter(day(dateTime) == 21) |> 
  filter(zenith <= 90.0) |>
  mutate(month = month(dateTime, label=TRUE))

sunpath
```

Now plot a simple sun path diagram.

```{r}
plot <- ggplot() + 
  geom_line(data = sunpath, aes(x = azimuth, y = zenith, group = month, colour = month)) + 
  scale_y_continuous(limits = c(0, 90), breaks = seq(0, 90, by=15)) + 
  scale_x_continuous(limits = c(0, 360), breaks = seq(0, 359, by=15)) +
  labs(title = "Sun path for Salzburg, Austria in 2023") +
  coord_polar()

plot
```

Let's see if we can add the analemma lines for each hour as well.

```{r}
hours <- sunpath |>
  filter(minute(dateTime) == 0, second(dateTime) == 0) |>
  mutate(hour = hour(dateTime))

plot <- plot + 
  geom_path(data = hours, aes(x = azimuth, y = zenith, group = hour), linetype = "dashed")

plot
```

Let's see if we can add labels to those. This requires a bit of manual adjustment.


```{r}
hour_ends <- hours |> filter(month == "Jun")

plot <- plot +
  geom_label(data = hour_ends, nudge_y = -6, label.padding = unit(0.1, "lines"), size = 2, aes(x = azimuth, y = zenith, label = hour))

plot
```