# activeer de renv omgeving; renv::activate('~/smartCane/r_app') renv::restore() library(here) library(sf) library(terra) library(tidyverse) library(lubridate) library(exactextractr) library(CIprep) laravel_storage_dir <- here("../laravel_app/storage/app") #preparing directories planet_tif_folder <- here(laravel_storage_dir, "chemba/merged_tif") merged_final <- here(laravel_storage_dir, "chemba/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(data_dir, "weekly_mosaic") daily_vrt <- here(data_dir, "vrt") harvest_dir <- here(data_dir, "HarvestData") 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) # Creating weekly mosaic dates <- date_list(0) #load pivot geojson pivot_sf_q <- st_read(here( "pivot_20210625.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% vect() #load and filter raster files raster_files <- list.files(planet_tif_folder,full.names = T, pattern = ".tif") #use pattern = '.tif$' or something else if you have multiple files in this folder filtered_files <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>% compact() %>% flatten_chr() #rasters_masked <- map(filtered_files, mask_raster, fields = pivot_sf_q) %>% set_names(filtered_files) # rasters_masked <- list() # Creëer een lege lijst om de resultaten op te slaan # for (i in seq_along(filtered_files[1])) { # file_name <- filtered_files[i] # result <- mask_raster(file_name, fields = pivot_sf_q) # rasters_masked[[file_name]] <- result # } create_mask_and_crop <- function(file, pivot_sf_q) { message("starting ", file) CI <- rast(file) # names(CI) <- c("green","nir") message("raster loaded") print(CI) CI <- CI[[2]]/CI[[1]]-1 # CI <- CI$nir/CI$green-1 message("CI calculated") CI <- terra::crop(CI, pivot_sf_q, mask = TRUE) #%>% CI_func() # v_crop[v_crop == 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(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) # v_crop <- mask_raster(v, pivot_sf_q) return(CI) } # rasters_masked <- map(filtered_files, create_mask_and_crop, pivot_sf_q) # list_global <- list() vrt_list <- list() for (file in filtered_files) { v_crop <- create_mask_and_crop(file, pivot_sf_q) 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){ # list_global[file] <- file vrt_list[vrt_file] <- vrt_file }else{ file.remove(vrt_file) # file.remove(file) # message(file, " removed") } # Save the processed raster to a new file # output_file <- here(data_dir, "vrt", paste0(basename(file), "_processed")) # terra::writeRaster(v_crop, output_file, overwrite = TRUE) message(file, " processed") gc() } # list_global <- list_global %>% flatten_chr() vrt_list <- vrt_list %>% flatten_chr() #testing writing raster # for(i in seq_along(filtered_files)){ # message("starting ", i) # x <- mask_raster(filtered_files[i], fields = pivot_sf_q) # writeRaster(x, filtered_files[i], overwrite=TRUE) # message("writing ", i) # } # rasters_masked[sapply(rasters_masked, is.null)] <- NULL # rasters_masked <- setNames(list_global, map_chr(names(list_global), date_extract)) # }) total_pix_area <- rast(vrt_list[1]) %>% terra::subset(1) %>% setValues(1) %>% crop(pivot_sf_q, mask = TRUE) %>% global(., fun="notNA") #%>% # as.matrix() %>% # `[`(1, 1) # total_pix_area <- rast(filtered_files[1]) %>% subset(1) %>% crop(pivot_sf_q, mask = TRUE) %>% freq(., usenames = TRUE) # list_global # rast(list_global[1])[[5]] %>% plot() # vrt_files <- list.files(here(data_dir, "vrt"),full.names = T) # vrt_days_filter <- tools::file_path_sans_ext(basename(list_global)) # vrt_list <- map(vrt_days_filter, ~ vrt_files[grepl(pattern = .x, x = vrt_files)]) %>% # compact() %>% # flatten_chr() 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 < 40) ) # cloud_perc_list <- freq(layer_5_list, usenames = TRUE) %>% # mutate(cloud_perc = (100 -((count/sum(total_pix_area$notNA))*100)), # cloud_thres_5perc = as.integer(cloud_perc < 5), # cloud_thres_40perc = as.integer(cloud_perc < 40)) %>% # rename(Date = layer) %>% select(-value, -count) 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" }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") }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_5perc] rsrc <- sprc(cloudy_rasters_list) x <- mosaic(rsrc, fun = "max") names(x) <- "CI" }else if(sum(missing_pixels_count$thres_40perc)==1){ message("Only 1 image available but contains clouds") x <- rast(vrt_list[index_5perc[1]]) names(x) <- c("CI") }else{ message("No cloud free images available") x <- rast(vrt_list[1]) %>% setValues(NA) names(x) <- c("CI") } plot(x$CI, main = paste("CI map", dates$week)) #plotRGB(x, main = paste("RGB image week", dates$week)) # writeRaster(x, here(weekly_CI_mosaic ,paste0("week_", dates$week, "_", dates$year, ".tif")), overwrite=TRUE) writeRaster(x, here(weekly_CI_mosaic ,paste0("week_", dates$week, "_", dates$year, ".tif")), overwrite=TRUE) message("raster written/ made") # 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, 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_20210625.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") 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= pivot_sf_q, quadrants = TRUE, daily_CI_vals_dir) message("after walk") 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") ) # pivots_dates_long <- pivots_dates0 %>% # select(c("pivot_quadrant", "season_start_2021", "season_end_2021", "season_start_2022", "season_end_2022", "season_start_2023", "season_end_2023")) %>% # pivot_longer(cols = c("season_start_2021", "season_end_2021", "season_start_2022", "season_end_2022", "season_start_2023", "season_end_2023")) %>% # separate(pivot_quadrant, into = c("name", "Year"), sep = "\\.") 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) # extracted_values <- list.files("C:\\Users\\timon\\Resilience BV\\4002 CMD App - General\\4002 CMD Team\\4002 TechnicalData\\04 WP2 technical\\DetectingSpotsR\\EcoFarm\\planet\\extracted", # pattern ="_quadrant", full.names = TRUE) extracted_values <- list.files(here(daily_CI_vals_dir), full.names = TRUE) #get CI values for this week only #extracted_values <- map(dates$days_filter, ~ extracted_values[grepl(pattern = .x, x = extracted_values)]) %>% # compact() %>% # flatten_chr() #combine them into one df pivot_stats <- extracted_values %>% map(readRDS) %>% list_rbind() %>% group_by(pivot_quadrant) %>% 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 #load historic CI data and update it with last week of CI data 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 <- purrr::map(list.files(here(daily_CI_vals_dir), full.names = TRUE, pattern = "quadrant"), readRDS) %>% list_rbind() %>% group_by(pivot_quadrant) %>% # summarise(across(everything(), ~ first(na.omit(.)))) pivots_data_present <- unique(pivots_dates0$pivot_quadrant) quadrant_list <- pivots_data_present # gather data into long format for easier calculation and visualisation pivot_stats_long <- pivot_stats2 %>% tidyr::gather("Date", value, -pivot_quadrant, -pivot ) %>% 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) # #2021 pivots_dates_Data_2021 <- pivots_dates0 %>% filter(!is.na(season_start_2021)) pivot_select_model_Data_2021 <- unique(pivots_dates_Data_2021$pivot_quadrant) # #2022 pivots_dates_Data_2022 <- pivots_dates0 %>% filter(!is.na(season_end_2022)) pivot_select_model_Data_2022 <- unique(pivots_dates_Data_2022$pivot_quadrant ) # #2023 pivots_dates_Data_2023 <- pivots_dates0 %>% filter(!is.na(season_start_2023)) pivot_select_model_Data_2023 <- unique(pivots_dates_Data_2023$pivot_quadrant) # #2024 pivots_dates_Data_2024 <- pivots_dates0 %>% filter(!is.na(season_start_2024)) pivot_select_model_Data_2024 <- unique(pivots_dates_Data_2024$pivot_quadrant) ## 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_2022, 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')