Merge branch 'dev' of https://github.com/KWB-R/lakeRS into dev

Author: JaKovonBER

R/dynamic_per_pixel.R

@@ -15,8 +15,12 @@
 #' @param threshold The mininum number of valid values for one pixel to be 
 #' processed (the higher the days around moving averages are, the lower this
 #' threshold can be)
-#' @param maxPixels Maximum number of pixels (randomly chosen from available pixels).
-#' If Inf, all pixels will be analysed.
+#' @param maxPixels If the number is below 1 it will be interpreted as a 
+#' proportion of all valid Pixels. Otherwise, maxPixels defines the number of 
+#' pixels used to calculate the overall dynamic.If it is smaller as the number of 
+#' valid pixels (revided by the threshold value), pixels will drawn in decreasing 
+#' order from most valid to lowest valid. If Inf, all pixels will be analysed. 
+#' MaxPixels is ignored if pixelFilter is defined.
 #' @param pixelFilter The ID of pixels to be used (The ID of a pixel is its 
 #' its number in the the original matrix of the netcdf).
 #' @param maxDataPoints The number of datapoints to be part of one matrix
@@ -31,8 +35,6 @@
 #' are available. If the number of images is not sufficient the days_around_ma
 #' need to be increased and the results will become smoother.
 #' 
