# Utils for ci extraction date_list <- function(end_date, offset, week){ 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))) # # } 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]) 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) } 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) # }