diff --git a/CMIP6VisR.Rproj b/CMIP6VisR.Rproj
index 2d53c6a..4ba5539 100644
--- a/CMIP6VisR.Rproj
+++ b/CMIP6VisR.Rproj
@@ -1,5 +1,4 @@
 Version: 1.0
-ProjectId: 04cd0a55-412a-4631-97ae-e0438f635192
 
 RestoreWorkspace: Default
 SaveWorkspace: Default
diff --git a/DESCRIPTION b/DESCRIPTION
index 96ddb39..838f9c1 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -24,7 +24,11 @@ Imports:
     ggplot2,
     ggpubr,
     scales,
-    magrittr
+    magrittr,
+    sf,
+    stars,
+    tidyr,
+    units
 Encoding: UTF-8
 LazyData: true
 RoxygenNote: 7.3.2
diff --git a/NAMESPACE b/NAMESPACE
index 26ea89f..610ce26 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/R/cv_basin_daily_precip.R b/R/cv_basin_daily_precip.R
index b62ffc7..800688b 100644
--- a/R/cv_basin_daily_precip.R
+++ b/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
diff --git a/R/cv_clip_basin.R b/R/cv_clip_basin.R
index 8476d3d..48a8f7b 100644
--- a/R/cv_clip_basin.R
+++ b/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
diff --git a/R/cv_plot_TS.R b/R/cv_plot_TS.R
index f4ac156..f2c7707 100644
--- a/R/cv_plot_TS.R
+++ b/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"),
diff --git a/R/cv_plot_map.R b/R/cv_plot_map.R
new file mode 100644
index 0000000..f3308ca
--- /dev/null
+++ b/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)
+    }
+  }
+}
diff --git a/R/cv_plot_prob.R b/R/cv_plot_prob.R
index 237b689..308c1e6 100644
--- a/R/cv_plot_prob.R
+++ b/R/cv_plot_prob.R
@@ -1,46 +1,60 @@
-#' Plot probability of exceedance for precipitation
+#' Plot Probability of Exceedance for Precipitation
 #'
-#' Plots the exceedance probabilities of non-zero values in a precipitation time 
-#' series as returned by `cv_basin_daily_precip()`, vs the precipitation value. 
-#' Uses a logarithmic scale for the y-axis.
+#' This function takes a data frame containing a precipitation time series, 
+#' filters out zero precipitation values, calculates the probability of exceedance, 
+#' and generates a scatter plot with a logarithmic y-axis.
 #'
-#' @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 probability of exceedance of nonzero 
-#' precipitation. The returned plots look best when saved at the size 
-#' 16.5 x 12 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 probability of exceedance of nonzero precipitation.
+#'          The returned plots look best when saved at the size 16.5 x 12 cm.
+#'          You can easily change the font sizes using theme().
 #' @examples
 #' cv_plot_prob(eg_TS)
 #' @import ggplot2
-#' @importFrom dplyr filter
-#' @importFrom magrittr %>%
-#' @importFrom scales label_comma
+#' @import dplyr
+#' @import lubridate
+#' @import scales
 #' @export
