# utils for mosaic creation date_list <- function(end_date, offset){ offset <- as.numeric(offset) - 1 end_date <- as.Date(end_date) start_date <- end_date - lubridate::days(offset) week <- week(start_date) year <- year(start_date) days_filter <- seq(from = start_date, to = end_date, by = "day") return(list("week" = week, "year" = year, "days_filter" = days_filter)) } # 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))) # # } 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) } 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}") } right = function(text, num_char) { substr(text, nchar(text) - (num_char-1), nchar(text)) } # 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) } #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) # }