Conflicts solved (the best from both changes)

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

# Conflicts:
#	R/dynamic_per_pixel.R

Author: ma-z-am

R/dynamic_per_pixel.R

@@ -27,7 +27,6 @@
 #' 
 #' @export
 #' 
-#' 
 dynamic_per_pixel <- function(
     ncImage, t_date, years, water_scenes_only = TRUE, days_around_ma = 20, 
     maxPixels = 1000, threshold = NULL, lakeInfo = c("", ""), maxDataPoints = 5000000
@@ -42,22 +41,31 @@ dynamic_per_pixel <- function(
 
   indexImages <- ncImage$RSindex[yearFilter]
   imageDOY <- imageDOY[yearFilter]
-  d <- dim(indexImages[[1]])
-  
-  sclImage <- ncImage$SCL[yearFilter]
+  sclImages <- ncImage$SCL[yearFilter]
 
-  cat("Image data is filtered pixel by pixel ... \n")
+  d <- dim(indexImages[[1]])
 
-  if(water_scenes_only){
-    indexImages <- lapply(seq_along(indexImages), function(i){
-      indexImages[[i]][sclImage[[i]] != 6] <- NA
-      indexImages[[i]]
-    })
-  } else {
-    indexImages <- lapply(seq_along(indexImages), function(i){
-      indexImages[[i]][sclImage[[i]] %in% c(8:11)] <- NA
-      indexImages[[i]]
-    })     
+  # rm(ncImage)
+  # gc()
+  # if(object.size(ncImage)/1E+09 > 10){
+  #   
+  # }
+  
+  cat("Image data is filtered for SCL categories pixel by pixel ... \n")
+  
+  for(i in seq_along(indexImages)){
+    indexImages[[i]] <- if(water_scenes_only){
+      scl_filter(
+        indexImage = indexImages[[i]], 
+        sclImage = sclImages[[i]], 
+        bands = 6, 
+        invert = TRUE)
+    } else {
+      scl_filter(
+        indexImage = indexImages[[i]], 
+        sclImage = sclImages[[i]], 
+        bands = c(8:11))    
+    }
   }
 
   t_doy <- imageDOY[order(imageDOY)]
@@ -81,6 +89,19 @@ dynamic_per_pixel <- function(
 
   valid_values <- apply(ts, 1 , function(p_ts){sum(!is.na(p_ts))})
   pixel_selection <- valid_values > threshold
+  
+  if(sum(pixel_selection) == 0){
+    return(
+      list("moving_averages" = NULL,
+           "raster_location" = data.frame(
+             "pixel" = NA,
+             "i_col" = 0,
+             "i_row" = 0,
+             "valid_values" = 0)
+      )
+    )
+  }
+  
   pixel_selection_i <- which(pixel_selection)
   nAvailable <- length(pixel_selection_i)
   if(nAvailable > maxPixels){
@@ -119,7 +140,7 @@ dynamic_per_pixel <- function(
     t_doy = t_doy,
     days_around_ma = days_around_ma
   )
-
+  
   cat(paste0("Processing ", nm+1, " matrices ... \n"))
   df_out <- lapply(ts_start, function(ts){
     cat(" | matrix done")
@@ -156,7 +177,6 @@ dynamic_per_pixel <- function(
   )
 }
 
-
 #' Subset of all images to be used for moving averages for each day
 #' 
 #' @param t_doy Numeric vector of the day of the year according to ts_pixel 
@@ -166,7 +186,7 @@ dynamic_per_pixel <- function(
 #' 
 images_per_ma <- function(t_doy, days_around_ma){
   lapply(1:365, function(doy){
-    d_range <- doy + c(-20, 20)
+    d_range <- doy + c(-days_around_ma, days_around_ma)
     d_range[d_range < 1] <- d_range[d_range < 1] + 365
     d_range[d_range > 365] <- d_range[d_range > 365] - 365
 
@@ -178,31 +198,18 @@ images_per_ma <- function(t_doy, days_around_ma){
   })
 }
 
-#' Moving average for a sorted pixel vector with according day of the year time 
-#' vector
+#' Set values NA by SCL
 #' 
-#' @param ts_pixel Data for one pixel of multiple images (ts)
-#' @param t_doy Numeric vector of the day of the year according to ts_pixel 
-#' (same length and order)
-#' @param days_around_ma The days around the moving average day (i.e. 10  means
-#' 10 days before and ten days after the actual day are used for averaging)
+#' @param indexImage,sclImage Index and SCL layer of the image
+#' @param bands Numeric vector specifiyng the SCL categories to filtered for
+#' @param invert If TRUE the filtering is inverted -> specefied bands are removed
 #' 
-#' 
-# moving_average <- function(ts_pixel, t_doy, days_around_ma){
-#   df_out <- data.frame("doy" = 1:365)
-#   
-#   df_out$ma <- sapply(df_out$doy, function(doy){
-#     d_range <- doy + c(-days_around_ma, days_around_ma)
-#     d_range[d_range < 1] <- d_range[d_range < 1] + 365
-#     d_range[d_range > 365] <- d_range[d_range > 365] - 365
-#     
-#     days_for_ma <- if(d_range[2] > d_range[1]){
-#       t_doy > d_range[1] & t_doy < d_range[2]
-#     } else {
-#       t_doy > d_range[1] | t_doy < d_range[2]
-#     }
-#     mean(ts_pixel[days_for_ma], na.rm = TRUE) 
-#   })
-#   df_out
-# }
-
+scl_filter <- function(indexImage, sclImage, bands, invert = FALSE){
+  if(invert){
+    indexImage[!(sclImage %in% bands)] <- NA
+  } else {
+    indexImage[!(sclImage %in% bands)] <- NA
+  }
+  indexImage
+}
+>>>>>>> c91d979d3f068e2d74c8eac59c874cbb27f57666

R/load_netcdf.R

@@ -25,19 +25,21 @@ load_netcdf <- function(
     start_xy = c(10,10), count_xy = c(NA,NA)
 ){
   ncMeta <- ncdump::NetCDF(x = file.path(filePath, fileName))
+  dimDf <- ncMeta$dimension
+  xLen <- dimDf$len[dimDf$name == "x"]
+  yLen <- dimDf$len[dimDf$name == "y"]
   if(any(is.na(count_xy))){
     all_values <- c("x", "y")[which(is.na(count_xy))]
     if("x" %in% all_values){
-      count_xy[1] <- ncMeta$dimension$len[ncMeta$dimension$name == "x"]
+      count_xy[1] <- xLen
     }
     if("y" %in% all_values){
-      count_xy[2] <- ncMeta$dimension$len[ncMeta$dimension$name == "y"]
+      count_xy[2] <-  yLen
     }
   }
 
-  dimDf <- ncMeta$dimension
-  to_many_counts <- c(dimDf$len[dimDf$name == "x"], dimDf$len[dimDf$name == "y"]) - 
-    (count_xy + start_xy - 1)
+ 
+  to_many_counts <- c(xLen, yLen) - (count_xy + start_xy - 1)
   if(any(to_many_counts < 0)){
     count_xy[to_many_counts < 0] <- count_xy[to_many_counts < 0] + to_many_counts[to_many_counts < 0] 
   }
@@ -56,8 +58,7 @@ load_netcdf <- function(
     vars <- availableVars
   }
 
-  if(count_xy[1] < dimDf$len[dimDf$name == "x"] |
-     count_xy[2] < dimDf$len[dimDf$name == "y"]){
+  if(count_xy[1] < xLen |count_xy[2] < yLen){
     crop_nc <- TRUE
   }
 
@@ -66,10 +67,12 @@ load_netcdf <- function(
   t_id <-  dimDf$id[dimDf$name == "t"]
 
   valDf <- ncMeta$dimension_values
+  initial_x <- valDf$vals[valDf$id == x_id]
+  initial_y <- valDf$vals[valDf$id == y_id]
   cat(paste0("Loading Coordinates and time variables ... \n"))
   nc <- list()
-  nc[["x"]] <- valDf$vals[valDf$id == x_id]
-  nc[["y"]] <- valDf$vals[valDf$id == y_id]
+  nc[["x"]] <- valDf$vals[valDf$id == x_id][start_xy[1]:(start_xy + count_xy - 1)[1]]
+  nc[["y"]] <- valDf$vals[valDf$id == y_id][start_xy[2]:(start_xy + count_xy - 1)[2]]
   nc[["t"]] <- valDf$vals[valDf$id == t_id]
   nc[["t_date"]] <- as.Date(nc[["t"]], origin = "1990-01-01")
 
@@ -91,8 +94,8 @@ load_netcdf <- function(
         )
       if(crop_nc){
         bb_crop <- c(
-          range(nc$x[v_crop[1:2]]),
-          range(nc$y[v_crop[3:4]])
+          range(initial_x[v_crop[1:2]]),
+          range(initial_y[v_crop[3:4]])
         )
         e <- raster::extent(bb_crop)
         nc[[varName]] <- raster::crop(nc[[varName]], e)