-#' @seealso \code{\link{cv_basin_daily_precip}} \code{\link{cv_plot_TS}} \code{\link{cv_plot_season}}
-
-cv_plot_prob <- function(data) {
+cv_plot_prob <- function(data, variable = "precipitation") {
   
-  precipitation <- p_e <- NULL
+  value <- precipitation <- p_e <- 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.")
   }
   
-  # Filter out zero precipitation values
-  data_nz <- data %>% filter(precipitation != 0)
+  #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"
   
-  # Calculate probability of exceedance using the 
-  data_nz$rank <- rank(data_nz$precipitation)
+  
+  # Filter out zero precipitation values
+  data_nz <- if (variable == "precipitation") {
+    filter(data, value != 0)
+  } else {
+    data
+  }
+  # Calculate probability of exceedance
+  data_nz$rank <- rank(data_nz$value)
   data_nz$p_ne <- data_nz$rank / (nrow(data_nz) + 1)
   data_nz$p_e <- 1 - data_nz$p_ne
   
   # Generate plot
   plot <- ggplot() +
-    geom_point(data = data_nz, aes(x = precipitation, y = p_e), size = 0.5) + 
+    geom_point(data = data_nz, aes(x = value, y = p_e), size = 0.5) + 
     scale_y_log10(labels = scales::label_comma()) + 
     ylab("Probability of exceedance") + 
-    xlab("Precipitation") +
+    xlab("Value") +
     theme(
       axis.title.x = element_text(size = 9, colour = "gray25"),
       axis.title.y = element_text(size = 9, colour = "gray25"),
diff --git a/R/cv_plot_season.R b/R/cv_plot_season.R
index 915755d..2811470 100644
--- a/R/cv_plot_season.R
+++ b/R/cv_plot_season.R
@@ -1,35 +1,43 @@
-#' Plot seasonal precipitation distribution
-#' @description
-#' Creates seasonal plots of a precipitation time series, as returned by `cv_basin_daily_precip()`. Assigns seasons based on the month, 
-#' computes the proportion of zero precipitation values (P0) and seasonal means for nonzero precipitation, 
-#' and generates a violin plot showing the seasonal distribution of precipitation.
+#' Plot Seasonal Distribution of a Hydroclimatic Variable
 #'
-#' @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 seasonal distribution of nonzero precipitation, 
-#' with mean values and P0 labelled. The returned plots look best when saved at 
-#' the size 16.5 x 14 cm. You can easily change the font sizes using theme().
+#' This function takes a data frame with a hydroclimatic time series (precipitation or temperature), assigns seasons based on months, 
+#' computes the proportion of zero precipitation values (P0) and seasonal means for nonzero values, 
+#' and generates a violin plot showing the seasonal distribution.
+#'
+#' @param data A data frame with two columns: "date" (date or datetime) and a variable column (e.g., "precipitation" or "temperature").
+#' @param variable A character string, either "precipitation" or "temperature". Determines plotting behavior and calculation of P0.
+#' @return A ggplot object showing seasonal distribution with summary values and optional P0.
 #' @examples
-#' cv_plot_season(eg_TS)
+#' cv_plot_season(data, variable = "precipitation")
 #' @import ggplot2
-#' @importFrom dplyr filter summarize mutate group_by case_when summarise n
-#' @importFrom lubridate month
-#' @importFrom magrittr %>%
+#' @importFrom dplyr filter mutate group_by summarise n
+#' @importFrom lubridate month year
 #' @importFrom ggpubr ggarrange
-#' @seealso \code{\link{cv_basin_daily_precip}} \code{\link{cv_plot_TS}} \code{\link{cv_plot_prob}}
-
 #' @export
-#' 
-cv_plot_season <- function(data) {
-  precipitation <- season <- mean_val <- p0 <- tot_val <- NULL
+cv_plot_season <- function(data, variable = "precipitation") {
+  
+  value <- precipitation <- season <- mean_val <- p0 <- tot_val <- 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.")
   }
   
-  # Extract month and assign seasons
+  #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"
   data$month <- month(data$date)
+  data$year <- year(data$date)
+  
   data <- data %>% 
     mutate(season = case_when(
       month %in% c(9, 10, 11) ~ "Fall",
@@ -38,37 +46,56 @@ cv_plot_season <- function(data) {
       month %in% c(6, 7, 8) ~ "Summer"
     ))
   
-  # Compute proportion of zero precipitation (P0)
-  p0_seas <- data %>% 
-    group_by(season) %>% 
-    summarise(p0 = sum(precipitation == 0, na.rm = TRUE) / n())
+  if (variable == "precipitation") {
+    p0_seas1 <- data %>%
+      group_by(season, year) %>%
+      summarise(p0 = sum(value == 0, na.rm = TRUE) / n(), .groups = "drop")
+    
+    p0_seas <- p0_seas1 %>%
+      group_by(season) %>%
+      summarise(p0 = mean(p0), .groups = "drop")
+  }
   
-  # Filter out zero precipitation values
-  data_nz <- data %>% filter(precipitation != 0)
+  data_nz <- if (variable == "precipitation") {
+    filter(data, value != 0)
+  } else {
+    data
+  }
   
-  # Compute seasonal mean for nonzero precipitation
-  mean_seas <- data_nz %>% 
-    group_by(season) %>% 
-    summarise(mean_val = mean(precipitation, na.rm = TRUE))
+  mean_seas <- data_nz %>%
+    group_by(season) %>%
+    summarise(mean_val = mean(value, na.rm = TRUE), .groups = "drop")
   
-  # Compute seasonal mean for nonzero precipitation
-  tot_mon <- data_nz %>% 
-    group_by(month) %>% 
-    summarise(tot_val = sum(precipitation, na.rm = TRUE))
+  tot_mon1 <- data_nz %>%
+    group_by(month, year) %>%
+    summarise(tot_val1 = if (variable == "precipitation") {
+      sum(value, na.rm = TRUE)
+    } else {
+      mean(value, na.rm = TRUE)
+    }, .groups = "drop")
+  
+  tot_mon <- tot_mon1 %>%
+    group_by(month) %>%
+    summarise(tot_val = mean(tot_val1, na.rm = TRUE), .groups = "drop")
   
-  # Generate plot
   plot1 <- ggplot() + 
-    geom_violin(data = data_nz, aes(x = season, y = precipitation, fill = season), linewidth = 0.2) + 
+    geom_violin(data = data_nz, aes(x = season, y = value, fill = season), linewidth = 0.2) + 
     geom_point(data = mean_seas, aes(x = season, y = mean_val), shape = 4) +
     geom_text(data = mean_seas, aes(x = season, y = mean_val, 
-                                    label = round(mean_val, 2)), vjust = -0.8, size = 2.2) + 
-    geom_label(data = p0_seas, aes(x = season, y = max(data_nz$precipitation), 
-                                   label = paste0("P0 = ", round(p0, 2))), size = 2.2) + 
-    scale_fill_manual(values = c("Summer" = "#3da83d",
-                                 "Spring" = "#FFC3A0",
-                                 "Fall" = "#9B2335",
-                                 "Winter" = "#B0E0E6")) +
-    ylab("Seasonal average nonzero precipitation") +
+                                    label = round(mean_val, 2)), vjust = -0.8, size = 2.2) +
+    (if (variable == "precipitation") {
+      list(geom_label(data = p0_seas, aes(x = season, y = max(data_nz$value), 
+                                          label = paste0("P0 = ", round(p0, 2))), size = 2.2))
+    } else {
+      list()
+    }) +
+    scale_fill_manual(values = c("Summer" = "#3da83d", "Spring" = "#FFC3A0",
+                                 "Fall" = "#9B2335", "Winter" = "#B0E0E6")) +
+    ylab(if (variable == "precipitation") {
+      "Seasonal average precipitation"
+    } else {
+      "Seasonal average temperature"
+    }) +
     xlab("Season") +
     theme(
       legend.text = element_text(size = 7),
@@ -98,12 +125,16 @@ cv_plot_season <- function(data) {
       text = element_text(color = "gray25"),
       legend.box.margin = margin(0.5, 0.5, 0.5, 0.5)
     )
-  plot1
+  
   
   plot2 <- ggplot() + 
     geom_col(data = tot_mon, aes(x = as.factor(month), y = tot_val),
              linewidth = 0.2, fill = "skyblue") + 
-    ylab("Monthly total nonzero precipitation") +
+    ylab(if (variable == "precipitation") {
+      "Monthly total nonzero precipitation"
+    } else {
+      "Monthly average temperature"
+    }) +
     xlab("Month") +
     theme(
       legend.text = element_text(size = 7),
@@ -133,10 +164,6 @@ cv_plot_season <- function(data) {
       text = element_text(color = "gray25"),
       legend.box.margin = margin(0.5, 0.5, 0.5, 0.5)
     )
-  plot2
   
-  gg <- ggarrange(plot1, plot2, nrow = 2, ncol = 1)
-  return(gg)
+  ggarrange(plot1, plot2, nrow = 2)
 }
-
-
diff --git a/R/cv_precip_map.R b/R/cv_precip_map.R
new file mode 100644
index 0000000..413cb10
--- /dev/null
+++ b/R/cv_precip_map.R
@@ -0,0 +1,164 @@
+#' Plot Spatial Precipitation Distribution
+#'
+#' This function takes a NetCDF file containing daily precipitation time series, 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 precipitation during the whole period given in the NC file.
+#'                  "monthly" provides the average total precipitation within every month during the whole duration.
+#'                  "annual" provides the annual maxima of daily precipitation during the whole duration.
+#' @return A ggplot object showing the specified precipitation distribution over space.Note: We are using the default "Virdis"
+#'          color scale which could 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
+#' plot_precip_map("Pincher_ck.nc", stat_type = "monthly")
+#' @import ggplot2
+#' @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
+plot_precip_map <- function(nc_file, stat_type = "daily") {
+  
+  stat_type <- tolower(stat_type[1])
+  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_nz <- sf_obj %>%
+      rowwise() %>%
+      mutate(non_zero_mean = mean(c_across(all_of(cols)) != 0, na.rm = TRUE)) %>%
+      ungroup() %>%
+      select(non_zero_mean, geometry)
+    
+    stars_daily_nz <- st_rasterize(daily_nz["non_zero_mean"])
+    
+    p <- ggplot() +
+      geom_stars(data = stars_daily_nz) +
+      coord_equal() +
+      scale_fill_viridis_c(
+        "Nonzero Average Daily Precipitation (mm)",
+        direction = -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_totals <- sf_long %>%
+        group_by(id, year, month) %>%
+        summarise(monthly_total = sum(value, na.rm = TRUE), .groups = "drop")
+      
+      monthly_avg <- monthly_totals %>%
+        group_by(id, month) %>%
+        summarise(avg_total = mean(monthly_total, na.rm = TRUE), .groups = "drop")
+      
+      sf_avg_monthly <- sf_obj %>%
+        mutate(id = row_number()) %>%
+        select(id, geometry) %>%
+        right_join(monthly_avg, 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))
+      )
+      
+      p <- ggplot() +
+        geom_stars(data = stars_monthly) +
+        coord_equal() +
+        facet_wrap(~month) +
+        scale_fill_viridis_c(
+          "Average Monthly Precipitation (mm)",
+          direction = -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_precip = max(value, na.rm = TRUE), .groups = "drop")
+      
+      avg_annual_max <- annual_max %>%
+        group_by(id) %>%
+        summarise(annual_max_avg = mean(max_precip, 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["annual_max_avg"])
+      
+      p <- ggplot() +
+        geom_stars(data = stars_max) +
+        coord_equal() +
+        scale_fill_viridis_c(
+          "Annual Maxima of daily precipitation (mm)",
+          direction = -1,
+          na.value = "transparent"
+        ) +
+        theme_minimal() + xlab("Longitude") + ylab("Latitude") + mytheme +
+        guides(fill = guide_colourbar(title.position = "top", title.hjust = 0.5))
+      
+      print(p)
+    }
+  }
+}
diff --git a/R/cv_temp_map.R b/R/cv_temp_map.R
new file mode 100644
index 0000000..699c54b
--- /dev/null
+++ b/R/cv_temp_map.R
@@ -0,0 +1,162 @@
+#' Plot Spatial Temperature Distribution
+#'
+#' This function takes a NetCDF file containing daily Temperature time series, 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 temperature during the whole period given in the NC file.
+#'                  "monthly" provides the average Temperature within every month during the whole duration.
+#'                  "annual" provides the annual maxima of daily average Temperature during the whole duration 
+#'                  (not commonly used in applications).
+#' @return A ggplot object showing the specified Temperature distribution over space.Note: We are using the default "plasma"
+#'          color scale which could 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
+#' plot_precip_map("Pincher_ck.nc", stat_type = "monthly")
+#' @import ggplot2
+#' @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
+plot_temp_map <- function(nc_file, stat_type = c("daily", "monthly", "annual")) {
+  stat_type <- tolower(stat_type[1])
+  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_nz <- sf_obj %>%
+      rowwise() %>%
+      mutate(temp_mean = mean(c_across(all_of(cols)), na.rm = TRUE)) %>%
+      ungroup() %>%
+      select(temp_mean, geometry)
+    
+    stars_daily_nz <- st_rasterize(daily_nz["temp_mean"])
+    
+    p <- ggplot() +
+      geom_stars(data = stars_daily_nz) +
+      coord_equal() +
+      scale_fill_viridis_c(
+        "Nonzero Average Daily Temperature (°C)",
+        option = "plasma",
+        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_avg <- 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_avg, 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))
+      )
+      
+      p <- ggplot() +
+        geom_stars(data = stars_monthly) +
+        coord_equal() +
+        facet_wrap(~month) +
+        scale_fill_viridis_c(
+          "Average Monthly Temperature (°C)",
+          option = "plasma",
+          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_precip = max(value, na.rm = TRUE), .groups = "drop")
+      
+      avg_annual_max <- annual_max %>%
+        group_by(id) %>%
+        summarise(annual_max_avg = mean(max_precip, 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["annual_max_avg"])
+      
+      p <- ggplot() +
+        geom_stars(data = stars_max) +
+        coord_equal() +
+        scale_fill_viridis_c(
+          "Annual Maxima of daily Temperature (°C)",
+          option = "plasma",
+          na.value = "transparent"
+        ) +
+        theme_minimal() + xlab("Longitude") + ylab("Latitude") + mytheme +
+        guides(fill = guide_colourbar(title.position = "top", title.hjust = 0.5))
+      
+      print(p)
+    }
+  }
+}
+
diff --git a/R/cv_zone_area_raster.R b/R/cv_zone_area_raster.R
index c498c34..eb1e750 100644
--- a/R/cv_zone_area_raster.R
+++ b/R/cv_zone_area_raster.R
@@ -9,7 +9,7 @@
 #' @returns Returns a \pkg{terra} \code{SpatRaster} object (415 rows x 883 columns x 2 layers) of all Canadian grid locations.
 #' @export
 #' @importFrom terra rast
