Merge pull request #29 from uace-azmet/sol-rad

Sol rad

.gitignore

@@ -2,3 +2,4 @@
 .Rhistory
 .RData
 .Ruserdata
+.Renviron

.renvignore

@@ -0,0 +1 @@
+notes/

app/R/calc_sol_rad_theoretical.R

@@ -0,0 +1,161 @@
+# library(TrenchR)
+
+# Wrapper to sum direct, diffuse, and reflected solar radiation output by `solar_radiation()`
+solar_radiation_total <- function(doy, psi, tau, elev, rho) {
+  purrr::pmap(list(doy, psi, tau, elev, rho), solar_radiation) |> 
+    purrr::map_dbl(sum) 
+}
+
+
+
+#' Calculate theoretical maximum solar radiation
+#'
+#' @param data data returned by either `az_daily()` or `az_hourly()`
+#' @param freq indicator of whether `data` is hourly or daily (default is hourly)
+#' @param tau passed to `TrenchR::solar_radiation()` (via
+#'   `solar_radiation_total()`).  Atmospheric sensitivity from 0-1.  While 1 is
+#'   the theoretical max, 0.7 seems like the max observed globaly.  E.g. in
+#'   https://www.mdpi.com/2072-4292/13/9/1716
+#' @param rho passed to `TrenchR::solar_radiation()` (via
+#'   `solar_radiation_total()`). Albedo ranging from 0-1. Typical albedo values:
+#'   desert sand, 0.4; concrete, 0.55; bare soil, 0.17; asphalt, 0.04 (values
+#'   from https://en.wikipedia.org/wiki/Albedo).
+#'
+#' @return `data`, but with a `sol_rad_est` column indicating the estimated
+#'   theoretical maximum solar radiation for that date or datetime
+#'   
+calc_sol_rad_theoretical <- function(data, freq = c("hourly", "daily"), tau = 0.7, rho = 0.4) {
+  # hourly and daily data have slightly different formats, so must deal with that
+  freq <- match.arg(freq)
+  if (freq == "hourly") {
+    date_col <- sym("date_datetime")
+  } else {
+    date_col <- sym("datetime")
+  }
+
+  # Join location data ------------------------------------------------------
+  data_joined <- 
+    left_join(data, station_info, by = join_by(meta_station_id, meta_station_name)) |> 
+    select(meta_station_id, meta_station_name, {{ date_col }}, latitude, longitude, elev_m, sol_rad_total)
+
+
+  # Expand grid -------------------------------------------------------------
+  # For hourly data, to get better resolution calculate instantaneous solar
+  # radiation every few minutes and then sum. For daily data, calculate hourly
+  # and sum.
+
+  if (freq == "hourly") {
+
+    # NOTE: this calculation of midpoint is NOT ROBUST. Decimal values will get
+    # coerced to integers and then be wrong. Edit with caution!
+    n_segments <- 10L
+    start <- 0 + (1/2 * 1/n_segments)
+    midpts <- seq(start, by = 1/n_segments, length.out = n_segments) * 60
+
+    data_expand <-
+      data_joined |> 
+      expand_grid(
+        minute = midpts
+      ) |>
+      mutate(grid_datetime = make_datetime(
+        year = year(date_datetime),
+        month = month(date_datetime),
+        day = day(date_datetime),
+        hour = hour(date_datetime),
+        min = as.integer(minute)
+      ))
+  } else {
+    data_expand <- 
+      data_joined |> 
+      expand_grid(
+        hour = 0:23
+      ) |> 
+      mutate(grid_datetime = make_datetime(
+        year = year(datetime),
+        month = month(datetime),
+        day = day(datetime),
+        hour = hour,
+        min = 30 #estimate at midpoint of hour
+      ))
+  }
+
+
+  # Calculate instantaneous solar radiation ---------------------------------
+  data_sol <- 
+    data_expand |> 
+    #Zenith angle
+    mutate(psi = zenith_angle(
+      doy = yday({{ date_col }}),
+      lat = latitude,
+      lon = longitude,
+      hour = as_hms(grid_datetime) |> as.duration() |> as.numeric("hours"),
+      offset = -7
+    )) |> 
+    #Total solar radiation
+    mutate(
+      sol_rad_est = solar_radiation_total(
+        doy = yday({{ date_col }}),
+        psi = psi * pi / 180,
+        tau = tau,
+        elev = elev_m,
+        rho = rho 
+      ) # in W/m^2
+    ) 
+
+
+  # Integrate instantaneous to total ----------------------------------------
+  # For hourly, need to multiply by 1/n() (∆x).  For daily, just sum hourly for each date
+
+  if (freq == "hourly") {
+    data_summed <- 
+      data_sol |> 
+      summarize(
+        sol_rad_est = sum(sol_rad_est * (1/n()), na.rm = TRUE),
+        .by = c(meta_station_id, meta_station_name, date_datetime)
+      ) |> 
+      # Sensor totals are for the previous hour, so lag estimates to match
+      mutate(sol_rad_est = lag(sol_rad_est), .by = meta_station_id)
+  } else {
+    data_summed <-
+      data_sol |> 
+      summarize(
+        sol_rad_est = sum(sol_rad_est, na.rm = TRUE),
+        .by = c(meta_station_id, meta_station_name, datetime)
+      )
+  }
+
+  data_summed |> 
+    mutate(
+      #convert to same units as azmet data: 1 W*hr/m^2 = 0.0036 MJ/m^2,
+      sol_rad_est = sol_rad_est * 0.0036
+    ) |> 
+    right_join(data, by = join_by(meta_station_id, meta_station_name, {{ date_col }}))
+}
+
+# hourly <- az_hourly(start = "2023-05-29 00", end = "2023-05-30 23")
+# hourly_calced <- calc_sol_rad_theoretical(hourly)
+# nrow(hourly) == nrow(hourly_calced)
+# 
+# hourly_calced |>
+#   filter(!lte(sol_rad_total, sol_rad_est, tol = 0.5)) |>
+#   select(date_datetime, sol_rad_est, sol_rad_total)
+# 
+# hourly_calced |>
+#   filter(meta_station_id %in% c("az01")) |>
+# ggplot(aes(x = date_datetime)) +
+#   geom_line(aes(y = sol_rad_total, color = "Observed")) +
+#   geom_line(aes(y = sol_rad_est, color = "Theoretical max")) +
+#   facet_wrap(~meta_station_id)
+# 
+# daily <- az_daily(start = "2022-10-10", end = "2023-10-10")
+# daily_calced <- calc_sol_rad_theoretical(daily, "daily")
+# daily_calced |>
+#   filter(!lte(sol_rad_total, sol_rad_est)) |>
+#   select(datetime, sol_rad_est, sol_rad_total)
+# 
+# daily_calced |>
+#   filter(meta_station_id %in% c("az01")) |>
+#   ggplot(aes(x = datetime)) +
+#   geom_line(aes(y = sol_rad_total, color = "Observed")) +
+#   geom_line(aes(y = sol_rad_est, color = "Theoretical max")) +
+#   facet_wrap(~meta_station_id)

