217 lines
6.2 KiB
R
217 lines
6.2 KiB
R
# 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
|
|
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("parameters_project.R")
|
|
source("utils_1.R")
|
|
|
|
dir.create(here(laravel_storage_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("20-6-2024")
|
|
week <- week(end_date)
|
|
|
|
|
|
#weeks_ago = 0
|
|
|
|
# Creating weekly mosaic
|
|
#dates <- date_list(weeks_ago)
|
|
dates <- date_list(end_date, offset)
|
|
print(dates)
|
|
#load pivot geojson
|
|
# pivot_sf_q <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% vect()
|
|
|
|
raster_files <- list.files(planet_tif_folder,full.names = T, pattern = ".tif")
|
|
|
|
filtered_files <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>%
|
|
compact() %>%
|
|
flatten_chr()
|
|
head(filtered_files)
|
|
|
|
# filtered_files <- raster_files #for first CI extraction
|
|
|
|
# create_mask_and_crop <- function(file, field_boundaries) {
|
|
# message("starting ", file)
|
|
# CI <- rast(file)
|
|
# # names(CI) <- c("green","nir")
|
|
# message("raster loaded")
|
|
#
|
|
# CI <- CI[[2]]/CI[[1]]-1
|
|
# # CI <- CI$nir/CI$green-1
|
|
# message("CI calculated")
|
|
# CI <- terra::crop(CI, field_boundaries, mask = TRUE) #%>% CI_func()
|
|
#
|
|
# new_file <- here(merged_final, paste0(tools::file_path_sans_ext(basename(file)), ".tif"))
|
|
# writeRaster(CI, new_file, overwrite = TRUE)
|
|
#
|
|
# vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt"))
|
|
# terra::vrt(new_file, vrt_file, overwrite = TRUE)
|
|
#
|
|
# return(CI)
|
|
# }
|
|
|
|
|
|
vrt_list <- list()
|
|
|
|
# for (file in raster_files) {
|
|
# v_crop <- create_mask_and_crop(file, field_boundaries)
|
|
# message(file, " processed")
|
|
# gc()
|
|
# }
|
|
|
|
for (file in filtered_files) {
|
|
v_crop <- create_mask_and_crop(file, field_boundaries)
|
|
emtpy_or_full <- global(v_crop, "notNA")
|
|
|
|
vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt"))
|
|
|
|
if(emtpy_or_full[1,] > 10000){
|
|
vrt_list[vrt_file] <- vrt_file
|
|
|
|
}else{
|
|
file.remove(vrt_file)
|
|
|
|
}
|
|
|
|
message(file, " processed")
|
|
gc()
|
|
}
|
|
|
|
# list_global <- list_global %>% flatten_chr()
|
|
vrt_list <- vrt_list %>% flatten_chr()
|
|
|
|
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 <- mosaic(rsrc, fun = "max")
|
|
|
|
# names(x) <- "CI"
|
|
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("CI")
|
|
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) <- "CI"
|
|
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("CI")
|
|
names(x) <- c("Red", "Green", "Blue", "NIR", "CI")
|
|
}else{
|
|
message("No cloud free images available, all images combined")
|
|
|
|
rsrc <- sprc(vrt_list)
|
|
x <- mosaic(rsrc, fun = "max")
|
|
# x <- rast(vrt_list[1]) %>% setValues(NA)
|
|
# names(x) <- c("CI")
|
|
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 ,paste0("week_", sprintf("%02d", dates$week), "_", dates$year, ".tif"))
|
|
writeRaster(x, file_path_tif, overwrite=TRUE)
|
|
message("Raster written/made at: ", file_path_tif)
|