-#' @author Heba Abdelmoaty Kevin Shook
+#' @author Heba Abdelmoaty Kevin SHook
 
 #' @examples
 #' zone_area_grid <- cv_zone_area_raster()
diff --git a/inst/maps_AM.jpg b/inst/maps_AM.jpg
new file mode 100644
index 0000000..0b5cacb
Binary files /dev/null and b/inst/maps_AM.jpg differ
diff --git a/inst/maps_daily.jpg b/inst/maps_daily.jpg
new file mode 100644
index 0000000..f475a8c
Binary files /dev/null and b/inst/maps_daily.jpg differ
diff --git a/inst/maps_monthly.jpg b/inst/maps_monthly.jpg
new file mode 100644
index 0000000..4e16c83
Binary files /dev/null and b/inst/maps_monthly.jpg differ
diff --git a/inst/plot_seasALL.jpg b/inst/plot_seasALL.jpg
deleted file mode 100644
index d454907..0000000
Binary files a/inst/plot_seasALL.jpg and /dev/null differ
diff --git a/man/cv_basin_daily_precip.Rd b/man/cv_basin_daily_precip.Rd
index 78484ad..0bfaa19 100644
--- a/man/cv_basin_daily_precip.Rd
+++ b/man/cv_basin_daily_precip.Rd
@@ -26,7 +26,7 @@ memory (which is what occurs if \code{temp_file = TRUE}), \emph{but} allows the
 to work with very large basins, which may require more memory than is available.}
 }
 \value{
-Returns a data frame with 2 columns: \code{date} and \code{precipitation}. The
+Reruns a data frame with 2 columns: \code{date} and \code{precipitation}. The
 \code{date} is a standard \R date over the interval 2015-01-01 to 2100-12-31, and the
 \code{precipitation} is the basin mean value.
 }
