Include temporal aggregtion per month prior to aggregation per season

Author: ma-z-am

R/seasonal_index_per_pixel.R

@@ -14,24 +14,30 @@
 #'   this mode should be checked before use.
 #' @param pixelQualityThreshold Numeric between 0 and 1. Minimum water-scene
 #'   proportion required for a pixel to be retained in the seasonal index.
+#' @param per_month_first Logical, If `TRUE`, the index is first averaged on a 
+#'   monthly basis and than averaged for the whole season.
 #'
 #' @return A list with `x`, `y`, `crs`, `year`, `QualityThreshold`, `RSindex`, and
 #'   `sceneProportions`. `RSindex` is a matrix of seasonal pixel index values;
 #'   pixels below the water-quality threshold are set to `NA`.
 #'
 #' @details Water-scene proportions are calculated by [waterscene_proportion()].
 #'   A warning is issued if no pixel meets `pixelQualityThreshold`.
+#'   
+#'   The idea of per_month_first is to create a balanced mean throughout the
+#'   season, independently of cloudy or sunny periods over a long time.
 #'
 #' @export
 seasonal_index_per_pixel <- function(
-    imageIndex, nc, water_scenes_only = TRUE, pixelQualityThreshold = 0.8
+    imageIndex, nc, water_scenes_only = TRUE, pixelQualityThreshold = 0.8,
+    per_month_first = TRUE
 ){
   output <- list()
   #output[["lakeInfo"]] <- lakeInfo
-  imageMonths <- unique(as.numeric(format(imageIndex$t_date, "%m")))
+  imageMonthsVector <- format(imageIndex$t_date, "%m")
+  imageMonths <- unique(as.numeric(imageMonthsVector))
   imageYear <- unique(as.numeric(format(imageIndex$t_date, "%Y")))
 
-  cat("Loading values of SCL Band ... \n")
   sclt <- load_BandLayer(
     nc = nc, 
     band = "SCL", 
@@ -65,10 +71,19 @@ seasonal_index_per_pixel <- function(
       )
     }
 
-  RSindex <- pixel_wise_average(
-    list_of_mats = imageIndex$RSindex, 
-    na.rm = TRUE
-  )
+  if(per_month_first){
+    monthlySplitIndex <- split(x = imageIndex$RSindex, f = imageMonthsVector)
+    RSindexMonth <- lapply(monthlySplitIndex, pixel_wise_average, na.rm = TRUE)
+    RSindex <- pixel_wise_average(
+      list_of_mats = RSindexMonth, 
+      na.rm = TRUE
+    )
+  } else {
+    RSindex <- pixel_wise_average(
+      list_of_mats = imageIndex$RSindex, 
+      na.rm = TRUE
+    )
+  }
 
   waterLayer <- lakeRS::waterscene_proportion(scl_image = sclt$band)