Task

Incident radiation at the ground surface depends on the local slope and aspect of the terrain and on the geographical latitude, and on the date and time of day. Annual sums of incident solar radiation (disregarding atmospheric effects, hence considering top-of-atmosphere radiation) thus have a dependency on slope, aspect, and latitude. This quantity should affect vegetation functioning and structure and be expressed by patterns in structure and functioning along topographical gradients. The ratio of the slope and aspect-dependent incident radiation versus the flat-surface radiation (annual sums) measures the influence of the local terrain on radiation availability for vegetation. The task here is to quantify these quantities and this ratio.

Load packages and data

# load libraries
library(geodata)
## Loading required package: terra
## terra 1.8.50
library(terra)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks terra::extract()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(tidyterra)
## 
## Attaching package: 'tidyterra'
## 
## The following object is masked from 'package:stats':
## 
##     filter
source(here::here("R/helpers.R"))
source(here::here("R/calc_sw_in.R"))

# download SRTM data
# This stores file srtm_38_03.tif in tempdir()
geodata::elevation_3s(
    lat = 46.6756,
    lon = 7.85480,
    path = here::here("data/")
  )
## class       : SpatRaster 
## dimensions  : 6000, 6000, 1  (nrow, ncol, nlyr)
## resolution  : 0.0008333333, 0.0008333333  (x, y)
## extent      : 5, 10, 45, 50  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source      : srtm_38_03.tif 
## name        : srtm_38_03
# read the downloaded data
# use file.path() to combine
# a directory path with a filename
dem <- terra::rast(
  file.path(
    here::here("data/elevation/"),
    "srtm_38_03.tif"
    )
)

# reduce extent to example region (around Aletsch glacier)
# extent <- ext(8.2, 8.3, 46.6, 46.7)
# extent <- ext(7.8, 8.3, 46.3, 46.7)
extent <- ext(7.8, 7.9, 46.4, 46.5)
dem_crop <- crop(dem, extent)

Calculate slope and aspect

# Calculate slope and aspect
slope <- terrain(dem_crop, v = "slope", unit = "degrees")
aspect <- terrain(dem_crop, v = "aspect", unit = "degrees")

Plot the results

# Plot original elevation
ggplot() +
  tidyterra::geom_spatraster(data = dem_crop) +
  labs(title = "Elevation (m)") +
  scale_fill_viridis_c()

# Plot slope
ggplot() +
  tidyterra::geom_spatraster(data = slope) +
  labs(title = "Slope (°)") +
  scale_fill_viridis_c()

# Plot aspect
ggplot() +
  tidyterra::geom_spatraster(data = aspect) +
  labs(title = "Aspect (°)") +
  scale_fill_viridis_c()

Extract data from raster (takes some time).

df <- as.data.frame(dem_crop, xy = TRUE) |> 
  left_join(
    as.data.frame(slope, xy = TRUE),
    by = c("x", "y")
  ) |> 
  left_join(
    as.data.frame(aspect, xy = TRUE),
    by = c("x", "y")
  ) |> 
  as_tibble() |> 
  drop_na()

Model radiation

tictoc::tic()
out <- calc_sw_in(
    lat = seq(0, 90, by = 0.1),
    year = 2025
    )
tictoc::toc()

tictoc::tic()
out <- calc_sw_in(
    lat = seq(0, 90, by = 0.1),
    slope = 30,
    aspect = 180,
    year = 2025
    )
tictoc::toc()

This is very slow, but works.

tmp <- df |> 
  rowwise() |> 
  mutate(r_toa = calc_sw_in(
    lat = y,
    slope = slope,
    aspect = aspect,
    year = 2025
    ))
ggplot() +
  geom_raster(aes(x, y, fill = r_toa), data = tmp) +
  labs(title = "Incident solar radiation") +
  coord_sf() +
  scale_fill_viridis_c(option = "B")

Parallel computing

library(multidplyr)

ncores <- parallel::detectCores() - 1 

cl <- multidplyr::new_cluster(ncores) |>
      multidplyr::cluster_assign(
        calc_sw_in = calc_sw_in,   # make the function known for each core
        calc_sw_in_daily = calc_sw_in_daily,
        julian_day = julian_day,
        berger_tls = berger_tls,
        dcos = dcos,
        dsin = dsin
      )

tmp <- df |>
  dplyr::group_by(id = row_number()) |>
  tidyr::nest(
    input = c(
      y,
      slope,
      aspect)
  ) |> 
  multidplyr::partition(cl) |> 
  dplyr::mutate(r_toa = purrr::map(input, 
                                  ~calc_sw_in(
                                    lat = .x$y[[1]], 
                                    slope = .x$slope[[1]], 
                                    aspect = .x$aspect[[1]],  # seems best to "turn" aspect clockwise
                                    year = 2025
                                    )
  )) |>
  dplyr::collect() |>
  dplyr::ungroup() |>
  tidyr::unnest(c(r_toa, input))

Looks plausible.

ggplot() +
  geom_raster(aes(x, y, fill = r_toa / (365*60*60*24)), data = tmp) +
  labs(title = "Incident solar radiation") +
  coord_sf() +
  scale_fill_viridis_c(option = "B")

Plausibility check

I would expect south-facing slopes to receive the most light. This would look like a hillshade with a light source (45 degrees zenith angle) at 180 degrees (from the south):

# this doesn't seem to work
hillshade <- shade(
  slope = slope * pi / 180, 
  aspect = aspect * pi / 180, 
  angle = 45, 
  direction = 180
  )

# Plot aspect
ggplot() +
  tidyterra::geom_spatraster(data = hillshade) +
  labs(title = "Hillshade") +
  scale_fill_viridis_c(option = "B")

tmp |> 
  ggplot(aes(aspect, r_toa / (365*60*60*24), color = slope)) +
  geom_point()

Schematic terrain

df_test <- expand.grid(
  aspect = seq(1, 360, by = 1),
  slope = seq(0, 45, by = 5)) |> 
  as_tibble() |> 
  mutate(lat = 45) |> 
  rowwise() |> 
  mutate(r_toa = calc_sw_in(
    lat = lat,
    slope = slope,
    aspect = aspect,
    year = 2025
    ))
df_test |> 
  ggplot(aes(aspect, r_toa / (365*60*60*24), color = slope)) +
  geom_point()

–>

–> –> –>

–> –>

–> –>

–> –>

–>