diff --git a/man/cv_plot_TS.Rd b/man/cv_plot_TS.Rd
index 01846b9..a92a872 100644
--- a/man/cv_plot_TS.Rd
+++ b/man/cv_plot_TS.Rd
@@ -2,27 +2,24 @@
 % Please edit documentation in R/cv_plot_TS.R
 \name{cv_plot_TS}
 \alias{cv_plot_TS}
-\title{Plot precipitation time series}
+\title{Plot Precipitation or Temperature Time Series}
 \usage{
-cv_plot_TS(data)
+cv_plot_TS(data, variable = "precipitation")
 }
 \arguments{
-\item{data}{A data frame with two columns: \code{date} (date or datetime) and
-\code{precipitation} (numeric) as returned by \code{cv_basin_daily_precip()}.}
+\item{data}{A data frame with two columns: "date" (date or datetime) and
+"precipitation" (numeric).}
 }
 \value{
-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().
+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().
 }
 \description{
-Creates a time series plot from a data frame of precipitation data
-as returned by \code{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.
 }
 \examples{
 cv_plot_TS(eg_TS)
 }
-\seealso{
-\code{\link{cv_basin_daily_precip}} \code{\link{cv_plot_prob}} \code{\link{cv_plot_season}}
-}
diff --git a/man/cv_plot_map.Rd b/man/cv_plot_map.Rd
new file mode 100644
index 0000000..47f0b9f
--- /dev/null
+++ b/man/cv_plot_map.Rd
@@ -0,0 +1,33 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/cv_plot_map.R
+\name{cv_plot_map}
+\alias{cv_plot_map}
+\title{Plot Spatial Precipitation or Temperature Distribution}
+\usage{
+cv_plot_map(nc_file, stat_type = "daily", variable = "precipitation")
+}
+\arguments{
+\item{nc_file}{A path to a NetCDF (.nc) file readable by the \code{stars} package.}
+
+\item{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.}
+
+\item{variable}{A string starting with 'p' (precipitation) or 't' (temperature).}
+}
+\value{
+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().
+}
+\description{
+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.
+}
+\examples{
+cv_plot_map("Pincher_ck.nc", stat_type = "monthly", variable = "temperature")
+cv_plot_map("Pincher_ck.nc", stat_type = "annual", variable = "precip")
+}
diff --git a/man/cv_plot_prob.Rd b/man/cv_plot_prob.Rd
index 640206e..a1ed3f3 100644
--- a/man/cv_plot_prob.Rd
+++ b/man/cv_plot_prob.Rd
@@ -2,27 +2,24 @@
 % Please edit documentation in R/cv_plot_prob.R
 \name{cv_plot_prob}
 \alias{cv_plot_prob}
-\title{Plot probability of exceedance for precipitation}
+\title{Plot Probability of Exceedance for Precipitation}
 \usage{
-cv_plot_prob(data)
+cv_plot_prob(data, variable = "precipitation")
 }
 \arguments{
-\item{data}{A data frame with two columns: \code{date} (date or datetime) and
-\code{precipitation} (numeric) as returned by \code{cv_basin_daily_precip()}.}
+\item{data}{A data frame with two columns: "date" (date or datetime) and
+"precipitation" (numeric).}
 }
 \value{
-A ggplot object displaying the probability of exceedance of nonzero
-precipitation. The returned plots look best when saved at the size
-16.5 x 12 cm. You can easily change the font sizes using theme().
+A ggplot object displaying the probability of exceedance of nonzero precipitation.
+The returned plots look best when saved at the size 16.5 x 12 cm.
+You can easily change the font sizes using theme().
 }
 \description{
-Plots the exceedance probabilities of non-zero values in a precipitation time
-series as returned by \code{cv_basin_daily_precip()}, vs the precipitation value.
-Uses a logarithmic scale for the y-axis.
+This function takes a data frame containing a precipitation time series,
+filters out zero precipitation values, calculates the probability of exceedance,
+and generates a scatter plot with a logarithmic y-axis.
 }
 \examples{
 cv_plot_prob(eg_TS)
 }
-\seealso{
-\code{\link{cv_basin_daily_precip}} \code{\link{cv_plot_TS}} \code{\link{cv_plot_season}}
-}
diff --git a/man/cv_plot_season.Rd b/man/cv_plot_season.Rd
index 204e20f..d60014a 100644
--- a/man/cv_plot_season.Rd
+++ b/man/cv_plot_season.Rd
@@ -2,27 +2,23 @@
 % Please edit documentation in R/cv_plot_season.R
 \name{cv_plot_season}
 \alias{cv_plot_season}
-\title{Plot seasonal precipitation distribution}
+\title{Plot Seasonal Distribution of a Hydroclimatic Variable}
 \usage{
-cv_plot_season(data)
+cv_plot_season(data, variable = "precipitation")
 }
 \arguments{
-\item{data}{A data frame with two columns: \code{date} (date or datetime) and
-\code{precipitation} (numeric) as returned by \code{cv_basin_daily_precip()}.}
+\item{data}{A data frame with two columns: "date" (date or datetime) and a variable column (e.g., "precipitation" or "temperature").}
+
+\item{variable}{A character string, either "precipitation" or "temperature". Determines plotting behavior and calculation of P0.}
 }
 \value{
-A ggplot object displaying the seasonal distribution of nonzero precipitation,
-with mean values and P0 labelled. The returned plots look best when saved at
-the size 16.5 x 14 cm. You can easily change the font sizes using theme().
+A ggplot object showing seasonal distribution with summary values and optional P0.
 }
 \description{
-Creates seasonal plots of a precipitation time series, as returned by \code{cv_basin_daily_precip()}. Assigns seasons based on the month,
-computes the proportion of zero precipitation values (P0) and seasonal means for nonzero precipitation,
-and generates a violin plot showing the seasonal distribution of precipitation.
+This function takes a data frame with a hydroclimatic time series (precipitation or temperature), assigns seasons based on months,
+computes the proportion of zero precipitation values (P0) and seasonal means for nonzero values,
+and generates a violin plot showing the seasonal distribution.
 }
 \examples{
-cv_plot_season(eg_TS)
-}
-\seealso{
-\code{\link{cv_basin_daily_precip}} \code{\link{cv_plot_TS}} \code{\link{cv_plot_prob}}
+cv_plot_season(data, variable = "precipitation")
 }
diff --git a/man/cv_zone_area_raster.Rd b/man/cv_zone_area_raster.Rd
index a7161b7..59d9abe 100644
--- a/man/cv_zone_area_raster.Rd
+++ b/man/cv_zone_area_raster.Rd
@@ -22,5 +22,5 @@ plot(zone_area_grid)
 
 }
 \author{
-Heba Abdelmoaty Kevin Shook
+Heba Abdelmoaty Kevin SHook
 }
