# activeer de renv omgeving; # renv::activate('~/smartCane/r_app') # renv::restore() library(here) library(sf) library(terra) library(tidyverse) library(lubridate) library(exactextractr) #funcion CI_prep package date_list <- function(weeks_in_the_paste){ week <- week(Sys.Date()- weeks(weeks_in_the_paste) ) year <- year(Sys.Date()- weeks(weeks_in_the_paste) ) days_filter <- Sys.Date() - weeks(weeks_in_the_paste) - days(0:6) return(c("week" = week, "year" = year, "days_filter" = list(days_filter))) } CI_func <- function(x, drop_layers = FALSE){ CI <- x[[4]]/x[[2]]-1 add(x) <- CI names(x) <- c("red", "green", "blue","nir", "cloud" ,"CI") if(drop_layers == FALSE){ return(x) }else{ return(x$CI) } } mask_raster <- function(raster, fields){ # x <- rast(filtered_files[1]) x <- rast(raster) emtpy_or_full <- global(x, sum) if(emtpy_or_full[1,] >= 2000000){ names(x) <- c("red", "green", "blue","nir", "cloud") cloud <- x$cloud cloud[cloud == 0 ] <- NA x_masked <- mask(x, cloud, inverse = T) %>% crop(.,fields, mask = TRUE ) x_masked <- x_masked %>% CI_func() message(raster, " masked") return(x_masked) } } date_extract <- function(file_path) { str_extract(file_path, "\\d{4}-\\d{2}-\\d{2}") } extract_rasters_daily <- function(file, field_geojson, quadrants = TRUE, save_dir) { # x <- rast(filtered_files[1])%>% CI_func(drop_layers = TRUE) # date <- date_extract(filtered_files[1]) # field_geojson <- sf::st_as_sf(pivot_sf_q) field_geojson <- sf::st_as_sf(field_geojson) x <- rast(file[1]) %>% CI_func(drop_layers = TRUE) date <- date_extract(file) pivot_stats <- cbind(field_geojson, mean_CI = round(exactextractr::exact_extract(x, field_geojson, fun = "mean"), 2)) %>% st_drop_geometry() %>% rename("{date}" := mean_CI) save_suffix <- if (quadrants){"quadrant"} else {"whole_field"} save_path <- here(save_dir, paste0("extracted_", date, "_", save_suffix, ".rds")) saveRDS(pivot_stats, save_path) } right = function(text, num_char) { substr(text, nchar(text) - (num_char-1), nchar(text)) } extract_CI_data <- function(field_names, harvesting_data, field_CI_data, season) { # field_names = "1.2A" # harvesting_data = harvesting_data # field_CI_data = pivot_stats_long # season= 2021 filtered_harvesting_data <- harvesting_data %>% filter(Year == season, Field %in% field_names) filtered_field_CI_data <- field_CI_data %>% filter(Field %in% field_names) # CI <- map_df(field_names, ~ { ApproxFun <- approxfun(x = filtered_field_CI_data$Date, y = filtered_field_CI_data$value) Dates <- seq.Date(ymd(min(filtered_field_CI_data$Date)), ymd(max(filtered_field_CI_data$Date)), by = 1) LinearFit <- ApproxFun(Dates) CI <- data.frame(Date = Dates, FitData = LinearFit) %>% left_join(., filtered_field_CI_data, by = "Date") %>% filter(Date > filtered_harvesting_data$Season_start & Date < filtered_harvesting_data$Season_end) %>% mutate(DOY = seq(1, n(), 1), model = paste0("Data", season, " : ", field_names), season = season, Field = field_names) # }) #%>% #{if (length(field_names) > 0) message("Done!")} return(CI) } load_fields <- function(geojson_path) { field_geojson <- st_read(geojson_path) %>% select(pivot, pivot_quadrant) %>% vect() return(field_geojson) } load_harvest_data <- function(havest_data_path){ harvest_data <- readRDS(havest_data_path) return(harvest_data) } load_rasters <- function(raster_path, dates) { raster_files <- list.files(raster_path, full.names = TRUE, pattern = ".tif") filtered_files <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>% compact() %>% flatten_chr() return(filtered_files) } mask_and_set_names <- function(filtered_files, fields) { rasters_masked <- map(filtered_files, mask_raster, fields = fields) %>% set_names(filtered_files) rasters_masked[sapply(rasters_masked, is.null)] <- NULL rasters_masked <- setNames(rasters_masked, map_chr(names(rasters_masked), date_extract)) return(rasters_masked) } calculate_total_pix_area <- function(filtered_files, fields_geojson) { # total_pix_area <- rast(filtered_files[1]) %>% # subset(1) %>% # crop(fields_geojson, mask = TRUE)%>% # global(.data, fun = "notNA") total_pix_area <- rast(filtered_files[1]) %>% subset(1) %>% crop(fields_geojson, mask = TRUE) %>% freq(., usenames = TRUE) return(total_pix_area) } cloud_layer_extract <- function(rasters_masked){ cloud_layer_rast <- map(rasters_masked, function(spatraster) { spatraster[[5]] }) %>% rast() return(cloud_layer_rast) } calculate_cloud_coverage <- function(cloud_layer_rast, total_pix_area) { cloud_perc_list <- freq(cloud_layer_rast, usenames = TRUE) %>% mutate(cloud_perc = (100 -((count/total_pix_area$count)*100)), cloud_thres_5perc = as.integer(cloud_perc < 5), cloud_thres_40perc = as.integer(cloud_perc < 40)) %>% rename(Date = layer) %>% select(-value, -count) cloud_index_5perc <- which(cloud_perc_list$cloud_thres_5perc == max(cloud_perc_list$cloud_thres_5perc)) cloud_index_40perc <- which(cloud_perc_list$cloud_thres_40perc == max(cloud_perc_list$cloud_thres_40perc)) return(list(cloud_perc_list = cloud_perc_list, cloud_index_5perc = cloud_index_5perc, cloud_index_40perc = cloud_index_40perc)) } process_cloud_coverage <- function(cloud_coverage, rasters_masked) { if (sum(cloud_coverage$cloud_perc_list$cloud_thres_5perc) > 1) { message("More than 1 raster without clouds (<5%), max mosaic made ") cloudy_rasters_list <- rasters_masked[cloud_coverage$cloud_index_5perc] rsrc <- sprc(cloudy_rasters_list) x <- mosaic(rsrc) names(x) <- c("red", "green", "blue", "nir", "cloud", "CI") } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_5perc) == 1) { message("Only 1 raster without clouds (<5%)") x <- rast(rasters_masked[cloud_coverage$cloud_index_5perc[1]]) names(x) <- c("red", "green", "blue", "nir", "cloud", "CI") } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_40perc) > 1) { message("More than 1 image contains clouds, composite made of <40% cloud cover images") cloudy_rasters_list <- rasters_masked[cloud_coverage$cloud_index_40perc] rsrc <- sprc(cloudy_rasters_list) x <- mosaic(rsrc) names(x) <- c("red", "green", "blue", "nir", "cloud", "CI") } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_40perc) == 0) { message("No cloud free images available") x <- rast(rasters_masked[1]) x[x] <- NA names(x) <- c("red", "green", "blue", "nir", "cloud", "CI") } return(x) } extract_rasters_daily_func <- function(daily_vals_dir, filtered_files, fields_geojson) { extracted_files <- walk(filtered_files, extract_rasters_daily, field_geojson = fields_geojson, quadrants = TRUE, daily_vals_dir) return(extracted_files) } CI_load <- function(daily_vals_dir, grouping_variable){ extracted_values <- list.files(here(daily_vals_dir), full.names = TRUE) field_CI_values <- extracted_values %>% map_dfr(readRDS) %>% group_by(.data[[grouping_variable]]) %>% summarise(across(everything(), ~ first(na.omit(.)))) return(field_CI_values) } CI_long <- function(field_CI_values, pivot_long_cols){ field_CI_long <- field_CI_values %>% gather("Date", value, -pivot_long_cols) %>% mutate(Date = right(Date, 8), Date = ymd(Date) ) %>% drop_na(c("value","Date")) %>% mutate(value = as.numeric(value))%>% filter_all(all_vars(!is.infinite(.)))%>% rename(Field = pivot_quadrant) return(field_CI_long) } process_year_data <- function(year, harvest_data, field_CI_long) { pivots_dates_year <- harvest_data %>% na.omit() %>% filter(Year == year) pivot_select_model_year <- unique(pivots_dates_year$Field) data <- map_dfr(pivot_select_model_year, ~ extract_CI_data(.x, harvest_data, field_CI_long, season = year)) return(data) } #functions for CI_data_prep create_mask_and_crop <- function(file, pivot_sf_q) { # file <- filtered_files[5] message("starting ", file) loaded_raster <- rast(file) names(loaded_raster) <- c("Red", "Green", "Blue", "NIR") # names(CI) <- c("green","nir") message("raster loaded") # CI <- CI[[2]] / CI[[1]] - 1 CI <- loaded_raster$NIR / loaded_raster$Green - 1 loaded_raster$CI <- CI # CI <- CI$nir/CI$green-1 message("CI calculated") loaded_raster <- terra::crop(loaded_raster, pivot_sf_q, mask = TRUE) #%>% CI_func() loaded_raster[loaded_raster == 0] <- NA # names(v_crop) <- c("red", "green", "blue","nir", "cloud" ,"CI") # v_crop$CI <- v_crop$CI - 1 new_file <- here(merged_final, paste0(tools::file_path_sans_ext(basename(file)), ".tif")) writeRaster(loaded_raster, 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) # v_crop <- mask_raster(v, pivot_sf_q) return(loaded_raster) } extract_rasters_daily <- function(file, field_geojson, quadrants = TRUE, save_dir) { # x <- rast(filtered_files[1])%>% CI_func(drop_layers = TRUE) # date <- date_extract(filtered_files[1]) # field_geojson <- sf::st_as_sf(pivot_sf_q) field_geojson <- sf::st_as_sf(field_geojson) x <- rast(file[1]) date <- date_extract(file) pivot_stats <- cbind(field_geojson, mean_CI = round(exactextractr::exact_extract(x$CI, field_geojson, fun = "mean"), 2)) %>% st_drop_geometry() %>% rename("{date}" := mean_CI) save_suffix <- if (quadrants){"quadrant"} else {"whole_field"} save_path <- here(save_dir, paste0("extracted_", date, "_", save_suffix, ".rds")) saveRDS(pivot_stats, save_path) } #functions for rmarkdown file create_RGB_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = TRUE, legend_is_portrait = FALSE, week, age, red = TRUE) { r <- if (red) 1 else 4 # Set r based on the value of red title <- if (red) paste0("RGB image of the fields") else paste0("False colour image of the fields") tm_shape(pivot_raster, unit = "m") + tm_rgb(r = r, g = 2, b = 3, max.value = 255) + tm_layout(main.title = title, main.title.size = 1) + tm_scale_bar(position = c("right", "top"), text.color = "black") + tm_compass(position = c("right", "top"), text.color = "black") + tm_shape(pivot_shape) + tm_borders(col = "gray") + tm_text("subField", size = 0.6, col = "gray") + tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha = 0.5) } create_CI_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = F, legend_is_portrait = F, week, age, legend_only = F){ tm_shape(pivot_raster, unit = "m")+ tm_raster(breaks = CI_breaks, palette = "RdYlGn",legend.is.portrait = legend_is_portrait ,midpoint = NA) + tm_layout(main.title = paste0("Max CI week ", week,"\n", age, " weeks old"), main.title.size = 1, legend.show = show_legend, legend.only = legend_only) + tm_shape(pivot_shape) + tm_borders(lwd = 3) + tm_text("subField", size = 1/2) + tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha=0.5) +tmap_options(check.and.fix = TRUE) } create_CI_diff_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = F, legend_is_portrait = F, week_1, week_2, age){ tm_shape(pivot_raster, unit = "m")+ tm_raster(breaks = CI_diff_breaks, palette = "PRGn",legend.is.portrait = legend_is_portrait ,midpoint = 0, title = "CI difference") + tm_layout(main.title = paste0("CI change week ", week_1, "- week ",week_2, "\n", age," weeks old"), main.title.size = 1, legend.show = show_legend) + tm_shape(pivot_shape) + tm_borders(lwd = 3) + tm_text("subField", size = 1/2) + tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha=0.5) } ci_plot <- function(pivotName){ # pivotName = "MV2B09" # pivotName = "1.1" pivotShape <- AllPivots_merged %>% terra::subset(Field %in% pivotName) %>% st_transform(crs(CI)) # age <- AllPivots %>% dplyr::filter(Field %in% pivotName) %>% st_drop_geometry() %>% dplyr::select(Age) %>% unique() %>% # mutate(Age = round(Age)) age <- AllPivots %>% group_by(Field) %>% filter(Season == max(Season, na.rm = TRUE), Field %in% pivotName) %>% dplyr::select(Age)%>% st_drop_geometry() %>% unique() AllPivots2 <- AllPivots0 %>% dplyr::filter(Field %in% pivotName) singlePivot <- CI %>% crop(., pivotShape) %>% mask(., pivotShape) singlePivot_m1 <- CI_m1 %>% crop(., pivotShape) %>% mask(., pivotShape) singlePivot_m2 <- CI_m2 %>% crop(., pivotShape) %>% mask(., pivotShape) # singlePivot_m3 <- CI_m3 %>% crop(., pivotShape) %>% mask(., pivotShape) singlePivot_RGB <- RGB_raster %>% crop(., pivotShape) %>% mask(., pivotShape) singlePivot_false <- RGB_raster_stretch %>% crop(., pivotShape) %>% mask(., pivotShape) abs_CI_last_week <- last_week_dif_raster_abs %>% crop(., pivotShape) %>% mask(., pivotShape) abs_CI_three_week <- three_week_dif_raster_abs %>% crop(., pivotShape) %>% mask(., pivotShape) # planting_date <- harvesting_data %>% dplyr::filter(Field %in% pivotName) %>% ungroup() %>% dplyr::select(planting_date) %>% unique() joined_spans2 <- joined_spans %>% st_transform(crs(pivotShape)) %>% dplyr::filter(Field %in% pivotName) %>% st_crop(., pivotShape) # CImap_m2 <- create_CI_map(singlePivot_m2, AllPivots2, joined_spans2, show_legend= T, legend_is_portrait = T, week = week_minus_2, age = age -2) Legend_map <- create_CI_map(singlePivot_m1, AllPivots2, joined_spans2, show_legend= T, legend_is_portrait =T, week = week_minus_1, age = age -1, legend_only = T) CImap_m1 <- create_CI_map(singlePivot_m1, AllPivots2, joined_spans2, show_legend= T, legend_is_portrait =T, week = week_minus_1, age = age -1) CImap <- create_CI_map(singlePivot, AllPivots2, joined_spans2, show_legend= F, legend_is_portrait = T, week = week, age = age ) RGBmap <- create_RGB_map(singlePivot_false, AllPivots2, joined_spans2, show_legend= F, week = week, age = age, red =T ) Falsemap <- create_RGB_map(singlePivot_false, AllPivots2, joined_spans2, show_legend= F, week = week, age = age, red =F ) CI_max_abs_last_week <- create_CI_diff_map(abs_CI_last_week,AllPivots2, joined_spans2, show_legend = T, legend_is_portrait = T, week_1 = week, week_2 = week_minus_1, age = age) CI_max_abs_three_week <- create_CI_diff_map(abs_CI_three_week, AllPivots2, joined_spans2, show_legend = F, legend_is_portrait = T, week_1 = week, week_2 = week_minus_3, age = age) # tst <- tmap_arrange(CImap_m2, CImap_m1, CImap,CI_max_abs_last_week, CI_max_abs_three_week, nrow = 1) tst <- tmap_arrange(RGBmap,Falsemap, CImap_m1, CImap, CI_max_abs_last_week, CI_max_abs_three_week, ncol = 2) cat(paste("## Field", pivotName, "-", age$Age[1], "weeks after planting/harvest", "\n")) # cat("\n") # cat('