-#' 
-#' 
 #' @export
 #' 
 dynamic_per_pixel <- function(
@@ -55,18 +57,19 @@ dynamic_per_pixel <- function(
 
   cat("Image data is filtered for SCL categories pixel by pixel ... \n")
   for(i in seq_along(indexList)){
-    indexList[[i]] <- if(water_scenes_only){
-      scl_filter(
-        indexImage = indexList[[i]], 
-        sclImage = sclList[[i]], 
-        bands = 6, 
-        invert = TRUE)
-    } else {
-      scl_filter(
-        indexImage = indexList[[i]], 
-        sclImage = sclList[[i]], 
-        bands = c(8:11))    
-    }
+    indexList[[i]] <- 
+      if(water_scenes_only){
+        scl_filter(
+          indexImage = indexList[[i]], 
+          sclImage = sclList[[i]], 
+          bands = 6, 
+          invert = TRUE)
+      } else {
+        scl_filter(
+          indexImage = indexList[[i]], 
+          sclImage = sclList[[i]], 
+          bands = c(8:11))    
+      }
   }
 
   imageDOY <- as.numeric(format(imageIndex$t_date, "%j"))
@@ -78,19 +81,16 @@ dynamic_per_pixel <- function(
   if(!is.null(pixelFilter)){
     pixel_selection_i <- pixelFilter
     nAvailable <- np <- length(pixelFilter)
-    maxPixels <- Inf
     filteredIndexList <- lapply(indexList, function(x){x[pixelFilter]})
     ts_select <- matrix(
       data = unlist(filteredIndexList), 
       nrow = np, 
       ncol = length(filteredIndexList)
     )
   } else {
-    ts <- matrix(
-      data = unlist(indexList), 
-      nrow = d[1] * d[2], 
-      ncol = length(indexList)
-    )
+    ts <- sapply(indexList, function(x){
+      unlist(x)
+    })
     valid_values <- apply(ts, 1 , function(p_ts){sum(!is.na(p_ts))})
     med_values <- median(valid_values[valid_values > 0])
     if(is.null(threshold)){
@@ -99,7 +99,17 @@ dynamic_per_pixel <- function(
     if(is.null(days_around_ma)){
       days_around_ma <- ceiling(365 / med_values * 3)
     }
-    pixel_selection <- valid_values > threshold
+    pixel_selection <- valid_values >= threshold
+    if(maxPixels <= 1){
+      maxPixels <- round(sum(pixel_selection) * maxPixels)
+    }
+    pixel_selection_i <- 
+      if(maxPixels < sum(pixel_selection)){
+        order(valid_values, decreasing = TRUE)[1:maxPixels]
+      } else {
+        which(pixel_selection)
+      }
+   
     if(sum(pixel_selection) == 0){
       return(
         list("moving_averages" = NULL,
@@ -111,12 +121,6 @@ dynamic_per_pixel <- function(
         )
       )
     }
-    pixel_selection_i <- which(pixel_selection)
-    nAvailable <- length(pixel_selection_i)
-    if(nAvailable > maxPixels){
-      pixel_selection_i <- sample(pixel_selection_i, maxPixels)
-      nAvailable <- length(pixel_selection_i)
-    }
     ts_select <- ts[pixel_selection_i,]
     rm(ts)
   }
@@ -146,7 +150,7 @@ dynamic_per_pixel <- function(
   }
 
   cat(paste0(
-    "Calculate moving averages of ", nAvailable, " pixel timesseries ... \n")
+    "Calculate moving averages of ", length(pixel_selection_i), " pixel timesseries ... \n")
   )
   ipa <- images_per_ma(
     t_doy = t_doy,

R/load_BandLayer.R

@@ -15,6 +15,7 @@
 load_BandLayer <- function(
     nc, band, monthFilter = NULL, yearFilter = NULL
 ){
+  cat(paste("Loading values of Band", band, "... \n"))
   if(!is.null(monthFilter)){
     mf <- as.integer(format(nc$t_date, "%m")) %in% as.integer(monthFilter)
     if(sum(mf) > 0){

R/ndi_per_image.R

@@ -39,15 +39,13 @@ ndi_per_image <- function(
   #   nc$SpRast[[IN]] <- 
   #     (nc$SpRast[[AN]] - nc$SpRast[[BN]]) / (nc$SpRast[[AN]] + nc$SpRast[[BN]])
   # }
-  cat(paste("Loading values of Band", bandNames[1], "... \n"))
   At <- load_BandLayer(
     nc = nc, 
     band = bandNames[1], 
     monthFilter = monthFilter, 
     yearFilter = year
   )
 
-  cat(paste("Loading values of Band", bandNames[2], "... \n"))
   Bt <- load_BandLayer(
     nc = nc, 
     band = bandNames[2],

inst/extdata/scripts/example.R

@@ -59,15 +59,35 @@ lakeRS::download_openEO_job(
 filePath <- "C:/Users/mzamzo/Documents/tmp/pcl/01_data/input/ndtri_openeo/Koerbaer Teich_lakeRS-example_2020-01-01_2021-01-01"
 filePath <- "C:/Users/mzamzo/Documents/tmp/pcl/01_data/input/ndtri_openeo/TrichonisLake_DS1_20190101-20200101"
 
+# further input
+lakeName <- 
+  "Lake Trichonis" 
+  # "Körbaer Teich"
+lakeID <- 
+  "DS1" 
+  # "800015388119" 
+indexName <- 
+  "NDTrI"
+  #"NDCI"
+  # "NDTI"
+  # "NDSSI"
+
+indexBands <- 
+  # c("B05", "B02")
+  # c("B05", "B03")
+  # c("B04", "B03")
+  # c("B02", "B08")
+
+#
 # open netcdf connection
 nc <- lakeRS::open_netcdf(filePath = filePath)
 
 # Index per image 
 range(nc$t_date)
 imageIndex <- lakeRS::ndi_per_image(
   nc = nc, 
-  year = 2020, 
-  bandNames = c("B05", "B02")
+  year = 2019, 
+  bandNames = indexBands
 )
 
 # yearly dynamic
@@ -83,18 +103,26 @@ imageIndex <- lakeRS::ndi_per_image(
     v_averageList = list(indexDynamic$lakeDynamic$q_0.5), 
     df_q50List = list(indexDynamic$lakeDynamic[,c("q_0.25", "q_0.75")]), 
     df_q95List = list(indexDynamic$lakeDynamic[,c("q_0.025", "q_0.975")]),
-    lakeName = "Körbaer Teich", 
-    ylab = "NDTrI per timestep")
+    lakeName = lakeName, 
+    ylab = paste0(indexName, " per timestep"))
+  
+  lakeRS::plot_dynamic(
+    v_averageList = list(indexDynamic$lakeDynamic$q_0.5), 
+    df_q50List = list(indexDynamic$lakeDynamic[,c("q_0.25", "q_0.75")]), 
+    df_q95List = NULL,
+    lakeName = lakeName, 
+    ylab = paste0(indexName, " per timestep"))
 
 
   lakeRS::best_nk(pixelDynamic = indexDynamic$pixelDynamics)
+  k <- 4
 
   pClusters <- lakeRS::pixel_clusters(
     pixelDynamic = indexDynamic$pixelDynamics, 
     nc = nc, correlate_first = TRUE,
-    k = 3)
+    k = k)
 
-  DynamicStatList <- lapply(1:3, function(CLUSTER){
+  DynamicStatList <- lapply(1:k, function(CLUSTER){
     lakeRS::dynamic_per_pixel(
       imageIndex = imageIndex, 
       nc = nc, 
@@ -109,13 +137,11 @@ imageIndex <- lakeRS::ndi_per_image(
   lakeRS::plot_dynamic(
     v_averageList = lapply(DynamicStatList, function(x){x$q_0.5}), 
     df_q50List = lapply(DynamicStatList, function(x){data.frame(x$q_0.25, x$q_0.75)}),
-    lakeName = "Körbaer Teich", 
-    TSnames = c("C1", "C2", "C3")
+    lakeName = lakeName, 
+    TSnames = paste0("C", 1:k)
   )
 
-    DynamicStatList = DynamicStatList, 
-    lakeName = "Körbaer Teich", 
-    smallBandOnly = TRUE)
+ 
 }
 
 # Temporal aggregation
@@ -140,7 +166,7 @@ if(FALSE){
 # spatial aggregation
 lakeIndex <- lakeRS::seasonal_index_per_lake(
   seasonIndex = seasonIndex, 
-  lakeName = "Koerbaer_teich",
+  lakeName = lakeName,
   lakeID = "800015388119"
 ) 
 
@@ -151,8 +177,8 @@ yearlyLakes <- lakeRS::combine_years_lakeData(lakeIndexList = list(lakeIndex))
 if(FALSE){
   lakeRS::plot_lake_index_histogram(
     lakeIndex = lakeIndex, 
-    lakeName = "Koerbaer_teich", 
-    indexName = "NDTrI"
+    lakeName = lakeName, 
+    indexName = indexName
   )
 
   # some plot as above
@@ -188,11 +214,23 @@ classData <- lakeRS::discreteClassAssessment(
     classColors = lakeRS::tenClassColors$color,
     plotLegend = FALSE) 
 }
-classData <- lakeRS::numericAssessment(
+
+numData <- lakeRS::numericAssessment(
   yearly_spread = yearlyPixels$indexTable
 )
 
 
+
+
+
+
+
+
+
+
+# -------------------------------------------------------
+# old
+
 xScene <- lakeRS::nc_scene_per_image(
   nc = nc, 
   scene = "water"

inst/extdata/scripts/example_manyYears.R

@@ -5,7 +5,6 @@ system.time(
   nc <- lakeRS::open_netcdf(filePath = filePath)
 )
 
-
 # Dynmaic
 DynamicList <- list()
 years <- 2018:2025