diff --git a/man/plot_precip_map.Rd b/man/plot_precip_map.Rd
new file mode 100644
index 0000000..13b2089
--- /dev/null
+++ b/man/plot_precip_map.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/cv_precip_map.R
+\name{plot_precip_map}
+\alias{plot_precip_map}
+\title{Plot Spatial Precipitation Distribution}
+\usage{
+plot_precip_map(nc_file, stat_type = "daily")
+}
+\arguments{
+\item{nc_file}{A path to a NetCDF (.nc) file readable by the \code{stars} package.}
+
+\item{stat_type}{One of "daily", "monthly", or "annual" to indicate which type of statistic to visualize.
+"daily" provides the nonzero average of daily precipitation during the whole period given in the NC file.
+"monthly" provides the average total precipitation within every month during the whole duration.
+"annual" provides the annual maxima of daily precipitation during the whole duration.}
+}
+\value{
+A ggplot object showing the specified precipitation distribution over space.Note: We are using the default "Virdis"
+color scale which could 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().
+}
+\description{
+This function takes a NetCDF file containing daily precipitation time series, extracts spatial and temporal features,
+and plots one of three options: daily nonzero frequency, average monthly totals, or average annual maxima.
+}
+\examples{
+plot_precip_map("Pincher_ck.nc", stat_type = "monthly")
+}
diff --git a/man/plot_temp_map.Rd b/man/plot_temp_map.Rd
new file mode 100644
index 0000000..a9d3756
--- /dev/null
+++ b/man/plot_temp_map.Rd
@@ -0,0 +1,30 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/cv_temp_map.R
+\name{plot_temp_map}
+\alias{plot_temp_map}
+\title{Plot Spatial Temperature Distribution}
+\usage{
+plot_temp_map(nc_file, stat_type = c("daily", "monthly", "annual"))
+}
+\arguments{
+\item{nc_file}{A path to a NetCDF (.nc) file readable by the \code{stars} package.}
+
+\item{stat_type}{One of "daily", "monthly", or "annual" to indicate which type of statistic to visualize.
+"daily" provides the nonzero average of daily temperature during the whole period given in the NC file.
+"monthly" provides the average Temperature within every month during the whole duration.
+"annual" provides the annual maxima of daily average Temperature during the whole duration
+(not commonly used in applications).}
+}
+\value{
+A ggplot object showing the specified Temperature distribution over space.Note: We are using the default "plasma"
+color scale which could 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().
+}
+\description{
+This function takes a NetCDF file containing daily Temperature time series, extracts spatial and temporal features,
+and plots one of three options: daily nonzero frequency, average monthly totals, or average annual maxima.
+}
+\examples{
+plot_precip_map("Pincher_ck.nc", stat_type = "monthly")
+}
diff --git a/vignettes/vignette_CMIP6VisR.R b/vignettes/vignette_CMIP6VisR.R
deleted file mode 100644
index 05e9ca7..0000000
--- a/vignettes/vignette_CMIP6VisR.R
+++ /dev/null
@@ -1,63 +0,0 @@
-## ----setup, include=FALSE-----------------------------------------------------
-knitr::opts_chunk$set(eval = TRUE,
-                      echo = TRUE,
-                      fig.width = 7, 
-                      warning = FALSE,
-                      message = FALSE)
-
-## ----warning=FALSE, message=FALSE---------------------------------------------
-library(CMIP6VisR)
-library(terra)
-
-## -----------------------------------------------------------------------------
-# Read the shapefile containing the basin
-fpath <- system.file("extdata", "07BF001.shp", package = "CMIP6VisR")
-basin_vector <- vect(fpath)
-plot(basin_vector)
-
-## ----warning=FALSE------------------------------------------------------------
-# convert the dataframe of zones to raster (built-in the package)
-zone_area_grid <- cv_zone_area_raster() 
-
-req_data <- cv_clip_basin(zone_area_grid, basin_vector)
-
-# printing the zones that are required for the targeted area
-print(req_data[["zone"]])
-
-## ----message=FALSE, warning=FALSE---------------------------------------------
-# Isolating the raster outputs
-raster_list <- req_data$raster
-
-terra::plot(raster_list[[1]], main = expression("Zone 4, Grid Areas in " ~ km^2))
-
-terra::plot(raster_list[[2]], main = expression("Zone 7, Grid Areas in " ~ km^2))
-
-# Merge the individual rasters using terra::merge
-merged_raster <- do.call(terra::merge, req_data$raster)
-
-# Plot the merged raster
-terra::plot(merged_raster, main = expression("Both zones, Grid Areas in " ~ km^2))
-
-
-## ----eval=FALSE---------------------------------------------------------------
-# values <- cv_basin_daily_precip(netcdf_directory = "./data/",
-#                                 scenario = "pr_day_ACCESS-CM2_ssp126_r2i1p1f1_gn_20150101-21001231_cannc_SPQM_",
-#                                 basin_zone_area = req_data,
-#                                 temp_file = FALSE)
-
-## ----include=FALSE------------------------------------------------------------
-data("eg_TS")
-values = eg_TS
-
-## -----------------------------------------------------------------------------
-head(values)
-
-## -----------------------------------------------------------------------------
-cv_plot_TS(values)
-
-## -----------------------------------------------------------------------------
-cv_plot_season(values)
-
-## -----------------------------------------------------------------------------
-cv_plot_prob(values)
-
diff --git a/vignettes/vignette_CMIP6VisR.Rmd b/vignettes/vignette_CMIP6VisR.Rmd
index aac9212..4ae83fa 100644
--- a/vignettes/vignette_CMIP6VisR.Rmd
+++ b/vignettes/vignette_CMIP6VisR.Rmd
@@ -129,7 +129,8 @@ The function `cv_basin_daily_precip` has several parameters:
 2. *scenario* the scenario argument which indicates the name of the simulation (the part of the string that is before the zone number). For example: "pr_day_ACCESS-CM2_ssp126_r2i1p1f1_gn_20150101-21001231_cannc_SPQM_"
 3. *basin_zone_area*: the output resulting from `cv_clip_basin` getting the grid area information of the target basin. 
 4. *temp_file*: if `FALSE` (the default), the function tries to complete all operations in memory, 
