Create vignette for brassens output analysis

R/analysis_plots.R

@@ -11,12 +11,12 @@ timeseries <- function(data, ts, te, var) {
   time <- site_id <- network <- NULL
   data_p <- data[which(between(data$time, ts, te)), ]
   min <- floor(quantile(data_p[[deparse(substitute(var))]],
-                        na.rm = TRUE,
-                        probs = 0.01
+    na.rm = TRUE,
+    probs = 0.0
   ))
   max <- ceiling(quantile(data_p[[deparse(substitute(var))]],
-                          na.rm = TRUE,
-                          probs = 0.997
+    na.rm = TRUE,
+    probs = 1
   ))
   p <- ggplot2::ggplot(data_p) +
     geom_line(aes(
@@ -31,6 +31,7 @@ timeseries <- function(data, ts, te, var) {
       date_minor_breaks = "1 hour",
       limits = c(ts, te)
     ) +
+    xlab("time (UTC)") +
     scale_y_continuous(
       limits = c(min, max),
       breaks = seq(min, max, 1)
@@ -64,12 +65,12 @@ hourly_boxplot <- function(data, ts, te, var) {
   network <- NULL
   data_p <- data[which(between(data$time, ts, te)), ]
   min <- floor(quantile(data_p[[deparse(substitute(var))]],
-                        na.rm = TRUE,
-                        probs = 0.01
+    na.rm = TRUE,
+    probs = 0.0
   ))
   max <- ceiling(quantile(data_p[[deparse(substitute(var))]],
-                          na.rm = TRUE,
-                          probs = 0.997
+    na.rm = TRUE,
+    probs = 1
   ))
   p <- ggplot2::ggplot(
     data = data_p,
@@ -111,16 +112,16 @@ hourly_boxplot <- function(data, ts, te, var) {
 #' @export
 #' @import ggplot2
 #' @author Eva Marques
-tile_ts <- function(data, ts, te, var, palname = "temp") {
+tile_ts <- function(data, ts, te, var, palname = "temp_ipcc") {
   time <- site_id <- NULL
   data_p <- data[which(between(data$time, ts, te)), ]
   min <- floor(quantile(data_p[[deparse(substitute(var))]],
-                        na.rm = TRUE,
-                        probs = 0.01
+    na.rm = TRUE,
+    probs = 0.0
   ))
   max <- ceiling(quantile(data_p[[deparse(substitute(var))]],
-                          na.rm = TRUE,
-                          probs = 0.997
+    na.rm = TRUE,
+    probs = 1
   ))
   p <- ggplot2::ggplot(
     data_p,
@@ -177,7 +178,7 @@ map_stations <- function(data,
                          datetime,
                          var,
                          imp,
-                         palname = "temp",
+                         palname = "temp_ipcc",
                          netw_shape = c(
                            "WU" = 17,
                            "PA" = 18,
@@ -186,18 +187,13 @@ map_stations <- function(data,
                          ),
                          title) {
   geometry <- network <- NULL
-  data_p <- data[which(between(data$time,
-                               datetime,
-                               datetime + lubridate::hours(1))), ]
-  min <- floor(quantile(data_p[[deparse(substitute(var))]],
-                        na.rm = TRUE,
-                        probs = 0.01
+  data_p <- data[which(data$time == datetime), ]
+  min <- floor(min(data_p[[deparse(substitute(var))]],
+    na.rm = TRUE
   ))
-  max <- ceiling(quantile(data_p[[deparse(substitute(var))]],
-                          na.rm = TRUE,
-                          probs = 0.997
+  max <- ceiling(max(data_p[[deparse(substitute(var))]],
+    na.rm = TRUE
   ))
-  cat("min: ", min, "max: ", max, "\n")
   p <- ggplot2::ggplot() +
     tidyterra::geom_spatraster(data = imp) +
     ggplot2::geom_sf(

R/clean_cws.R

@@ -57,7 +57,6 @@ clean_cws <- function(x) {
       data,
       m5_radius = 10000,
       m5_n_buddies = 5,
-      m5_keep_isolated = TRUE
     )
   }
   col_rm <- c(

R/format_observations.R

@@ -10,7 +10,7 @@
 #' @author Eva Marques
 #' @importFrom lubridate floor_date
 #' @importFrom stats median
-#' @importFrom dplyr group_by summarise ungroup
+#' @importFrom dplyr group_by summarise ungroup rename
 #' @importFrom data.table as.data.table
 #' @export
 summarize_hourly_temp <- function(x, time, temp, lat, lon) {
@@ -27,15 +27,15 @@ summarize_hourly_temp <- function(x, time, temp, lat, lon) {
   )
   hourly_avg <- x |>
     data.table::as.data.table() |>
-    dplyr::rename("lat" = lat) |>
-    dplyr::rename("lon" = lon) |>
-    dplyr::rename("time" = time) |>
-    dplyr::rename("temp" = temp)
+    rename("lat" = lat) |>
+    rename("lon" = lon) |>
+    rename("time" = time) |>
+    rename("temp" = temp)
   hourly_avg$time <- lubridate::floor_date(hourly_avg$time, "hour")
   hourly_avg <- hourly_avg |>
-    dplyr::group_by(lat, lon, time) |>
-    dplyr::summarise(temp = median(temp, na.rm = TRUE)) |>
-    dplyr::ungroup() |>
+    group_by(lat, lon, time) |>
+    summarise(temp = median(temp, na.rm = TRUE)) |>
+    ungroup() |>
     as.data.frame()
   # remove duplicates
   hourly_avg <- unique(hourly_avg)

R/utils.R

@@ -39,11 +39,41 @@ convert_temp <- function(x, from, to) {
 my_pal <- function(name) {
   stopifnot(
     "name should be one of \"temp\", \"sw\", \"reds\", \"prior\", or \"uhi\"" =
-      name %in% c("temp", "sw", "reds", "prior", "uhi")
+      name %in% c("temp", "sw", "reds", "prior", "uhi", "rh_ipcc", "temp_ipcc")
   )
   # -- define palettes
   if (name == "temp") {
     return(fields::tim.colors(n = 64, alpha = 1.0))
+  } else if (name == "rh_ipcc") {
+    pal_ipcc <- list(c(140, 81, 10),
+      c(191, 129, 45),
+      c(223, 194, 125),
+      c(246, 232, 195),
+      c(245, 245, 245),
+      c(199, 234, 229),
+      c(128, 205, 193),
+      c(53, 151, 143),
+      c(1, 102, 94)
+    ) |>
+      lapply(function(x) rgb(x[1], x[2], x[3], maxColorValue = 255)) |>
+      rev()
+    return(pal_ipcc)
+  } else if (name == "temp_ipcc") {
+    pal_ipcc <- list(c(103, 0, 31),
+      c(178, 24, 43),
+      c(214, 96, 77),
+      c(244, 165, 130),
+      c(253, 219, 199),
+      c(247, 247, 247),
+      c(209, 229, 240),
+      c(146, 197, 222),
+      c(67, 147, 195),
+      c(33, 102, 172),
+      c(5, 48, 97)
+    ) |>
+      lapply(function(x) rgb(x[1], x[2], x[3], maxColorValue = 255)) |>
+      rev()
+    return(pal_ipcc)
   } else if (name == "sw") {
     return(RColorBrewer::brewer.pal(10, "RdYlBu"))
   } else if (name == "reds") {

tests/testthat/test-utils.R

@@ -26,6 +26,8 @@ testthat::test_that("my_pal works well", {
   expect_no_error(pal3 <- my_pal("reds"))
   expect_no_error(pal4 <- my_pal("prior"))
   expect_no_error(pal5 <- my_pal("uhi"))
+  expect_no_error(pal5 <- my_pal("rh_ipcc"))
+  expect_no_error(pal5 <- my_pal("temp_ipcc"))
   expect_true(all(are_colors(pal1)))
   expect_true(all(are_colors(pal2)))
   expect_true(all(are_colors(pal3)))

vignettes/case_study_triangle_area.Rmd

@@ -0,0 +1,91 @@
+
+---
+title: "Brassens output analysis"
+format: html
+editor: visual
+author: "Eva Marques"
+date: "2021-07-23"
+---
+
+```{r}
+devtools::load_all()
+```
+
+### Some configuration parameters
+
+```{r}
+ts <- as.POSIXct("2021-07-22 00:00:00", tz = "UTC")
+te <- as.POSIXct("2021-07-23 23:59:59", tz = "UTC")
+timestamp <- as.POSIXct("2021-07-23 06:00:00", tz = "America/New_York")
+imp <- terra::rast("../case_study_triangle/imp_triangle.tif") |>
+  terra::project("EPSG:4326")
+```
+
+
+## Open datasets in .rds file
+
+```{r}
+cws_raw <- readRDS(paste0(
+  "../case_study_triangle/cws_triangle_raw_",
+  lubridate::date(ts),
+  "_",
+  lubridate::date(te),
+  ".rds"
+))
+cws_qc <- readRDS(paste0(
+  "../case_study_triangle/cws_triangle_cleaned_",
+  lubridate::date(ts),
+  "_",
+  lubridate::date(te),
+  ".rds"
+))
+cws_final <- readRDS(paste0(
+  "../case_study_triangle/cws_triangle_cleaned_calibrated_",
+  lubridate::date(ts),
+  "_",
+  lubridate::date(te),
+  ".rds"
+))
+```
+
+## Analysis
+
+### Timeseries at each step of brassens pipeline
+
+```{r}
+timeseries(cws_raw, ts = ts, te = te, var = temp)
+timeseries(cws_qc, ts = ts, te = te, var = temp)
+timeseries(cws_final, ts = ts, te = te, var = temp)
+```
+
+### Hourly boxplots at each step of brassens pipeline
+
+```{r}
+hourly_boxplot(cws_raw, ts = ts, te = te, temp)
+hourly_boxplot(cws_qc, ts = ts, te = te, temp)
+hourly_boxplot(cws_final, ts = ts, te = te, temp)
+```
+
+### Tile plots of temperature through time ar each step of brassens pipeline
+
+```{r}
+tile_ts(cws_raw, ts = ts, te = te, var = temp)
+tile_ts(cws_qc, ts = ts, te = te, var = temp)
+tile_ts(cws_final, ts = ts, te = te, var = temp)
+```
+
+### Maps of stations at each step of brassens pipeline
+
+```{r}
+map_stations(cws_raw, timestamp, var = temp, imp, title = "Raw data")
+map_stations(cws_qc,
+             timestamp,
+             var = temp,
+             imp,
+             title = "Cleaned with CrowdQC+")
+map_stations(cws_final,
+             timestamp,
+             var = temp,
+             imp,
+             title = "Cleaned with CrowdQC+ and calibrated with GHCNh")
+```