wip

laravel_app/app/Jobs/ProjectDownloadTiffJob.php

@@ -53,6 +53,7 @@ public function handle(): void
         if (!$process->isSuccessful()) {
             logger('error', $process->getErrorOutput());
         }
+        logger($process->getOutput());
 
         $this->download->setStatusSuccess();
     }

r_app/2_CI_data_prep.R

@@ -3,14 +3,453 @@
 renv::activate('~/smartCane/r_app')
 renv::restore()
 
+
+
 library(here)
 library(sf)
 library(terra)
 library(tidyverse)
 library(lubridate)
 library(exactextractr)
+#funcion CI_prep package
+
+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]) %>% CI_func(drop_layers = TRUE)
+  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)
+}
+
+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)
+}
+
+#functions for rmarkdown file
+
+
+create_RGB_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = TRUE, legend_is_portrait = FALSE, week, age, red = TRUE) {
+  r <- if (red) 1 else 4  # Set r based on the value of red
+  title <- if (red) paste0("RGB image of the fields") else paste0("False colour image of the fields")
+  
+  tm_shape(pivot_raster, unit = "m") + tm_rgb(r = r, g = 2, b = 3, max.value = 255) +
+    tm_layout(main.title = title, 
+              main.title.size = 1) +
+    tm_scale_bar(position = c("right", "top"), text.color = "black") +
+    tm_compass(position = c("right", "top"), text.color = "black") +
+    tm_shape(pivot_shape) + tm_borders(col = "gray") +
+    tm_text("subField", size = 0.6, col = "gray") +
+    tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha = 0.5)
+}
+
+create_CI_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = F, legend_is_portrait = F, week, age, legend_only = F){
+  tm_shape(pivot_raster, unit = "m")+
+    tm_raster(breaks =  CI_breaks, palette = "RdYlGn",legend.is.portrait = legend_is_portrait ,midpoint = NA) +
+    tm_layout(main.title = paste0("Max CI week ", week,"\n",  age, " weeks old"),
+              main.title.size = 1, legend.show = show_legend, legend.only = legend_only) +
+    tm_shape(pivot_shape) +
+    tm_borders(lwd = 3) +  tm_text("subField", size = 1/2) +
+    tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha=0.5) +tmap_options(check.and.fix = TRUE)
+}
+
+create_CI_diff_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = F, legend_is_portrait = F, week_1, week_2, age){
+  tm_shape(pivot_raster, unit = "m")+
+    tm_raster(breaks = CI_diff_breaks, palette = "PRGn",legend.is.portrait = legend_is_portrait ,midpoint = 0, title = "CI difference") +
+    tm_layout(main.title =  paste0("CI change week ", week_1, "- week ",week_2, "\n", age," weeks old"),
+              main.title.size = 1, legend.show = show_legend) +
+    tm_shape(pivot_shape) +
+    tm_borders(lwd = 3) +  tm_text("subField", size = 1/2) +
+    tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha=0.5)
+}
+
+ci_plot <- function(pivotName){
+  # pivotName = "MV2B09"
+  # pivotName = "1.1"
+  
+  pivotShape <- AllPivots_merged %>% terra::subset(Field %in% pivotName)  %>% st_transform(crs(CI))
+  # age <- AllPivots %>% dplyr::filter(Field %in% pivotName) %>% st_drop_geometry() %>% dplyr::select(Age) %>% unique() %>%
+  #   mutate(Age = round(Age))
+  
+  age <- AllPivots %>%
+    group_by(Field) %>%
+    filter(Season == max(Season, na.rm = TRUE), Field %in% pivotName) %>%
+    dplyr::select(Age)%>% st_drop_geometry()  %>% unique() 
+  
+  AllPivots2 <- AllPivots0 %>% dplyr::filter(Field %in% pivotName)
+  
+  singlePivot <- CI %>% crop(., pivotShape) %>% mask(., pivotShape)
+  singlePivot_m1 <- CI_m1 %>% crop(., pivotShape) %>% mask(., pivotShape)
+  singlePivot_m2 <- CI_m2 %>% crop(., pivotShape) %>% mask(., pivotShape)
+  # singlePivot_m3 <- CI_m3 %>% crop(., pivotShape) %>% mask(., pivotShape)
+  
+  singlePivot_RGB <- RGB_raster %>% crop(., pivotShape) %>% mask(., pivotShape)
+  singlePivot_false <- RGB_raster_stretch %>% crop(., pivotShape) %>% mask(., pivotShape)
+  
+  abs_CI_last_week <- last_week_dif_raster_abs %>% crop(., pivotShape) %>% mask(., pivotShape)
+  abs_CI_three_week <- three_week_dif_raster_abs %>% crop(., pivotShape) %>% mask(., pivotShape)
+  
+  # planting_date <- harvesting_data %>% dplyr::filter(Field %in% pivotName) %>% ungroup()  %>% dplyr::select(planting_date)  %>% unique()
+  
+  joined_spans2 <- joined_spans %>% st_transform(crs(pivotShape))  %>%  dplyr::filter(Field %in% pivotName) %>% st_crop(., pivotShape)
+  
+  # CImap_m2 <- create_CI_map(singlePivot_m2, AllPivots2,  joined_spans2, show_legend= T, legend_is_portrait = T, week = week_minus_2, age = age  -2)
+  Legend_map <- create_CI_map(singlePivot_m1, AllPivots2,  joined_spans2, show_legend= T, legend_is_portrait =T, week = week_minus_1, age = age  -1, legend_only = T)
+  CImap_m1 <- create_CI_map(singlePivot_m1, AllPivots2,  joined_spans2, show_legend= T, legend_is_portrait =T, week = week_minus_1, age = age  -1)
+  CImap <- create_CI_map(singlePivot, AllPivots2,  joined_spans2, show_legend= F, legend_is_portrait = T, week = week, age = age )
+  RGBmap <- create_RGB_map(singlePivot_false, AllPivots2,  joined_spans2, show_legend= F,  week = week, age = age, red =T )
+  Falsemap <- create_RGB_map(singlePivot_false, AllPivots2,  joined_spans2, show_legend= F,  week = week, age = age, red =F )
+  
+  
+  CI_max_abs_last_week <- create_CI_diff_map(abs_CI_last_week,AllPivots2,  joined_spans2,  show_legend = T, legend_is_portrait = T, week_1 = week, week_2 = week_minus_1, age = age)
+  CI_max_abs_three_week <- create_CI_diff_map(abs_CI_three_week, AllPivots2,  joined_spans2, show_legend = F, legend_is_portrait = T, week_1 = week, week_2 = week_minus_3, age = age)
+  
+  # tst <- tmap_arrange(CImap_m2, CImap_m1, CImap,CI_max_abs_last_week, CI_max_abs_three_week, nrow = 1)
+  tst <- tmap_arrange(RGBmap,Falsemap,
+                      CImap_m1, CImap, 
+                      CI_max_abs_last_week, CI_max_abs_three_week, 
+                      ncol = 2)
+  
+  cat(paste("## Field", pivotName, "-", age$Age[1], "weeks after planting/harvest", "\n"))
+  # cat("\n")
+  # cat('<h2> Pivot', pivotName, '- week', week, '-', age$Age, 'weeks after planting/harvest <h2>')
+  # cat(paste("# Pivot",pivots$pivot[i],"\n"))
+  print(tst)
+  
+}
+
+
+subchunkify <- function(g, fig_height=7, fig_width=5) {
+  g_deparsed <- paste0(deparse(
+    function() {g}
+  ), collapse = '')
+  
+  sub_chunk <- paste0("```{r sub_chunk_", floor(runif(1) * 10000), ", fig.height=", fig_height, ", fig.width=", fig_width, ", echo=FALSE}",
+                      "\n(",
+                      g_deparsed
+                      , ")()",
+                      "\n```
+  ")
+  
+  cat(knitr::knit(text = knitr::knit_expand(text = sub_chunk), quiet = TRUE))
+}
+
+cum_ci_plot <- function(pivotName){
+  
+  # pivotName = "2.1"
+  
+  # Check if pivotName exists in the data
+  if (!pivotName %in% unique(CI_quadrant$Field)) {
+    # message("PivotName '", pivotName, "' not found. Plotting empty graph.")
+    g <- ggplot() +   labs(title = "Empty Graph - Yield dates missing")
+    return(
+      subchunkify(g, fig_height=2, fig_width = 10)
+    )
+  } else {
+    # message("PivotName '", pivotName, "'  found. Plotting normal graph.")
+    data_ci <- CI_quadrant %>% filter(Field %in% pivotName)
+    
+    
+    
+    data_ci2 <- data_ci  %>% mutate(CI_rate = cumulative_CI/DOY,
+                                    week = week(Date))%>% group_by(subField) %>%
+      mutate(mean_rolling10 = rollapplyr(CI_rate  , width = 10, FUN = mean, partial = TRUE))
+    
+    # date_preperation_perfect_pivot <- data_ci2 %>% group_by(season)  %>% summarise(min_date = min(Date),
+    #                                                                                max_date = max(Date),
+    #                                                                                days = max_date - min_date)
+    
+    # Identify unique seasons
+    filtered_data <- data_ci2 %>%
+      group_by(season) %>%
+      mutate(rank = dense_rank(desc(season))) %>%
+      filter(rank <= 2) %>%
+      ungroup() %>%
+      dplyr::select(-rank)
+    
+    
+    # g <- ggplot(data= data_ci2 %>% filter(season %in% unique_seasons)) +
+    g <- ggplot(data= filtered_data ) +
+      # geom_line(aes(Date, mean_rolling10, col = subField)) +
+      geom_line(aes(Date, CI_rate, col = subField)) +
+      facet_wrap(~season, scales = "free_x") +
+      # geom_line(data= perfect_pivot, aes(Date , mean_rolling10, col = "Model CI (p5.1 Data 2022, \n date x axis is fictive)"), lty="11",size=1) +
+      labs(title = paste("CI rate - Field", pivotName),
+           y = "CI rate (cumulative CI / Age)")+
+      # scale_y_continuous(limits=c(0.5,3), breaks = seq(0.5, 3, 0.5))+
+      scale_x_date(date_breaks = "1 month", date_labels =  "%m-%Y") +
+      theme(axis.text.x = element_text(angle = 60, hjust = 1),
+            legend.justification=c(1,0), legend.position = c(1, 0),
+            legend.title = element_text(size = 8),
+            legend.text = element_text(size = 8)) +
+      guides(color = guide_legend(nrow = 2, byrow = TRUE))
+    subchunkify(g, fig_height=6, fig_width = 10)
+  }
+}
+
 
-library(CIprep)
 
 # Vang alle command line argumenten op
 args <- commandArgs(trailingOnly = TRUE)
@@ -56,7 +495,7 @@ weekly_CI_mosaic <- here(laravel_storage_dir, "weekly_mosaic")
 daily_vrt <- here(data_dir, "vrt")
 harvest_dir <-  here(data_dir, "HarvestData")
 
-source(here("r_app", "parameters_project.R"))
+source("parameters_project.R")
 
 dir.create(here(laravel_storage_dir))
 dir.create(here(data_dir))
@@ -69,19 +508,19 @@ dir.create(merged_final)
 dir.create(harvest_dir)
 
 
-
+weeks_ago = 0
 # Creating weekly mosaic
 dates <- date_list(weeks_ago)
-
+print(dates)
 #load pivot geojson
 # pivot_sf_q <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% vect()
 
 raster_files <- list.files(planet_tif_folder,full.names = T, pattern = ".tif")
-head(raster_files)
 
 filtered_files  <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>%
   compact() %>%
   flatten_chr()
+head(filtered_files)
 
 # filtered_files <- raster_files #for first CI extraction