-which is relatively fast. However, if you are using a very large basin, then you might run of memory. In that case, you would need to specify `temp_file = TRUE`.
+which is relatively fast. However, if you are using a very large basin, then you might run of memory.
+In this case, you need to specify `temp_file = TRUE`.
 
 The function uses the areas for each cell to weight the precipitation from the NetCDF file(s).
 
@@ -146,7 +147,7 @@ data("eg_TS")
 values = eg_TS
 ```
 
-This data should be in a data frame containing two columns, named "date" and "precipitation".
+This data should be in dataframe format that contains two columns, named "date" and "precipitation"
 
 ```{r}
 head(values)
@@ -155,34 +156,64 @@ head(values)
 
 ## Visualization functions
 
-The visualization functions are used to generate precipitation time series, seasonal and monthly precipitation distributions, and probability exceedance plots. All functions returns a `ggplot2` object which can be altered.
+we introduce three visualization functions used to generate precipitation time series, seasonal and monthly precipitation distributions, and probability exceedance plots. 
+
 
 
 1. Precipitation Time Series
-For analyzing precipitation trends over long periods, a time series plot is used. This visualization enables an understanding of interannual variability and extreme precipitation events over time. The mean, standard deviation and p0 (probability of zero values) are
-listed in the title. 
+For analyzing precipitation trends over long periods, a time series plot is used. This visualization enables an understanding of interannual variability and extreme precipitation events over time.
 
 ```{r}
-cv_plot_TS(values)
+cv_plot_TS(values, variable = "precipitation")
 ```
 
 2. Seasonal and Monthly Precipitation Distribution
