Plotting functions changes

CMIP6VisR.Rproj

@@ -1,5 +1,4 @@
 Version: 1.0
-ProjectId: 04cd0a55-412a-4631-97ae-e0438f635192
 
 RestoreWorkspace: Default
 SaveWorkspace: Default

DESCRIPTION

@@ -24,7 +24,11 @@ Imports:
     ggplot2,
     ggpubr,
     scales,
-    magrittr
+    magrittr,
+    sf,
+    stars,
+    tidyr,
+    units
 Encoding: UTF-8
 LazyData: true
 RoxygenNote: 7.3.2

NAMESPACE

@@ -4,24 +4,38 @@ export(cv_basin_daily_precip)
 export(cv_clip_basin)
 export(cv_find_point)
 export(cv_plot_TS)
+export(cv_plot_map)
 export(cv_plot_prob)
 export(cv_plot_season)
 export(cv_point_daily_precip)
 export(cv_zone_area_raster)
+export(plot_precip_map)
+export(plot_temp_map)
+import(dplyr)
 import(ggplot2)
-importFrom(dplyr,case_when)
+import(lubridate)
+import(scales)
 importFrom(dplyr,filter)
 importFrom(dplyr,group_by)
 importFrom(dplyr,mutate)
 importFrom(dplyr,n)
+importFrom(dplyr,right_join)
+importFrom(dplyr,rowwise)
+importFrom(dplyr,select)
 importFrom(dplyr,summarise)
-importFrom(dplyr,summarize)
+importFrom(dplyr,ungroup)
 importFrom(geosphere,distVincentyEllipsoid)
 importFrom(geosphere,distVincentySphere)
+importFrom(ggplot2,theme)
 importFrom(ggpubr,ggarrange)
 importFrom(lubridate,month)
-importFrom(magrittr,"%>%")
-importFrom(scales,label_comma)
+importFrom(lubridate,year)
+importFrom(sf,st_as_sf)
+importFrom(sf,st_drop_geometry)
+importFrom(stars,geom_stars)
+importFrom(stars,read_stars)
+importFrom(stars,st_rasterize)
+importFrom(stars,st_set_dimensions)
 importFrom(stats,na.omit)
 importFrom(stats,sd)
 importFrom(stringr,str_sub)
@@ -35,3 +49,5 @@ importFrom(terra,subset)
 importFrom(terra,trim)
 importFrom(terra,values)
 importFrom(terra,vect)
+importFrom(tidyr,pivot_longer)
+importFrom(units,drop_units)

R/cv_basin_daily_precip.R

@@ -26,8 +26,7 @@
 #' @importFrom stringr str_sub
 #' @importFrom terra crop
 #' @importFrom terra global
-#' @importFrom terra rast
-#' @returns Returns a data frame with 2 columns: `date` and `precipitation`. The
+#' @returns Reruns a data frame with 2 columns: `date` and `precipitation`. The
 #' `date` is a standard \R date over the interval 2015-01-01 to 2100-12-31, and the 
 #' `precipitation` is the basin mean value.
 #' @export

R/cv_clip_basin.R

