218 lines
7.2 KiB
R
218 lines
7.2 KiB
R
# 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)
|
|
# }
|
|
|
|
|