-The function `cv_plot_season` is useful to analyze precipitation variations across seasons and months. The functio returns a plot with two sub-plots:
-
-1. violin plots of the seasonal mean nonzero precipitation, with the coorresponding p0 value displayed, and
-2. a bar graph of the total precipitation by month.
+To analyze precipitation variations across seasons and months, the following function is used. These plots highlight seasonal trends and monthly variations in precipitation data.
 
 ```{r}
-cv_plot_season(values)
+cv_plot_season(values, variable = "precipitation")
 ```
 
 
 3. Probability of Exceedance Plot
-The probability of exceedance plot is useful for assessing extreme precipitation events. The precipitation data are sorted in descending order and the the exceedance probabilities
-are calculated from the rankings. The log-scale for probability ensures better visualization of extreme events.
+The probability of exceedance plot is a critical tool for assessing extreme precipitation events. This plot is generated using the following function. This function sorts the precipitation data in descending order and calculates exceedance probabilities. The log-scale for probability ensures better visualization of extreme events.
 
 ```{r}
-cv_plot_prob(values)
+
+cv_plot_prob(values, variable = "precipitation")
+```
+
+4. Spatial distribution
+
+4.1. Precipitation
+The following plot shows a map of the average nonzero daily, average total monthly, or annual maxima of daily precipitation.
+```{r eval=FALSE}
+cv_plot_map("ncpath/ncfile.nc", stat_type = "daily")
+```
+
+```{r eval=FALSE}
+cv_plot_map("ncpath/ncfile.nc", stat_type = "monthly")
+```
+
+```{r eval=FALSE}
+cv_plot_map("ncpath/ncfile.nc", stat_type = "annual")
 ```
 
 
+4.2. Temperature
+The following plot shows a map of the average daily, average  monthly, or annual maxima of daily temperature
+
+```{r eval=FALSE}
+cv_plot_map("ncpath/ncfile.nc", stat_type = "daily")
+```
+
+
+```{r eval=FALSE}
+cv_plot_map("ncpath/ncfile.nc", stat_type = "monthly")
+```
+
+
+```{r eval=FALSE}
+cv_plot_map("ncpath/ncfile.nc", stat_type = "annual")
+```
+