@@ -39,7 +39,7 @@ cv_clip_basin <- function(za_rast, basin) {
 
   # Ensure zones is a vector before applying unique()
   zone_values <- terra::values(zones)  # Extract numeric values from raster
-  zone_values <- stats::na.omit(zone_values)  # Remove NA values
+  zone_values <- na.omit(zone_values)  # Remove NA values
   zone_values <- unique(zone_values)   # Get unique zone numbers
 
   # Apply unique() safely

R/cv_plot_TS.R

@@ -1,43 +1,59 @@
-#' Plot precipitation time series
+#' Plot Precipitation or Temperature Time Series
 #'
-#' Creates a time series plot from a data frame of precipitation data 
-#' as returned by `cv_basin_daily_precip()`. Basic statistics (mean, standard deviation, and proportion of zero values),
-#' are added to the title of the plot.
+#' This function takes a data frame containing a precipitation time series, 
+#' calculates basic statistics (mean, standard deviation, and proportion of zero values),
+#' and generates a ggplot visualization.
 #'
-#' @param data A data frame with two columns: \code{date} (date or datetime) and  
-#' \code{precipitation} (numeric) as returned by `cv_basin_daily_precip()`.
-#' @return A ggplot object displaying the precipitation time series with summary 
-#' statistics in the title. The returned plots look best when saved at the size 
-#' 16.5 x 6 cm You can easily change the font sizes using theme().
+#' @param data A data frame with two columns: "date" (date or datetime) and 
+#'             "precipitation" (numeric).
+#' @return A ggplot object displaying the precipitation time series with summary statistics in the title.
+#'          The returned plots look best when saved at the size 16.5 x 6 cm.
+#'          You can easily change the font sizes using theme().
 #' @examples
 #' cv_plot_TS(eg_TS)
 #' @import ggplot2
 #' @importFrom stats sd
-#' @seealso \code{\link{cv_basin_daily_precip}} \code{\link{cv_plot_prob}} \code{\link{cv_plot_season}}
 #' @export
-cv_plot_TS <- function(data) {
+cv_plot_TS <- function(data, variable = "precipitation") {
 
-  precipitation <- mean_val <- sd_val <-  p0  <- NULL
+  value <- precipitation <- mean_val <- sd_val <-  p0  <- NULL
 
 
-  # Ensure correct column names
-  if (!all(c("date", "precipitation") %in% colnames(data))) {
-    stop("Data frame must contain columns named 'date' and 'precipitation'.")
+  #search for the variable
+  var_lower <- tolower(variable)
+  if (!startsWith(var_lower, "p") && !startsWith(var_lower, "t")) {
+    stop("The variable argument must start with 'p' for precipitation or 't' for temperature.")
   }
 
+  #bring the variable to the required name
+  if (startsWith(var_lower, "p")) {
+    variable <- "precipitation"
+  } else if (startsWith(var_lower, "t")) {
+    variable <- "temperature"
+  } else {
+    stop("The variable argument must start with 'p' for precipitation or 't' for temperature.")
+  }
+
+
+  colnames(data)[2] <- "value"
   # Compute statistics
-  mean_val <- mean(data$precipitation, na.rm = TRUE)
-  sd_val <- sd(data$precipitation, na.rm = TRUE)
-  p0 <- sum(data$precipitation == 0, na.rm = TRUE) / nrow(data)
+  mean_val <- mean(data$value, na.rm = TRUE)
+  sd_val <- sd(data$value, na.rm = TRUE)
+  p0 <- sum(data$value == 0, na.rm = TRUE) / nrow(data)
 
   # Generate plot
-  plot <- ggplot(data, aes(x = date, y = precipitation), linewidth = 0.05) + 
+  plot <- ggplot(data, aes(x = date, y = value), linewidth = 0.05) + 
     geom_line(color = "gray20") + 
-    ylab("Precipitation") + 
+    ylab("Value") + 
     xlab("Date") +
-    ggtitle(paste0("mean = ", round(mean_val, 2), 
-                   ", sd = ", round(sd_val, 2), 
-                   ", and P0 = ", round(p0, 2))) + 
+    (if (variable == "precipitation") {
+      list(ggtitle(paste0("mean = ", round(mean_val, 2), 
+                          ", sd = ", round(sd_val, 2), 
+                          ", and P0 = ", round(p0, 2))))
+    } else {
+      list(ggtitle(paste0("mean = ", round(mean_val, 2), 
+                          ", sd = ", round(sd_val, 2))))
+    }) +
     theme(
       axis.title.x = element_text(size = 9, colour = "gray25"),
       axis.title.y = element_text(size = 9, colour = "gray25"),

R/cv_plot_map.R

@@ -0,0 +1,152 @@
+#' Plot Spatial Precipitation or Temperature Distribution
+#'
+#' This function takes a NetCDF file containing daily time series of precipitation or temperature,
+#' extracts spatial and temporal features, and plots one of three options: daily nonzero frequency,
+#' average monthly totals, or average annual maxima.
+#'
+#' @param nc_file A path to a NetCDF (.nc) file readable by the `stars` package.
+#' @param stat_type One of "daily", "monthly", or "annual" to indicate which type of statistic to visualize.
+#'                  "daily" provides the nonzero average of daily values during the whole period given in the NC file.
+#'                  "monthly" provides the average monthly totals (precipitation) or means (temperature).
+#'                  "annual" provides the annual maxima of daily values during the whole duration.
+#' @param variable A string starting with 'p' (precipitation) or 't' (temperature).
+#' @return A ggplot object showing the specified variable's distribution over space.
+#'         Note: We are using the default Viridis color scale which can be overwritten. 
+#'         The returned plots look best when saved at the size 14 x 12 cm.
+#'         You can easily change the font sizes using theme().
+#' @examples
+#' cv_plot_map("Pincher_ck.nc", stat_type = "monthly", variable = "temperature")
+#' cv_plot_map("Pincher_ck.nc", stat_type = "annual", variable = "precip")
+#' @importFrom ggplot2 theme
+#' @importFrom stars read_stars st_rasterize st_set_dimensions geom_stars
+#' @importFrom sf  st_as_sf st_drop_geometry
+#' @importFrom dplyr filter summarise mutate group_by select right_join rowwise ungroup
+#' @importFrom tidyr pivot_longer
+#' @importFrom lubridate year month
+#' @importFrom units drop_units
+#' @export
+cv_plot_map <- function(nc_file, stat_type = "daily", variable = "precipitation") {
+  stat_type <- tolower(stat_type[1])
+  var_lower <- tolower(variable)
+
+  if (!startsWith(var_lower, "p") && !startsWith(var_lower, "t")) {
+    stop("The variable argument must start with 'p' for precipitation or 't' for temperature.")
+  }
+
+  variable <- if (startsWith(var_lower, "p")) "precipitation" else "temperature"
+  valid_types <- c("daily", "monthly", "annual")
+
+  if (!stat_type %in% valid_types) {
+    stop("Invalid 'stat_type'. Please choose one of: 'daily', 'monthly', or 'annual'.")
+  }
+
+  mytheme <- theme(
+    legend.text = element_text(size = 7),
+    legend.title = element_text(size = 8),
+    axis.text.y = element_text(size = 5.5),
+    axis.text.x = element_text(size = 5.5, hjust = 0.4),
+    axis.title = element_text(size = 8),
+    legend.text.align = 0.5,
+    legend.title.align = 0.5,
+    legend.box.just = "center",
+    legend.justification = "center",
+    legend.position = 'bottom',
+    legend.key.size = unit(0.2, "cm"),
+    legend.key.width = unit(1.2, "cm"),
+    text = element_text(size = 7, color = gray(0.25))
+  )
+
+  nc_data <- read_stars(nc_file)
+  sf_obj <- st_as_sf(nc_data, as_points = FALSE, merge = FALSE)
+  sf_obj[] <- lapply(sf_obj, function(col) if (inherits(col, "units")) drop_units(col) else col)
+  cols <- names(sf_obj)[sapply(sf_obj, is.numeric)]
+
+  if (stat_type == "daily") {
+    daily_stat <- if (variable == "precipitation") {
+      sf_obj %>% rowwise() %>% mutate(result = mean(c_across(all_of(cols)) != 0, na.rm = TRUE))
+    } else {
+      sf_obj %>% rowwise() %>% mutate(result = mean(c_across(all_of(cols)), na.rm = TRUE))
+    }
+    daily_stat <- daily_stat %>% ungroup() %>% select(result, geometry)
+    stars_raster <- st_rasterize(daily_stat["result"])
+
+    fill_label <- if (variable == "precipitation") "Nonzero Avg Daily Precipitation (mm)" else "Nonzero Avg Daily Temperature (°C)"
+
+    p <- ggplot() +
+      geom_stars(data = stars_raster) +
+      coord_equal() +
+      scale_fill_viridis_c(
+        fill_label,
+        direction = if (variable == "precipitation") -1 else 1,
+        trans = scales::pseudo_log_trans(),
+        na.value = "transparent"
+      ) +
+      theme_minimal() + xlab("Longitude") + ylab("Latitude") + mytheme +
+      guides(fill = guide_colourbar(title.position = "top", title.hjust = 0.5))
+
+    print(p)
+
+  } else {
+    sf_long <- sf_obj %>%
+      st_drop_geometry() %>%
+      mutate(id = row_number()) %>%
+      pivot_longer(cols = matches("^\\d{4}-\\d{2}-\\d{2}$"), names_to = "date", values_to = "value") %>%
+      mutate(date = as.Date(date), year = year(date))
+
+    if (stat_type == "monthly") {
+      sf_long <- sf_long %>% mutate(month = month(date, label = TRUE, abbr = TRUE))
+
+      monthly_data <- if (variable == "precipitation") {
+        sf_long %>% group_by(id, year, month) %>% summarise(total = sum(value, na.rm = TRUE), .groups = "drop") %>%
+          group_by(id, month) %>% summarise(avg_total = mean(total, na.rm = TRUE), .groups = "drop")
+      } else {
+        sf_long %>% group_by(id, month) %>% summarise(avg_total = mean(value, na.rm = TRUE), .groups = "drop")
+      }
+
+      sf_avg_monthly <- sf_obj %>% mutate(id = row_number()) %>% select(id, geometry) %>% right_join(monthly_data, by = "id") %>% st_as_sf()
+      raster_list <- sf_avg_monthly %>% split(.$month) %>% lapply(function(sf_month) st_rasterize(sf_month["avg_total"]))
+
+      stars_monthly <- do.call(c, c(raster_list, along = "month"))
+      stars_monthly <- st_set_dimensions(stars_monthly, "month", values = levels(factor(sf_avg_monthly$month, levels = month.abb)))
+
+      fill_label <- if (variable == "precipitation") "Avg Monthly Precipitation (mm)" else "Avg Monthly Temperature (°C)"
+
+      p <- ggplot() +
+        geom_stars(data = stars_monthly) +
+        coord_equal() +
+        facet_wrap(~month) +
+        scale_fill_viridis_c(
+          fill_label,
+          direction = if (variable == "precipitation") -1 else 1,
+          trans = scales::pseudo_log_trans(),
+          na.value = "transparent"
+        ) +
+        theme_minimal() + xlab("Longitude") + ylab("Latitude") + mytheme +
+        theme(panel.spacing.x = unit(1, "lines")) +
+        guides(fill = guide_colourbar(title.position = "top", title.hjust = 0.5))
+
+      print(p)
+
+    } else if (stat_type == "annual") {
+      annual_max <- sf_long %>% group_by(id, year) %>% summarise(max_val = max(value, na.rm = TRUE), .groups = "drop")
+      avg_annual_max <- annual_max %>% group_by(id) %>% summarise(result = mean(max_val, na.rm = TRUE), .groups = "drop")
+      sf_annual_max <- sf_obj %>% mutate(id = row_number()) %>% select(id, geometry) %>% inner_join(avg_annual_max, by = "id") %>% st_as_sf()
+      stars_max <- st_rasterize(sf_annual_max["result"])
+
+      fill_label <- if (variable == "precipitation") "Annual Max of Daily Precip (mm)" else "Annual Max of Daily Temp (°C)"
+
+      p <- ggplot() +
+        geom_stars(data = stars_max) +
+        coord_equal() +
+        scale_fill_viridis_c(
+          fill_label,
+          direction = if (variable == "precipitation") -1 else 1,
+          na.value = "transparent"
+        ) +
+        theme_minimal() + xlab("Longitude") + ylab("Latitude") + mytheme +
+        guides(fill = guide_colourbar(title.position = "top", title.hjust = 0.5))
+
+      print(p)
+    }
+  }
+}