app/R/check_daily.R

@@ -11,6 +11,8 @@
 check_daily <- function(daily) {
   report <- data.validator::data_validation_report()
 
+  daily <- calc_sol_rad_theoretical(daily, freq = "daily") 
+
   data.validator::validate(daily, name = "Daily Data") |>
     # Internal consistency checks from 'NWS (1994) TSP 88-21-R2':
     data.validator::validate_if(gte(temp_air_meanC, dwpt_mean, na_pass = TRUE),
@@ -32,6 +34,8 @@ check_daily <- function(daily) {
       ),
       description = "`wind_spd_mean_mps` ≤ `wind_spd_max_mps`"
     ) |>
+    validate_if(lte(sol_rad_total, sol_rad_est, na_pass = TRUE),
+                "`sol_rad_total` ≤ theoretical max") |> 
     data.validator::validate_if(
       btwn(
         relative_humidity_mean,

app/R/check_hourly.R

@@ -9,6 +9,7 @@
 #' check_hourly(hourly)
 check_hourly <- function(hourly) {
   hourly <- hourly |>
+    calc_sol_rad_theoretical() |> 
     dplyr::group_by(meta_station_name) |>
     dplyr::arrange(dplyr::desc(date_datetime)) |> #make sure arranged in datetime order
     # Create a few new columns for temporal consistency checks from 'NWS (1994) TSP 88-21-R2':
@@ -43,13 +44,15 @@ check_hourly <- function(hourly) {
   report <- data.validator::data_validation_report()
   data.validator::validate(hourly, name = "Hourly Data") |>
     data.validator::validate_if(# within rounding error
-      gte(temp_airC, dwpt - 0.2, na_pass = TRUE),
+      gte(temp_airC, na_pass = TRUE, tol = 0.2),
       "`temp_airC` ≥ `dwpt`") |>
     data.validator::validate_if(lte(wind_spd_mps, wind_spd_max_mps, na_pass = TRUE),
                                 "`wind_spd` ≤ `wind_spd_max`") |>
     validate_if(temp_airC_delta, "|∆`temp_airC`| ≤ 19.4") |> 
     validate_if(relative_humidity_delta, "|∆`relative_humidity`| ≤ 50") |>
     validate_if(wind_spd_delta, "|∆`wind_spd_mps`| ≤ 10.3") |> 
+    validate_if(lte(sol_rad_total, sol_rad_est, na_pass = TRUE, tol = 0.5),
+                "`sol_rad_total` ≤ theoretical max") |> 
     #TODO would be nice if CSV would contain all 14+ hours in a row maybe?
     validate_if(sol_rad_total_14 | is.na(sol_rad_total_14),
                 "`sol_rad_total` not < 0.1 for more than 14 hrs") |> 

app/R/helpers.R

@@ -1,14 +1,18 @@
 #' @param na_pass logical; do NAs count as passing the validation step?
-gte <- function(x,y, na_pass = FALSE) {
-  res <- x >= y
+#' @param tol tolerance; passed on to `all.equal()` to evaluate the "equal" part
+#'   of greater than or equal to
+gte <- function(x,y, na_pass = FALSE, tol = sqrt(.Machine$double.eps)) {
+  eq <- purrr::map2_lgl(x, y, \(x,y) isTRUE(all.equal(x,y, tolerance = tol, scale = 1)))
+  res <- x > y | eq
   if(isTRUE(na_pass)) {
     res <- ifelse(is.na(res), TRUE, res)
   }
   res
 }
 
-lte <- function(x,y, na_pass = FALSE) {
-  res <- x <= y
+lte <- function(x,y, na_pass = FALSE, tol = sqrt(.Machine$double.eps)) {
+  eq <- purrr::map2_lgl(x, y, \(x,y) isTRUE(all.equal(x,y, tolerance = tol, scale = 1)))
+  res <- x < y | eq
   if(isTRUE(na_pass)) {
     res <- ifelse(is.na(res), TRUE, res)
   }

app/app.R

@@ -15,6 +15,8 @@ library(plotly)
 library(slider)
 library(units)
 library(patchwork)
+library(hms)
+library(TrenchR)
 
 # Everything in R/ should be sourced automatically, but I found I needed to add
 # these to get it to work.  Not sure why.

app/testdata_hourly.csv

@@ -26,7 +26,7 @@ meta_bat_volt,meta_needs_review,meta_station_id,meta_station_name,meta_version,d
 12.68,0,az02,Yuma Valley,1,2023-07-02T07:00:00Z,183,0700,2023,9.3,48.7,0,0,22.8,73,0,0,37,0,9.8,25.2,77.4,27.1,80.8,26.5,79.7,1.17,2.05,0.7,0.3,0,0,2023-07-02T06:00:10-07:00,18,1.8,0.8,0,0,317,3,0,0
 13.01,0,az02,Yuma Valley,1,2023-07-02T08:00:00Z,183,0800,2023,6.2,43.2,0.2,0.01,26,78.7,0,0,22,0,25.33,30.2,86.4,26.9,80.4,26.5,79.7,0.95,3.33,3.8,1.7,0.9,0.4,2023-07-02T08:00:00-07:00,189,5.1,2.3,0.9,0.4,159,18,0.3,0.7
 13.61,0,az02,Yuma Valley,1,2023-07-02T09:00:00Z,183,0900,2023,7.3,45.1,0.5,0.02,27.2,81,0,0,22,0,41.59,32,89.6,27.2,81,26.5,79.7,1.02,3.73,11.4,5.1,8.1,3.6,2023-07-02T08:43:30-07:00,177,14.3,6.4,8.1,3.6,166,14,3.5,7.8
-12.84,0,az02,Yuma Valley,1,2023-07-02T21:00:00Z,183,2100,2023,11.2,52.1,0.1,0,29.9,85.7,0,0,26,0,0,33.2,91.8,33.8,92.8,26.4,79.5,1.33,3.76,4,1.8,2,0.9,2023-07-02T20:41:20-07:00,269,5.6,2.5,2,0.9,241,21,0.8,1.8
+12.84,0,az02,Yuma Valley,1,2023-07-02T21:00:00Z,183,2100,2023,11.2,52.1,0.1,0,29.9,85.7,0,0,26,2.3,0,33.2,91.8,33.8,92.8,26.4,79.5,1.33,3.76,4,1.8,2,0.9,2023-07-02T20:41:20-07:00,269,5.6,2.5,2,0.9,241,21,0.8,1.8
 12.8,0,az02,Yuma Valley,1,2023-07-02T22:00:00Z,183,2200,2023,12.2,53.9,0.1,0,28,82.4,0,0,32,0,0,31.1,88,33,91.4,26.4,79.5,1.42,3.08,3.1,1.4,1.1,0.5,2023-07-02T21:02:00-07:00,223,4,1.8,1.1,0.5,237,13,0.5,1.1
 12.76,0,az02,Yuma Valley,1,2023-07-02T23:00:00Z,183,2300,2023,13.1,55.6,0,0,24.6,76.3,0,0,41,0,0,27.6,81.7,32.2,90,26.4,79.5,1.51,2.18,2,0.9,0,0,2023-07-02T22:16:40-07:00,87,2.7,1.2,0,0,82,3,0,0
 12.76,0,az04,Safford,1,2023-07-01T21:00:00Z,182,2100,2023,-0.5,31.1,0.2,0.01,28.7,83.6,0,0,11,0,0,33.5,92.3,34.3,93.7,29.4,84.9,0.59,4.57,10.1,4.5,7.6,3.4,2023-07-01T20:54:30-07:00,315,12.8,5.7,0,3.4,314,13,3.3,7.4