Pivot', pivotName, '- week', week, '-', age$Age, 'weeks after planting/harvest

') # cat(paste("# Pivot",pivots$pivot[i],"\n")) print(tst) } subchunkify <- function(g, fig_height=7, fig_width=5) { g_deparsed <- paste0(deparse( function() {g} ), collapse = '') sub_chunk <- paste0("```{r sub_chunk_", floor(runif(1) * 10000), ", fig.height=", fig_height, ", fig.width=", fig_width, ", echo=FALSE}", "\n(", g_deparsed , ")()", "\n``` ") cat(knitr::knit(text = knitr::knit_expand(text = sub_chunk), quiet = TRUE)) } cum_ci_plot <- function(pivotName){ # pivotName = "2.1" # Check if pivotName exists in the data if (!pivotName %in% unique(CI_quadrant$Field)) { # message("PivotName '", pivotName, "' not found. Plotting empty graph.") g <- ggplot() + labs(title = "Empty Graph - Yield dates missing") return( subchunkify(g, fig_height=2, fig_width = 10) ) } else { # message("PivotName '", pivotName, "' found. Plotting normal graph.") data_ci <- CI_quadrant %>% filter(Field %in% pivotName) data_ci2 <- data_ci %>% mutate(CI_rate = cumulative_CI/DOY, week = week(Date))%>% group_by(subField) %>% mutate(mean_rolling10 = rollapplyr(CI_rate , width = 10, FUN = mean, partial = TRUE)) # date_preperation_perfect_pivot <- data_ci2 %>% group_by(season) %>% summarise(min_date = min(Date), # max_date = max(Date), # days = max_date - min_date) # Identify unique seasons filtered_data <- data_ci2 %>% group_by(season) %>% mutate(rank = dense_rank(desc(season))) %>% filter(rank <= 2) %>% ungroup() %>% dplyr::select(-rank) # g <- ggplot(data= data_ci2 %>% filter(season %in% unique_seasons)) + g <- ggplot(data= filtered_data ) + # geom_line(aes(Date, mean_rolling10, col = subField)) + geom_line(aes(Date, CI_rate, col = subField)) + facet_wrap(~season, scales = "free_x") + # geom_line(data= perfect_pivot, aes(Date , mean_rolling10, col = "Model CI (p5.1 Data 2022, \n date x axis is fictive)"), lty="11",size=1) + labs(title = paste("CI rate - Field", pivotName), y = "CI rate (cumulative CI / Age)")+ # scale_y_continuous(limits=c(0.5,3), breaks = seq(0.5, 3, 0.5))+ scale_x_date(date_breaks = "1 month", date_labels = "%m-%Y") + theme(axis.text.x = element_text(angle = 60, hjust = 1), legend.justification=c(1,0), legend.position = c(1, 0), legend.title = element_text(size = 8), legend.text = element_text(size = 8)) + guides(color = guide_legend(nrow = 2, byrow = TRUE)) subchunkify(g, fig_height=6, fig_width = 10) } } # 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 weeks_ago <- as.numeric(args[1]) # Controleer of weeks_ago een geldig getal is if (is.na(weeks_ago)) { # stop("Het argument is geen geldig getal") weeks_ago <- 0 } # Converteer het tweede argument naar een string waarde project_dir <- as.character(args[2]) # 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") 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) #weeks_ago = 0 # Creating weekly mosaic dates <- date_list(weeks_ago) 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) # } create_mask_and_crop <- function(file, field_boundaries) { # file <- filtered_files[5] message("starting ", file) loaded_raster <- rast(file) names(loaded_raster) <- c("Red", "Green", "Blue", "NIR") # names(CI) <- c("green","nir") message("raster loaded") # CI <- CI[[2]] / CI[[1]] - 1 CI <- loaded_raster$NIR / loaded_raster$Green - 1 loaded_raster$CI <- CI # CI <- CI$nir/CI$green-1 message("CI calculated") loaded_raster <- terra::crop(loaded_raster, field_boundaries, mask = TRUE) #%>% CI_func() loaded_raster[loaded_raster == 0] <- NA # names(v_crop) <- c("red", "green", "blue","nir", "cloud" ,"CI") # v_crop$CI <- v_crop$CI - 1 new_file <- here(merged_final, paste0(tools::file_path_sans_ext(basename(file)), ".tif")) writeRaster(loaded_raster, 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) # v_crop <- mask_raster(v, pivot_sf_q) return(loaded_raster) } 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) # Extracting CI extract_rasters_daily <- function(file, field_geojson, quadrants = TRUE, save_dir) { # x <- rast(filtered_files[1])%>% CI_func(drop_layers = TRUE) # date <- date_extract(filtered_files[1]) # field_geojson <- sf::st_as_sf(pivot_sf_q) field_geojson <- sf::st_as_sf(field_geojson) x <- rast(file[1]) date <- date_extract(file) pivot_stats <- cbind(field_geojson, mean_CI = round(exactextractr::exact_extract(x$CI, field_geojson, fun = "mean"), 2)) %>% st_drop_geometry() %>% rename("{date}" := mean_CI) save_suffix <- if (quadrants){"quadrant"} else {"whole_field"} save_path <- here(save_dir, paste0("extracted_", date, "_", save_suffix, ".rds")) saveRDS(pivot_stats, save_path) } # pivot_sf_q <- st_read(here("..", "Data", "pivot_20210625.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% vect() # pivot_sf <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% group_by(pivot) %>% summarise() %>% vect() # message("pivot loaded") raster_files_NEW <- list.files(merged_final,full.names = T, pattern = ".tif") # pivots_dates0 <- readRDS(here(harvest_dir, "harvest_data_new")) %>% filter( # pivot %in% c("1.1", "1.2", "1.3", "1.4", "1.6", "1.7", "1.8", "1.9", "1.10", "1.11", "1.12", "1.13", # "1.14" , "1.16" , "1.17" , "1.18" ,"2.1", "2.2", "2.3" , "2.4", "2.5", "3.1", "3.2", "3.3", # "4.1", "4.2", "4.3", "4.4", "4.5", "4.6", "5.1" ,"5.2", "5.3", "5.4", "6.1", "6.2", "DL1.1", "DL1.3") # ) # harvesting_data <- pivots_dates0 %>% # select(c("pivot_quadrant", "season_start_2021", "season_end_2021", "season_start_2022", "season_end_2022", "season_start_2023", "season_end_2023", "season_start_2024", "season_end_2024")) %>% # pivot_longer(cols = starts_with("season"), names_to = "Year", values_to = "value") %>% # separate(Year, into = c("name", "Year"), sep = "(?<=season_start|season_end)\\_", remove = FALSE) %>% # mutate(name = str_to_title(name)) %>% # pivot_wider(names_from = name, values_from = value) %>% # rename(Field = pivot_quadrant) #If run for the firsttime, it will extract all data since the start and put into a table.rds. otherwise it will only add on to the existing table. if (new_project_question == TRUE) { print("combined_CI_data.rds does not exist. Preparing combined_CI_data.rds file for all available images.") walk(raster_files_NEW, extract_rasters_daily, field_geojson= field_boundaries, quadrants = F, daily_CI_vals_dir) extracted_values <- list.files(here(daily_CI_vals_dir), full.names = TRUE) pivot_stats <- extracted_values %>% map(readRDS) %>% list_rbind() %>% group_by(subField) %>% summarise(across(everything(), ~ first(na.omit(.)))) saveRDS(pivot_stats, here(cumulative_CI_vals_dir,"combined_CI_data.rds")) #used to save the rest of the data into one file pivot_stats2 <- pivot_stats print("All CI values extracted from all historic images") } else { print("combined_CI_data.rds exists, adding the latest image data to the table.") filtered_files <- map(dates$days_filter, ~ raster_files_NEW[grepl(pattern = .x, x = raster_files_NEW)]) %>% compact() %>% flatten_chr() walk(filtered_files, extract_rasters_daily, field_geojson= field_boundaries, quadrants = TRUE, daily_CI_vals_dir) extracted_values <- list.files(daily_CI_vals_dir, full.names = TRUE) extracted_values <- map(dates$days_filter, ~ extracted_values[grepl(pattern = .x, x = extracted_values)]) %>% compact() %>% flatten_chr() pivot_stats <- extracted_values %>% map(readRDS) %>% list_rbind() %>% group_by(subField) %>% summarise(across(everything(), ~ first(na.omit(.)))) combined_CI_data <- readRDS(here(cumulative_CI_vals_dir,"combined_CI_data.rds")) #%>% drop_na(pivot_quadrant) pivot_stats2 <- bind_rows(pivot_stats, combined_CI_data) # pivot_stats2 <- combined_CI_data print("All CI values extracted from latest 7 images.") saveRDS(combined_CI_data, here(cumulative_CI_vals_dir,"combined_CI_data.rds")) #used to save the rest of the data into one file } # gather data into long format for easier calculation and visualisation pivot_stats_long <- pivot_stats2 %>% tidyr::gather("Date", value, -Field, -subField, -sub_area ) %>% mutate(Date = right(Date, 8), Date = lubridate::ymd(Date) ) %>% drop_na(c("value","Date")) %>% mutate(value = as.numeric(value))%>% filter_all(all_vars(!is.infinite(.))) %>% # rename(Field = pivot_quadrant, # subField = Field) %>% unique() # #2021 # pivot_select_model_Data_2021 <- harvesting_data %>% filter(Year == 2021) %>% pull(Field) # # pivot_select_model_Data_2022 <- harvesting_data %>% filter(Year == 2022) %>% pull(Field) pivot_select_model_Data_2023 <- harvesting_data %>% filter(Year == 2023) %>% filter(!is.na(Season_start)) %>% pull(subField) pivot_select_model_Data_2024 <- harvesting_data %>% filter(Year == 2024)%>% filter(!is.na(Season_start)) %>% pull(subField) # pivots_dates_Data_2022 <- pivots_dates0 %>% filter(!is.na(season_end_2022)) # pivot_select_model_Data_2022 <- unique(pivots_dates_Data_2022$pivot_quadrant ) # # pivots_dates_Data_2023 <- pivots_dates0 %>% filter(!is.na(season_start_2023)) # pivot_select_model_Data_2023 <- unique(pivots_dates_Data_2023$pivot_quadrant) # # pivots_dates_Data_2024 <- pivots_dates0 %>% filter(!is.na(season_start_2024)) # pivot_select_model_Data_2024 <- unique(pivots_dates_Data_2024$pivot_quadrant) extract_CI_data <- function(field_names, harvesting_data, field_CI_data, season) { # field_names = "AG1D06P" # field_names = "1.13A" # harvesting_data = harvesting_data # field_CI_data = pivot_stats_long # season= 2023 filtered_harvesting_data <- harvesting_data %>% filter(Year == season, subField %in% field_names) filtered_field_CI_data <- field_CI_data %>% filter(subField %in% field_names) # CI <- map_df(field_names, ~ { ApproxFun <- approxfun(x = filtered_field_CI_data$Date, y = filtered_field_CI_data$value) Dates <- seq.Date(ymd(min(filtered_field_CI_data$Date)), ymd(max(filtered_field_CI_data$Date)), by = 1) LinearFit <- ApproxFun(Dates) CI <- data.frame(Date = Dates, FitData = LinearFit) %>% left_join(., filtered_field_CI_data, by = "Date") %>% filter(Date > filtered_harvesting_data$Season_start & Date < filtered_harvesting_data$Season_end) %>% mutate(DOY = seq(1, n(), 1), model = paste0("Data", season, " : ", field_names), season = season, subField = field_names) # }) #%>% #{if (length(field_names) > 0) message("Done!")} return(CI) } ## Extracting the correct CI values # Data_2021 <- map(pivot_select_model_Data_2021, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2021)) %>% list_rbind() # message('2021') # Data_2022 <- map(pivot_select_model_Data_2022, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2022)) %>% list_rbind() # message('2022') Data_2023 <- map(pivot_select_model_Data_2023, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2023)) %>% list_rbind() message('2023') Data_2024 <- map(pivot_select_model_Data_2024, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2024)) %>% list_rbind() message('2024') CI_all <- rbind(Data_2023, Data_2024) message('CI_all created') #CI_all <- Data_2023 CI_all <- CI_all %>% group_by(model) %>% mutate(CI_per_day = FitData - lag(FitData), cumulative_CI = cumsum(FitData)) message('CI_all cumulative') head(CI_all) message('show head') saveRDS(CI_all, here(cumulative_CI_vals_dir,"All_pivots_Cumulative_CI_quadrant_year_v2.rds")) message('rds saved')