# activeer de renv omgeving; # renv::activate('~/smartCane/r_app') # renv::restore() library(here) library(sf) library(terra) library(tidyverse) library(lubridate) # library(exactextractr) # library(readxl) #funcion CI_prep package # Vang alle command line argumenten op args <- commandArgs(trailingOnly = TRUE) # Controleer of er ten minste één argument is doorgegeven if (length(args) == 0) { stop("Geen argumenten doorgegeven aan het script") } # Converteer het eerste argument naar een numerieke waarde end_date <- as.Date(args[1]) offset <- as.numeric(args[2]) # Controleer of weeks_ago een geldig getal is if (is.na(offset)) { # stop("Het argument is geen geldig getal") offset <- 7 } # Converteer het tweede argument naar een string waarde project_dir <- as.character(args[3]) # Controleer of data_dir een geldige waarde is if (!is.character(project_dir)) { project_dir <- "chemba" } laravel_storage_dir <- here("laravel_app/storage/app", project_dir) #preparing directories reports_dir <- here(laravel_storage_dir, "reports") planet_tif_folder <- here(laravel_storage_dir, "merged_tif") merged_final <- here(laravel_storage_dir, "merged_final_tif") new_project_question = FALSE planet_tif_folder <- here(laravel_storage_dir, "merged_tif") merged_final <- here(laravel_storage_dir, "merged_final_tif") data_dir <- here(laravel_storage_dir, "Data") extracted_CI_dir <- here(data_dir, "extracted_ci") daily_CI_vals_dir <- here(extracted_CI_dir, "daily_vals") cumulative_CI_vals_dir <- here(extracted_CI_dir, "cumulative_vals") weekly_CI_mosaic <- here(laravel_storage_dir, "weekly_mosaic") daily_vrt <- here(data_dir, "vrt") harvest_dir <- here(data_dir, "HarvestData") source(here("r_app/parameters_project.R")) source(here("r_app/mosaic_creation_utils.R")) source("parameters_project.R") source("mosaic_creation_utils.R") # dir.create(here(laravel_storage_dir)) # dir.create(here(reports_dir)) # dir.create(here(data_dir)) # dir.create(here(extracted_CI_dir)) # dir.create(here(daily_CI_vals_dir)) # dir.create(here(cumulative_CI_vals_dir)) # dir.create(here(weekly_CI_mosaic)) # dir.create(here(daily_vrt)) # dir.create(merged_final) # dir.create(harvest_dir) # end_date <- lubridate::dmy("28-8-2024") week <- week(end_date) #weeks_ago = 0 # Creating weekly mosaic #dates <- date_list(weeks_ago) dates <- date_list(end_date, offset) file_name_tif <- as.character(args[4]) if (is.na(file_name_tif)) { file_name_tif <- paste0("week_", sprintf("%02d", dates$week), "_", dates$year, ".tif") } print(dates) print(file_name_tif) #load pivot geojson # pivot_sf_q <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% vect() vrt_files <- list.files(here(daily_vrt),full.names = T) vrt_list <- map(dates$days_filter, ~ vrt_files[grepl(pattern = .x, x = vrt_files)]) %>% compact() %>% flatten_chr() raster_files_final <- list.files(merged_final,full.names = T, pattern = ".tif") if (length(vrt_list) > 0) { print("vrt list made, preparing mosaic creation by counting cloud cover") total_pix_area <- rast(vrt_list[1]) %>% terra::subset(1) %>% setValues(1) %>% crop(field_boundaries, mask = TRUE) %>% global(., fun = "notNA") #%>% layer_5_list <- purrr::map(vrt_list, function(vrt_list) { rast(vrt_list[1]) %>% terra::subset(1) }) %>% rast() missing_pixels_count <- layer_5_list %>% global(., fun = "notNA") %>% mutate( total_pixels = total_pix_area$notNA, missing_pixels_percentage = round(100 - (( notNA / total_pix_area$notNA ) * 100)), thres_5perc = as.integer(missing_pixels_percentage < 5), thres_40perc = as.integer(missing_pixels_percentage < 45) ) index_5perc <- which(missing_pixels_count$thres_5perc == max(missing_pixels_count$thres_5perc)) index_40perc <- which(missing_pixels_count$thres_40perc == max(missing_pixels_count$thres_40perc)) ## Create mosaic if (sum(missing_pixels_count$thres_5perc) > 1) { message("More than 1 raster without clouds (<5%), max composite made") cloudy_rasters_list <- vrt_list[index_5perc] rsrc <- sprc(cloudy_rasters_list) x <- terra::mosaic(rsrc, fun = "max") names(x) <- c("Red", "Green", "Blue", "NIR", "CI") } else if (sum(missing_pixels_count$thres_5perc) == 1) { message("Only 1 raster without clouds (<5%)") x <- rast(vrt_list[index_5perc[1]]) names(x) <- c("Red", "Green", "Blue", "NIR", "CI") } else if (sum(missing_pixels_count$thres_40perc) > 1) { message("More than 1 image contains clouds, composite made of <40% cloud cover images") cloudy_rasters_list <- vrt_list[index_40perc] rsrc <- sprc(cloudy_rasters_list) x <- mosaic(rsrc, fun = "max") names(x) <- c("Red", "Green", "Blue", "NIR", "CI") } else if (sum(missing_pixels_count$thres_40perc) == 1) { message("Only 1 image available but contains clouds") x <- rast(vrt_list[index_40perc[1]]) names(x) <- c("Red", "Green", "Blue", "NIR", "CI") } else if (sum(missing_pixels_count$thres_40perc) == 0) { message("No cloud free images available, all images combined") rsrc <- sprc(vrt_list) x <- mosaic(rsrc, fun = "max") names(x) <- c("Red", "Green", "Blue", "NIR", "CI") } } else{ message("No images available this week, empty mosaic created") x <- rast(raster_files_final[1]) %>% setValues(0) %>% crop(field_boundaries, mask = TRUE) names(x) <- c("Red", "Green", "Blue", "NIR", "CI") } plot(x$CI, main = paste("CI map ", dates$week)) plotRGB(x, main = paste("RGB map ", dates$week)) file_path_tif <- here(weekly_CI_mosaic ,file_name_tif) writeRaster(x, file_path_tif, overwrite=TRUE) message("Raster written/made at: ", file_path_tif)