285 lines
10 KiB
R
285 lines
10 KiB
R
# Utils for ci extraction
|
|
|
|
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)))
|
|
#
|
|
# }
|
|
|
|
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)
|
|
# }
|
|
|