[200~Improve CI report visualization and migrate to terra package
- Replace raster package with terra throughout the codebase - Update map visualizations with better layout and legends - Add descriptive headers to report sections - Improve map legend positioning and sizing - Enhance error handling for missing data - Remove redundant legends in field-specific visualizations - Optimize figure dimensions to prevent page overflow - Expand documentation of CI index and report components - Update package dependencies in packages.
This commit is contained in:
parent
fbc8f5e11b
commit
07aee7bed1
|
|
@ -17,7 +17,9 @@
|
|||
class ProjectDownloadRDSJob implements ShouldQueue
|
||||
{
|
||||
use Batchable, Dispatchable, InteractsWithQueue, Queueable, SerializesModels;
|
||||
protected Project $project;
|
||||
|
||||
public $timeout = 800;
|
||||
protected Project $project;
|
||||
protected Carbon $date;
|
||||
protected int $offset;
|
||||
|
||||
|
|
@ -41,7 +43,7 @@ public function handle(): void
|
|||
];
|
||||
|
||||
$process = new Process($command);
|
||||
$process->setTimeout(600);
|
||||
$process->setTimeout(800);
|
||||
$process->run();
|
||||
|
||||
if (!$process->isSuccessful()) {
|
||||
|
|
|
|||
|
|
@ -94,23 +94,12 @@ public static function handleFor(Project $project, Carbon $endDate, int $offset)
|
|||
if (Carbon::parse($project->mail_day)->dayOfWeek < $endDate->dayOfWeek) {
|
||||
$endDate->next($project->mail_day);
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* @var ProjectMosaic $mosaic
|
||||
*/
|
||||
$mosaic = $project->mosaics()->updateOrCreate(
|
||||
[
|
||||
'name' => sprintf('Week_%s_%s', $endDate->week, $endDate->year),
|
||||
],
|
||||
[
|
||||
'name' => sprintf('Week_%s_%s', $endDate->week, $endDate->year),
|
||||
'path' => sprintf('%s/%s/%s',
|
||||
$project->download_path,
|
||||
'mosaics',
|
||||
sprintf('week_%s_%s.tif', $endDate->week, $endDate->year)
|
||||
),
|
||||
'end_date' => $endDate->format('Y-m-d'),
|
||||
'offset' => $offset,
|
||||
]);
|
||||
$mosaic = $project->createOrUpdateMosaic($endDate, $offset);
|
||||
return new self($mosaic);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -35,7 +35,7 @@ public function saveMosaic()
|
|||
]);
|
||||
|
||||
$mosaic = $this->project->mosaics()->updateOrCreate([
|
||||
'name' => sprintf('Week %s, %s', $this->formData['week'], $this->formData['year']),
|
||||
'name' => sprintf('Week %s, %s', str_pad($this->formData['week'], 2, '0', STR_PAD_LEFT), $this->formData['year']),
|
||||
'year' => $this->formData['year'],
|
||||
'week' => $this->formData['week'],
|
||||
], [
|
||||
|
|
|
|||
|
|
@ -6,6 +6,7 @@
|
|||
use App\Models\ProjectEmailRecipient;
|
||||
use App\Rules\HarvestFile;
|
||||
use Illuminate\Support\Facades\Validator;
|
||||
use Illuminate\Support\Str;
|
||||
use Illuminate\Validation\Rule;
|
||||
use Livewire\Component;
|
||||
use Livewire\WithFileUploads;
|
||||
|
|
@ -87,11 +88,19 @@ public function createProject()
|
|||
])->validate();
|
||||
$project = Project::create([
|
||||
'name' => $this->formData['name'],
|
||||
'download_path' => $this->formData['name']
|
||||
'download_path' => $this->makeValidDirectoryName($this->formData['name'])
|
||||
]);
|
||||
return redirect()->route('project.show',[$project->name,'settings']);
|
||||
}
|
||||
|
||||
private function makeValidDirectoryName($string)
|
||||
{
|
||||
return Str::of($string)
|
||||
->replaceMatches('/[^a-zA-Z0-9_-]/', '_') // Replace invalid characters
|
||||
->trim() // Remove leading/trailing spaces
|
||||
->lower(); // Convert to lowercase
|
||||
}
|
||||
|
||||
public function saveProject()
|
||||
{
|
||||
$this->resetErrorBag();
|
||||
|
|
|
|||
|
|
@ -565,4 +565,24 @@ public function createAllMosaicsInDatabaseAndUpdateStatusForAllDownloadedImages(
|
|||
})
|
||||
->unique();
|
||||
}
|
||||
|
||||
public function createOrUpdateMosaic(Carbon $date, $offset)
|
||||
{
|
||||
$lpadWeek = str_pad($date->isoWeek(), 2, '0', STR_PAD_LEFT);
|
||||
$name = sprintf('Week_%s_%s', $lpadWeek, $date->isoWeekYear());
|
||||
return $this->mosaics()->updateOrCreate(
|
||||
[
|
||||
'name' => $name,
|
||||
],
|
||||
[
|
||||
'name' => $name,
|
||||
'path' => sprintf('%s/%s/%s',
|
||||
$this->download_path,
|
||||
'mosaics',
|
||||
sprintf('week_%s_%s.tif', $lpadWeek, $date->isoWeekYear())
|
||||
),
|
||||
'end_date' => $date->format('Y-m-d'),
|
||||
'offset' => $offset,
|
||||
]);
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -17,353 +17,497 @@ editor_options:
|
|||
chunk_output_type: console
|
||||
---
|
||||
|
||||
```{r setup, include=FALSE}
|
||||
#set de filename van de output
|
||||
# knitr::opts_chunk$set(echo = TRUE)
|
||||
# output_file <- params$output_file
|
||||
```{r setup_parameters, include=FALSE}
|
||||
# Set up basic report parameters from input values
|
||||
report_date <- params$report_date
|
||||
mail_day <- params$mail_day
|
||||
borders <- params$borders
|
||||
|
||||
|
||||
borders <- params$borders
|
||||
#
|
||||
#
|
||||
# Environment setup notes (commented out)
|
||||
# # Activeer de renv omgeving
|
||||
# renv::activate()
|
||||
# renv::deactivate()
|
||||
# Optioneel: Herstel de omgeving als dat nodig is
|
||||
# Je kunt dit commentaar geven als je het normaal niet wilt uitvoeren
|
||||
# # Optioneel: Herstel de omgeving als dat nodig is
|
||||
# # Je kunt dit commentaar geven als je het normaal niet wilt uitvoeren
|
||||
# renv::restore()
|
||||
```
|
||||
|
||||
```{r libraries, message=FALSE, warning=FALSE, include=FALSE}
|
||||
```{r load_libraries, message=FALSE, warning=FALSE, include=FALSE}
|
||||
# Configure knitr options
|
||||
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
|
||||
|
||||
# Path management
|
||||
library(here)
|
||||
|
||||
# Spatial data libraries
|
||||
library(sf)
|
||||
library(tidyverse)
|
||||
library(terra)
|
||||
library(exactextractr)
|
||||
# library(raster) - Removed as it's no longer maintained
|
||||
|
||||
# Data manipulation and visualization
|
||||
library(tidyverse) # Includes dplyr, ggplot2, etc.
|
||||
library(tmap)
|
||||
library(lubridate)
|
||||
library(exactextractr)
|
||||
library(zoo)
|
||||
library(raster)
|
||||
library(terra)
|
||||
|
||||
# Machine learning
|
||||
library(rsample)
|
||||
library(caret)
|
||||
library(randomForest)
|
||||
library(CAST)
|
||||
|
||||
source("report_utils.R")
|
||||
# source(here("r_app", "report_utils.R"))
|
||||
|
||||
# Load custom utility functions
|
||||
tryCatch({
|
||||
source("report_utils.R")
|
||||
}, error = function(e) {
|
||||
message(paste("Error loading report_utils.R:", e$message))
|
||||
# Try alternative path if the first one fails
|
||||
tryCatch({
|
||||
source(here::here("r_app", "report_utils.R"))
|
||||
}, error = function(e) {
|
||||
stop("Could not load report_utils.R from either location: ", e$message)
|
||||
})
|
||||
})
|
||||
```
|
||||
|
||||
```{r directories, message=FALSE, warning=FALSE, include=FALSE}
|
||||
```{r initialize_project_config, message=FALSE, warning=FALSE, include=FALSE}
|
||||
# Set the project directory from parameters
|
||||
project_dir <- params$data_dir
|
||||
source(here("r_app", "parameters_project.R"))
|
||||
|
||||
log_message("Starting the R Markdown script")
|
||||
log_message(paste("mail_day params:", params$mail_day))
|
||||
log_message(paste("report_date params:", params$report_date))
|
||||
log_message(paste("mail_day variable:", mail_day))
|
||||
# Source project parameters with error handling
|
||||
tryCatch({
|
||||
source(here::here("r_app", "parameters_project.R"))
|
||||
}, error = function(e) {
|
||||
stop("Error loading parameters_project.R: ", e$message)
|
||||
})
|
||||
|
||||
# s2_dir <- "C:/Users/timon/Resilience BV/4002 CMD App - General/4002 CMD Team/4002 TechnicalData/04 WP2 technical/python/chemba_S2/"
|
||||
# Log initial configuration
|
||||
safe_log("Starting the R Markdown script")
|
||||
safe_log(paste("mail_day params:", params$mail_day))
|
||||
safe_log(paste("report_date params:", params$report_date))
|
||||
safe_log(paste("mail_day variable:", mail_day))
|
||||
```
|
||||
|
||||
|
||||
```{r week, message=FALSE, warning=FALSE, include=FALSE}
|
||||
```{r calculate_dates_and_weeks, message=FALSE, warning=FALSE, include=FALSE}
|
||||
# Set locale for consistent date formatting
|
||||
Sys.setlocale("LC_TIME", "C")
|
||||
|
||||
# Initialize date variables from parameters
|
||||
today <- as.character(report_date)
|
||||
mail_day_as_character <- as.character(mail_day)
|
||||
|
||||
report_date_as_week_day <- weekdays(ymd(today))
|
||||
|
||||
# Calculate week days
|
||||
report_date_as_week_day <- weekdays(lubridate::ymd(today))
|
||||
days_of_week <- c("Sunday", "Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday")
|
||||
#als de index of report_date_as_week_day groter dan de index van de mail_day dan moet de week + 1
|
||||
week <- week(today)
|
||||
log_message(paste("week", week, "today", today))
|
||||
today_minus_1 <- as.character(ymd(today) - 7)
|
||||
today_minus_2 <- as.character(ymd(today) - 14)
|
||||
today_minus_3 <- as.character(ymd(today) - 21)
|
||||
|
||||
log_message(paste("report_date_as_week_day", report_date_as_week_day))
|
||||
log_message(paste("which(days_of_week == report_date_as_week_day)", which(days_of_week == report_date_as_week_day)))
|
||||
log_message(paste("mail_day_as_character", mail_day_as_character))
|
||||
log_message(paste(" which(days_of_week == mail_day_as_character)", which(days_of_week == mail_day_as_character)))
|
||||
# Calculate initial week number
|
||||
week <- lubridate::week(today)
|
||||
safe_log(paste("Initial week calculation:", week, "today:", today))
|
||||
|
||||
if (which(days_of_week == report_date_as_week_day) > which(days_of_week == mail_day_as_character)){
|
||||
log_message("adjusting weeks because of mail day")
|
||||
week <- week(today) + 1
|
||||
today_minus_1 <- as.character(ymd(today))
|
||||
today_minus_2 <- as.character(ymd(today) - 7)
|
||||
today_minus_3 <- as.character(ymd(today) - 14)
|
||||
}
|
||||
# Calculate previous dates for comparisons
|
||||
today_minus_1 <- as.character(lubridate::ymd(today) - 7)
|
||||
today_minus_2 <- as.character(lubridate::ymd(today) - 14)
|
||||
today_minus_3 <- as.character(lubridate::ymd(today) - 21)
|
||||
|
||||
# week <- week(today)
|
||||
# Log the weekday calculations for debugging
|
||||
safe_log(paste("Report date weekday:", report_date_as_week_day))
|
||||
safe_log(paste("Weekday index:", which(days_of_week == report_date_as_week_day)))
|
||||
safe_log(paste("Mail day:", mail_day_as_character))
|
||||
safe_log(paste("Mail day index:", which(days_of_week == mail_day_as_character)))
|
||||
|
||||
# week <- 25
|
||||
# today = "2024-06-21"
|
||||
# Adjust week calculation based on mail day
|
||||
if (which(days_of_week == report_date_as_week_day) > which(days_of_week == mail_day_as_character)) {
|
||||
safe_log("Adjusting weeks because of mail day")
|
||||
week <- lubridate::week(today) + 1
|
||||
today_minus_1 <- as.character(lubridate::ymd(today))
|
||||
today_minus_2 <- as.character(lubridate::ymd(today) - 7)
|
||||
today_minus_3 <- as.character(lubridate::ymd(today) - 14)
|
||||
}
|
||||
|
||||
|
||||
#today = as.character(Sys.Date())
|
||||
#week = lubridate::week(Sys.time())
|
||||
## week = 26
|
||||
#title_var <- paste0("CI dashboard week ", week, " - all pivots dashboard using 3x3 meter resolution")
|
||||
# Generate subtitle for report
|
||||
subtitle_var <- paste("Report generated on", Sys.Date())
|
||||
|
||||
# Calculate week numbers for previous weeks
|
||||
week_minus_1 <- week - 1
|
||||
week_minus_2 <- week - 2
|
||||
week_minus_3 <- week - 3
|
||||
|
||||
week_minus_1 <- week -1 # sprintf("%02d", week(today_minus_1))
|
||||
week_minus_2 <- week -2 # sprintf("%02d", week(today_minus_2))
|
||||
week_minus_3 <- week -3 # sprintf("%02d", week(today_minus_3))
|
||||
# Format current week with leading zeros
|
||||
week <- sprintf("%02d", week)
|
||||
|
||||
|
||||
year = year(today)
|
||||
year_1 = year(today_minus_1)
|
||||
year_2 = year(today_minus_2)
|
||||
year_3 = year(today_minus_3)
|
||||
|
||||
# Get years for each date
|
||||
year <- lubridate::year(today)
|
||||
year_1 <- lubridate::year(today_minus_1)
|
||||
year_2 <- lubridate::year(today_minus_2)
|
||||
year_3 <- lubridate::year(today_minus_3)
|
||||
```
|
||||
|
||||
`r subtitle_var`
|
||||
|
||||
\pagebreak
|
||||
# Explanation of the maps
|
||||
# Explanation of the Report
|
||||
|
||||
This PDF-dashboard shows the status of your fields on a weekly basis. We will show this in different ways:
|
||||
This report provides a detailed analysis of your sugarcane fields based on satellite imagery, helping you monitor crop health and development throughout the growing season. The data is processed weekly to give you timely insights for optimal farm management decisions.
|
||||
|
||||
1) The first way is with a general overview of field heterogeneity using ‘variation’ – a higher number indicates more differences between plants in the same field.
|
||||
2) The second map shows a normal image of the latest week in colour, of the demo fields.
|
||||
3) Then come the maps per field, which show the status for three weeks ago, two weeks ago, last week, and this week, as well as a percentage difference map between last week and this week. The percentage difference maps shows the relative difference in growth over the last week, with positive numbers showing growth, and negative numbers showing decline.
|
||||
4) Below the maps are graphs that show how each pivot quadrant is doing, measured through the chlorophyll index.
|
||||
## What is the Chlorophyll Index (CI)?
|
||||
|
||||
```{r data, message=TRUE, warning=TRUE, include=FALSE}
|
||||
# get latest CI index
|
||||
The **Chlorophyll Index (CI)** is a vegetation index that measures the relative amount of chlorophyll in plant leaves. Chlorophyll is the green pigment responsible for photosynthesis in plants. Higher CI values indicate:
|
||||
|
||||
# remove_pivots <- c("1.1", "1.12", "1.8", "1.9", "1.11", "1.14")
|
||||
CI_quadrant <- readRDS(here(cumulative_CI_vals_dir,"All_pivots_Cumulative_CI_quadrant_year_v2.rds"))# %>%
|
||||
# rename(pivot_quadrant = field)
|
||||
* Greater photosynthetic activity
|
||||
* Healthier plant tissue
|
||||
* Better nitrogen uptake
|
||||
* More vigorous crop growth
|
||||
|
||||
CI values typically range from 0 (bare soil or severely stressed vegetation) to 7+ (very healthy, dense vegetation). For sugarcane, values between 3-7 generally indicate good crop health, depending on the growth stage.
|
||||
|
||||
path_to_week_current = here(weekly_CI_mosaic, paste0("week_",week, "_", year, ".tif"))
|
||||
path_to_week_minus_1 = here(weekly_CI_mosaic, paste0("week_",week_minus_1, "_", year_1, ".tif"))
|
||||
path_to_week_minus_2 = here(weekly_CI_mosaic, paste0("week_",week_minus_2, "_", year_2, ".tif"))
|
||||
path_to_week_minus_3 = here(weekly_CI_mosaic, paste0("week_",week_minus_3, "_", year_3, ".tif"))
|
||||
## What You'll Find in This Report:
|
||||
|
||||
log_message("required mosaic paths")
|
||||
log_message(paste("path to week current",path_to_week_current))
|
||||
log_message(paste("path to week minus 1",path_to_week_minus_1))
|
||||
log_message(paste("path to week minus 2",path_to_week_minus_2))
|
||||
log_message(paste("path to week minus 3",path_to_week_minus_3))
|
||||
1. **Chlorophyll Index Overview Map**: A comprehensive view of all your fields showing current CI values. This helps identify which fields are performing well and which might need attention.
|
||||
|
||||
CI <- brick(path_to_week_current) %>% subset("CI")
|
||||
CI_m1 <- brick(path_to_week_minus_1) %>% subset("CI")
|
||||
CI_m2 <- brick(path_to_week_minus_2) %>% subset("CI")
|
||||
CI_m3 <- brick(path_to_week_minus_3) %>% subset("CI")
|
||||
2. **Weekly Difference Map**: Shows changes in CI values over the past week. Positive values (green) indicate improving crop health, while negative values (red) may signal stress or decline.
|
||||
|
||||
3. **Field-by-Field Analysis**: Detailed maps for each field showing:
|
||||
* CI values for the current week and two previous weeks
|
||||
* Week-to-week changes in CI values
|
||||
* Three-week change in CI values to track longer-term trends
|
||||
|
||||
# last_week_dif_raster <- ((CI - CI_m1) / CI_m1) * 100
|
||||
last_week_dif_raster_abs <- (CI - CI_m1)
|
||||
```
|
||||
```{r data_129, message=TRUE, warning=TRUE, include=FALSE}
|
||||
three_week_dif_raster_abs <- (CI - CI_m3)
|
||||
```
|
||||
```{r data_132, message=TRUE, warning=TRUE, include=FALSE}
|
||||
4. **Growth Trend Graphs**: Time-series visualizations showing how CI values have changed throughout the growing season for each section of your fields.
|
||||
|
||||
# AllPivots0 <-st_read(here(data_dir_project, "pivot.geojson"))
|
||||
5. **Yield Prediction**: For mature crops (over 300 days), we provide estimated yield predictions based on historical data and current CI measurements.
|
||||
|
||||
# AllPivots0$pivot <- factor(AllPivots0$pivot, levels = c("1.1", "1.2", "1.3", "1.4", "1.6", "1.7", "1.8", "1.9", "1.10", "1.11", "1.12", "1.13", "1.14" , "1.16" , "1.17" , "1.18" ,"2.1", "2.2", "2.3" , "2.4", "2.5", "3.1", "3.2", "3.3", "4.1", "4.2", "4.3", "4.4", "4.5", "4.6", "5.1" ,"5.2", "5.3", "5.4", "6.1", "6.2", "DL1.1", "DL1.3"))
|
||||
AllPivots0 <- field_boundaries_sf
|
||||
Use these insights to identify areas that may need irrigation, fertilization, or other interventions, and to track the effectiveness of your management practices over time.
|
||||
|
||||
# joined_spans <-st_read(here(data_dir_project, "span.geojson")) %>% st_transform(crs(AllPivots0))
|
||||
|
||||
# pivots_dates <- readRDS(here(harvest_dir, "harvest_data_new"))
|
||||
# pivots_dates$pivot <- factor(pivots_dates$pivot, levels = c("1.1", "1.2", "1.3", "1.4", "1.6", "1.7", "1.8", "1.9", "1.10", "1.11", "1.12", "1.13", "1.14" , "1.16" , "1.17" , "1.18" ,"2.1", "2.2", "2.3" , "2.4", "2.5", "3.1", "3.2", "3.3", "4.1", "4.2", "4.3", "4.4", "4.5", "4.6", "5.1" ,"5.2", "5.3", "5.4", "6.1", "6.2", "DL1.1", "DL1.3"))
|
||||
|
||||
#AllPivots <- merge(AllPivots0, harvesting_data, by = c("field", "sub_field")) #%>%
|
||||
#rename(field = pivot, sub_field = pivot_quadrant) #%>% select(-pivot.y)
|
||||
#head(AllPivots)
|
||||
|
||||
#AllPivots_merged <- AllPivots %>% #dplyr::select(field, sub_field, sub_area) %>% unique() %>%
|
||||
# group_by(field) %>% summarise(sub_area = first(sub_area))
|
||||
|
||||
#AllPivots_merged <- st_transform(AllPivots_merged, crs = proj4string(CI))
|
||||
|
||||
#pivot_names <- unique(CI_quadrant$field)
|
||||
|
||||
```
|
||||
|
||||
|
||||
```{r eval=FALSE, fig.height=7.2, fig.width=10, message=FALSE, warning=FALSE, include=FALSE}
|
||||
RGB_raster <- list.files(paste0(s2_dir,week),full.names = T, pattern = ".tiff", recursive = TRUE)[1] #use pattern = '.tif$' or something else if you have multiple files in this folder
|
||||
|
||||
|
||||
RGB_raster <- brick(RGB_raster)
|
||||
# RGB_raster <- brick(here("planet", paste0("week_",week, ".tif")))
|
||||
tm_shape(RGB_raster, unit = "m") + tm_rgb(r=1, g=2, b=3, max.value = 255) +
|
||||
tm_layout(main.title = paste0("RGB image of the fields - week ", week),
|
||||
main.title.position = 'center') +
|
||||
tm_scale_bar(position = c("right", "top"), text.color = "white") +
|
||||
|
||||
tm_compass(position = c("right", "top"), text.color = "white") +
|
||||
tm_shape(AllPivots0)+ tm_borders( col = "white") +
|
||||
tm_text("pivot_quadrant", size = .6, col = "white")
|
||||
|
||||
```
|
||||
\newpage
|
||||
|
||||
```{r ci_overzicht_kaart, echo=FALSE, fig.height=7.3, fig.width=9, message=FALSE, warning=FALSE}
|
||||
tm_shape(CI, unit = "m")+
|
||||
tm_raster(breaks = c(0,0.5,1,2,3,4,5,6,7,Inf), palette = "RdYlGn", midpoint = NA,legend.is.portrait = F) +
|
||||
tm_layout(legend.outside = TRUE,legend.outside.position = "bottom",legend.show = T, main.title = "Overview all fields (CI)")+
|
||||
tm_scale_bar(position = c("right", "top"), text.color = "black") +
|
||||
|
||||
tm_compass(position = c("right", "top"), text.color = "black") +
|
||||
tm_shape(AllPivots0)+ tm_borders( col = "black") +
|
||||
tm_text("sub_field", size = .6, col = "black")
|
||||
|
||||
```
|
||||
\newpage
|
||||
|
||||
```{r ci_diff_kaart, echo=FALSE, fig.height=7.3, fig.width=9, message=FALSE, warning=FALSE}
|
||||
|
||||
tm_shape(last_week_dif_raster_abs, unit = "m")+
|
||||
tm_raster(breaks = c(-3,-2,-1,0,1,2, 3), palette = "RdYlGn", midpoint = NA,legend.is.portrait = F) +
|
||||
tm_layout(legend.outside = TRUE,legend.outside.position = "bottom",legend.show = T, main.title = "Overview all fields - CI difference")+
|
||||
tm_scale_bar(position = c("right", "top"), text.color = "black") +
|
||||
|
||||
tm_compass(position = c("right", "top"), text.color = "black") +
|
||||
tm_shape(AllPivots0)+ tm_borders( col = "black") +
|
||||
tm_text("sub_field", size = .6, col = "black")
|
||||
|
||||
```
|
||||
|
||||
\newpage
|
||||
```{r plots_ci_estate, eval=TRUE, fig.height=3.8, fig.width=10, message=FALSE,echo=FALSE, warning=FALSE, include=TRUE, results='asis'}
|
||||
AllPivots_merged <- AllPivots0 %>% dplyr::group_by(field) %>% summarise()
|
||||
|
||||
walk(AllPivots_merged$field, ~ {
|
||||
cat("\n") # Add an empty line for better spacing
|
||||
ci_plot(.x)
|
||||
cat("\n")
|
||||
cum_ci_plot(.x)
|
||||
\pagebreak
|
||||
```{r load_ci_quadrant_data, message=TRUE, warning=TRUE, include=FALSE}
|
||||
# Load CI index data with error handling
|
||||
tryCatch({
|
||||
CI_quadrant <- readRDS(here::here(cumulative_CI_vals_dir, "All_pivots_Cumulative_CI_quadrant_year_v2.rds"))
|
||||
safe_log("Successfully loaded CI quadrant data")
|
||||
}, error = function(e) {
|
||||
stop("Error loading CI quadrant data: ", e$message)
|
||||
})
|
||||
```
|
||||
|
||||
```{r looping_over_sub_area, echo=FALSE, fig.height=3.8, fig.width=10, message=FALSE, warning=FALSE, results='asis', eval=FALSE}
|
||||
pivots_grouped <- AllPivots0
|
||||
```{r load_raster_mosaics, message=TRUE, warning=TRUE, include=FALSE}
|
||||
# Get file paths for different weeks using the utility function
|
||||
tryCatch({
|
||||
path_to_week_current = get_week_path(weekly_CI_mosaic, today, 0)
|
||||
path_to_week_minus_1 = get_week_path(weekly_CI_mosaic, today, -1)
|
||||
path_to_week_minus_2 = get_week_path(weekly_CI_mosaic, today, -2)
|
||||
path_to_week_minus_3 = get_week_path(weekly_CI_mosaic, today, -3)
|
||||
|
||||
# Iterate over each subgroup
|
||||
for (subgroup in unique(pivots_grouped$sub_area)) {
|
||||
cat("# HELLO!!!")
|
||||
print(" PRINT")
|
||||
# cat("\n")
|
||||
# cat("# Subgroup: ", subgroup, "\n") # Add a title for the subgroup
|
||||
subset_data <- filter(pivots_grouped, sub_area == subgroup)
|
||||
# cat("\\pagebreak")
|
||||
walk(subset_data$field, ~ {
|
||||
# cat("\n") # Add an empty line for better spacing
|
||||
ci_plot(.x)
|
||||
# cat("\n")
|
||||
cum_ci_plot(.x)
|
||||
# cat("\n")
|
||||
# Log the calculated paths
|
||||
safe_log("Required mosaic paths:")
|
||||
safe_log(paste("Path to current week:", path_to_week_current))
|
||||
safe_log(paste("Path to week minus 1:", path_to_week_minus_1))
|
||||
safe_log(paste("Path to week minus 2:", path_to_week_minus_2))
|
||||
safe_log(paste("Path to week minus 3:", path_to_week_minus_3))
|
||||
|
||||
# Validate that files exist
|
||||
if (!file.exists(path_to_week_current)) warning("Current week mosaic file does not exist: ", path_to_week_current)
|
||||
if (!file.exists(path_to_week_minus_1)) warning("Week minus 1 mosaic file does not exist: ", path_to_week_minus_1)
|
||||
if (!file.exists(path_to_week_minus_2)) warning("Week minus 2 mosaic file does not exist: ", path_to_week_minus_2)
|
||||
if (!file.exists(path_to_week_minus_3)) warning("Week minus 3 mosaic file does not exist: ", path_to_week_minus_3)
|
||||
|
||||
# Load raster data with terra functions
|
||||
CI <- terra::rast(path_to_week_current)$CI
|
||||
CI_m1 <- terra::rast(path_to_week_minus_1)$CI
|
||||
CI_m2 <- terra::rast(path_to_week_minus_2)$CI
|
||||
CI_m3 <- terra::rast(path_to_week_minus_3)$CI
|
||||
|
||||
}, error = function(e) {
|
||||
stop("Error loading raster data: ", e$message)
|
||||
})
|
||||
```
|
||||
|
||||
```{r calculate_difference_rasters, message=TRUE, warning=TRUE, include=FALSE}
|
||||
# Calculate difference rasters for comparisons
|
||||
tryCatch({
|
||||
# Calculate weekly difference
|
||||
last_week_dif_raster_abs <- (CI - CI_m1)
|
||||
safe_log("Calculated weekly difference raster")
|
||||
|
||||
# Calculate three-week difference
|
||||
three_week_dif_raster_abs <- (CI - CI_m3)
|
||||
safe_log("Calculated three-week difference raster")
|
||||
}, error = function(e) {
|
||||
safe_log(paste("Error calculating difference rasters:", e$message), "ERROR")
|
||||
# Create placeholder rasters if calculations fail
|
||||
if (!exists("last_week_dif_raster_abs")) {
|
||||
last_week_dif_raster_abs <- CI * 0
|
||||
}
|
||||
if (!exists("three_week_dif_raster_abs")) {
|
||||
three_week_dif_raster_abs <- CI * 0
|
||||
}
|
||||
})
|
||||
```
|
||||
|
||||
```{r load_field_boundaries, message=TRUE, warning=TRUE, include=FALSE}
|
||||
# Load field boundaries from parameters
|
||||
tryCatch({
|
||||
AllPivots0 <- field_boundaries_sf
|
||||
safe_log("Successfully loaded field boundaries")
|
||||
}, error = function(e) {
|
||||
stop("Error loading field boundaries: ", e$message)
|
||||
})
|
||||
```
|
||||
|
||||
# Chlorophyll Index (CI) Overview Map - Current Week
|
||||
```{r render_ci_overview_map, echo=FALSE, fig.height=6.8, fig.width=9, message=FALSE, warning=FALSE}
|
||||
# Create overview chlorophyll index map
|
||||
tryCatch({
|
||||
tmap::tm_shape(CI, unit = "m") +
|
||||
tmap::tm_raster(breaks = c(0,0.5,1,2,3,4,5,6,7,Inf),
|
||||
palette = "RdYlGn",
|
||||
midpoint = NA,
|
||||
legend.is.portrait = FALSE,
|
||||
title = "Chlorophyll Index (CI)") +
|
||||
tmap::tm_layout(legend.outside = TRUE,
|
||||
legend.outside.position = "bottom",
|
||||
legend.show = TRUE,
|
||||
scale.position = c("right", "bottom"),
|
||||
compass.position = c("right", "bottom")) +
|
||||
tmap::tm_scale_bar(position = "outside", text.color = "black") +
|
||||
tmap::tm_compass(position = "outside", text.color = "black") +
|
||||
tmap::tm_shape(AllPivots0) +
|
||||
tmap::tm_borders(col = "black") +
|
||||
tmap::tm_text("sub_field", size = 0.6, col = "black")
|
||||
}, error = function(e) {
|
||||
safe_log(paste("Error creating CI overview map:", e$message), "ERROR")
|
||||
plot(1, type="n", axes=FALSE, xlab="", ylab="")
|
||||
text(1, 1, "Error creating CI overview map", cex=1.5)
|
||||
})
|
||||
```
|
||||
\newpage
|
||||
|
||||
# Weekly Chlorophyll Index Difference Map
|
||||
```{r render_ci_difference_map, echo=FALSE, fig.height=6.8, fig.width=9, message=FALSE, warning=FALSE}
|
||||
# Create chlorophyll index difference map
|
||||
tryCatch({
|
||||
tmap::tm_shape(last_week_dif_raster_abs, unit = "m") +
|
||||
tmap::tm_raster(breaks = c(-3,-2,-1,0,1,2,3),
|
||||
palette = "RdYlGn",
|
||||
midpoint = NA,
|
||||
legend.is.portrait = FALSE,
|
||||
title = "Chlorophyll Index (CI) Change") +
|
||||
tmap::tm_layout(legend.outside = TRUE,
|
||||
legend.outside.position = "bottom",
|
||||
legend.show = TRUE,
|
||||
scale.position = c("right", "bottom"),
|
||||
compass.position = c("right", "bottom")) +
|
||||
tmap::tm_scale_bar(position = "outside", text.color = "black") +
|
||||
tmap::tm_compass(position = "outside", text.color = "black") +
|
||||
tmap::tm_shape(AllPivots0) +
|
||||
tmap::tm_borders(col = "black") +
|
||||
tmap::tm_text("sub_field", size = 0.6, col = "black")
|
||||
}, error = function(e) {
|
||||
safe_log(paste("Error creating CI difference map:", e$message), "ERROR")
|
||||
plot(1, type="n", axes=FALSE, xlab="", ylab="")
|
||||
text(1, 1, "Error creating CI difference map", cex=1.5)
|
||||
})
|
||||
```
|
||||
\newpage
|
||||
|
||||
```{r generate_field_visualizations, eval=TRUE, fig.height=3.8, fig.width=10, message=FALSE,echo=FALSE, warning=FALSE, include=TRUE, results='asis'}
|
||||
# Generate detailed visualizations for each field
|
||||
tryCatch({
|
||||
# Merge field polygons for processing
|
||||
AllPivots_merged <- AllPivots0 %>%
|
||||
dplyr::group_by(field) %>%
|
||||
dplyr::summarise(.groups = 'drop')
|
||||
|
||||
# Generate plots for each field
|
||||
purrr::walk(AllPivots_merged$field, function(field_name) {
|
||||
tryCatch({
|
||||
cat("\n") # Add an empty line for better spacing
|
||||
ci_plot(field_name)
|
||||
cat("\n")
|
||||
cum_ci_plot(field_name)
|
||||
}, error = function(e) {
|
||||
safe_log(paste("Error generating plots for field", field_name, ":", e$message), "ERROR")
|
||||
cat(paste("## Error generating plots for field", field_name, "\n"))
|
||||
cat(paste(e$message, "\n"))
|
||||
})
|
||||
})
|
||||
}
|
||||
}, error = function(e) {
|
||||
safe_log(paste("Error in field visualization section:", e$message), "ERROR")
|
||||
cat("Error generating field plots. See log for details.\n")
|
||||
})
|
||||
```
|
||||
|
||||
```{r generate_subarea_visualizations, echo=FALSE, fig.height=3.8, fig.width=10, message=FALSE, warning=FALSE, results='asis', eval=FALSE}
|
||||
# Alternative visualization grouped by sub-area (disabled by default)
|
||||
tryCatch({
|
||||
# Group pivots by sub-area
|
||||
pivots_grouped <- AllPivots0
|
||||
|
||||
# Iterate over each subgroup
|
||||
for (subgroup in unique(pivots_grouped$sub_area)) {
|
||||
# Add subgroup heading
|
||||
cat("\n")
|
||||
cat("## Subgroup: ", subgroup, "\n")
|
||||
|
||||
# Filter data for current subgroup
|
||||
subset_data <- dplyr::filter(pivots_grouped, sub_area == subgroup)
|
||||
|
||||
# Generate visualizations for each field in the subgroup
|
||||
purrr::walk(subset_data$field, function(field_name) {
|
||||
cat("\n")
|
||||
ci_plot(field_name)
|
||||
cat("\n")
|
||||
cum_ci_plot(field_name)
|
||||
cat("\n")
|
||||
})
|
||||
|
||||
# Add page break after each subgroup
|
||||
cat("\\pagebreak\n")
|
||||
}
|
||||
}, error = function(e) {
|
||||
safe_log(paste("Error in subarea visualization section:", e$message), "ERROR")
|
||||
cat("Error generating subarea plots. See log for details.\n")
|
||||
})
|
||||
```
|
||||
|
||||
# Yield prediction
|
||||
The below table shows estimates of the biomass if you would harvest them now.
|
||||
|
||||
```{r message=FALSE, warning=FALSE, include=FALSE}
|
||||
|
||||
#CI_quadrant <- readRDS(here(cumulative_CI_vals_dir,"All_pivots_Cumulative_CI_quadrant_year_v2.rds"))
|
||||
CI_quadrant <- readRDS(here(cumulative_CI_vals_dir,"All_pivots_Cumulative_CI_quadrant_year_v2.rds")) %>%
|
||||
group_by(model) %>% # Group by model or other relevant columns
|
||||
fill(field, sub_field, .direction = "downup") %>% # Fill down then up within each group
|
||||
ungroup()
|
||||
|
||||
# Check if tonnage_ha is empty
|
||||
if (all(is.na(CI_quadrant$tonnage_ha))) {
|
||||
log_message("Lacking historic harvest data, please provide for yield prediction calculation")
|
||||
knitr::knit_exit() # Exit the chunk if tonnage_ha is empty
|
||||
}
|
||||
|
||||
harvesting_data <- harvesting_data %>% rename(season = year)
|
||||
|
||||
CI_and_yield <- left_join(CI_quadrant , harvesting_data, by = c("field", "sub_field", "season")) %>% #filter(!is.na(tonnage_ha)) %>%
|
||||
group_by(sub_field, season) %>% slice(which.max(DOY)) %>%
|
||||
dplyr::select(field, sub_field, tonnage_ha, cumulative_CI, DOY, season, sub_area) %>%
|
||||
mutate(CI_per_day = cumulative_CI/ DOY)
|
||||
|
||||
predictors <- c( "cumulative_CI" , "DOY" ,"CI_per_day" )
|
||||
response <- "tonnage_ha"
|
||||
# CI_and_yield_test <- as.data.frame(CI_and_yield_test)
|
||||
CI_and_yield_test <- CI_and_yield %>% as.data.frame() %>% filter(!is.na(tonnage_ha))
|
||||
CI_and_yield_validation <- CI_and_yield_test
|
||||
prediction_yields <- CI_and_yield %>% as.data.frame() %>% filter(is.na(tonnage_ha))
|
||||
|
||||
ctrl <- trainControl(method="cv",
|
||||
savePredictions = TRUE,
|
||||
allowParallel= TRUE,
|
||||
number = 5,
|
||||
verboseIter = TRUE)
|
||||
|
||||
set.seed(202)
|
||||
model_ffs_rf <- ffs( CI_and_yield_test[,predictors],
|
||||
CI_and_yield_test[,response],
|
||||
method="rf" ,
|
||||
trControl=ctrl,
|
||||
importance=TRUE,
|
||||
withinSE = TRUE,
|
||||
tuneLength = 5,
|
||||
na.rm = TRUE
|
||||
)
|
||||
|
||||
|
||||
# Function to prepare predictions
|
||||
prepare_predictions <- function(predictions, newdata) {
|
||||
return(predictions %>%
|
||||
as.data.frame() %>%
|
||||
rename(predicted_Tcha = ".") %>%
|
||||
mutate(sub_field = newdata$sub_field,
|
||||
field = newdata$field,
|
||||
Age_days = newdata$DOY,
|
||||
total_CI = round(newdata$cumulative_CI, 0),
|
||||
predicted_Tcha = round(predicted_Tcha, 0),
|
||||
season = newdata$season) %>%
|
||||
dplyr::select(field, sub_field, Age_days, total_CI, predicted_Tcha, season) %>%
|
||||
left_join(., newdata, by = c("field", "sub_field", "season")))
|
||||
}
|
||||
# Predict yields for the validation dataset
|
||||
pred_ffs_rf <- prepare_predictions(predict(model_ffs_rf, newdata = CI_and_yield_validation), CI_and_yield_validation)
|
||||
|
||||
# Predict yields for the current season
|
||||
pred_rf_current_season <- prepare_predictions(predict(model_ffs_rf, newdata = prediction_yields), prediction_yields) %>%
|
||||
filter(Age_days > 300) %>%
|
||||
mutate(CI_per_day = round(total_CI / Age_days, 1))
|
||||
```{r yield_data_training, message=FALSE, warning=FALSE, include=FALSE}
|
||||
# Load and prepare yield prediction data with error handling
|
||||
tryCatch({
|
||||
# Load CI quadrant data and fill missing values
|
||||
CI_quadrant <- readRDS(here::here(cumulative_CI_vals_dir, "All_pivots_Cumulative_CI_quadrant_year_v2.rds")) %>%
|
||||
dplyr::group_by(model) %>%
|
||||
tidyr::fill(field, sub_field, .direction = "downup") %>%
|
||||
dplyr::ungroup()
|
||||
|
||||
# Check if tonnage_ha is empty
|
||||
if (all(is.na(CI_quadrant$tonnage_ha))) {
|
||||
safe_log("Lacking historic harvest data, please provide for yield prediction calculation", "WARNING")
|
||||
knitr::knit_exit() # Exit the chunk if tonnage_ha is empty
|
||||
}
|
||||
|
||||
# Rename year column to season for consistency
|
||||
harvesting_data <- harvesting_data %>% dplyr::rename(season = year)
|
||||
|
||||
# Join CI and yield data
|
||||
CI_and_yield <- dplyr::left_join(CI_quadrant, harvesting_data, by = c("field", "sub_field", "season")) %>%
|
||||
dplyr::group_by(sub_field, season) %>%
|
||||
dplyr::slice(which.max(DOY)) %>%
|
||||
dplyr::select(field, sub_field, tonnage_ha, cumulative_CI, DOY, season, sub_area) %>%
|
||||
dplyr::mutate(CI_per_day = cumulative_CI / DOY)
|
||||
|
||||
# Define predictors and response variables
|
||||
predictors <- c("cumulative_CI", "DOY", "CI_per_day")
|
||||
response <- "tonnage_ha"
|
||||
|
||||
# Prepare test and validation datasets
|
||||
CI_and_yield_test <- CI_and_yield %>%
|
||||
as.data.frame() %>%
|
||||
dplyr::filter(!is.na(tonnage_ha))
|
||||
|
||||
CI_and_yield_validation <- CI_and_yield_test
|
||||
|
||||
# Prepare prediction dataset (fields without harvest data)
|
||||
prediction_yields <- CI_and_yield %>%
|
||||
as.data.frame() %>%
|
||||
dplyr::filter(is.na(tonnage_ha))
|
||||
|
||||
# Configure model training parameters
|
||||
ctrl <- caret::trainControl(
|
||||
method = "cv",
|
||||
savePredictions = TRUE,
|
||||
allowParallel = TRUE,
|
||||
number = 5,
|
||||
verboseIter = TRUE
|
||||
)
|
||||
|
||||
# Train the model with feature selection
|
||||
set.seed(202) # For reproducibility
|
||||
model_ffs_rf <- CAST::ffs(
|
||||
CI_and_yield_test[, predictors],
|
||||
CI_and_yield_test[, response],
|
||||
method = "rf",
|
||||
trControl = ctrl,
|
||||
importance = TRUE,
|
||||
withinSE = TRUE,
|
||||
tuneLength = 5,
|
||||
na.rm = TRUE
|
||||
)
|
||||
|
||||
# Function to prepare predictions with consistent naming and formatting
|
||||
prepare_predictions <- function(predictions, newdata) {
|
||||
return(predictions %>%
|
||||
as.data.frame() %>%
|
||||
dplyr::rename(predicted_Tcha = ".") %>%
|
||||
dplyr::mutate(
|
||||
sub_field = newdata$sub_field,
|
||||
field = newdata$field,
|
||||
Age_days = newdata$DOY,
|
||||
total_CI = round(newdata$cumulative_CI, 0),
|
||||
predicted_Tcha = round(predicted_Tcha, 0),
|
||||
season = newdata$season
|
||||
) %>%
|
||||
dplyr::select(field, sub_field, Age_days, total_CI, predicted_Tcha, season) %>%
|
||||
dplyr::left_join(., newdata, by = c("field", "sub_field", "season"))
|
||||
)
|
||||
}
|
||||
|
||||
# Predict yields for the validation dataset
|
||||
pred_ffs_rf <- prepare_predictions(stats::predict(model_ffs_rf, newdata = CI_and_yield_validation), CI_and_yield_validation)
|
||||
|
||||
# Predict yields for the current season (focus on mature fields over 300 days)
|
||||
pred_rf_current_season <- prepare_predictions(stats::predict(model_ffs_rf, newdata = prediction_yields), prediction_yields) %>%
|
||||
dplyr::filter(Age_days > 300) %>%
|
||||
dplyr::mutate(CI_per_day = round(total_CI / Age_days, 1))
|
||||
|
||||
safe_log("Successfully completed yield prediction calculations")
|
||||
|
||||
}, error = function(e) {
|
||||
safe_log(paste("Error in yield prediction:", e$message), "ERROR")
|
||||
# Create empty dataframes to prevent errors in subsequent chunks
|
||||
pred_ffs_rf <- data.frame()
|
||||
pred_rf_current_season <- data.frame()
|
||||
})
|
||||
```
|
||||
|
||||
```{r yield_plaatjes, echo=FALSE }
|
||||
ggplot(pred_ffs_rf, aes(y = predicted_Tcha, x = tonnage_ha)) +
|
||||
geom_point(size = 2, alpha = 0.6) + # Adjust point size and transparency
|
||||
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "red") + # Reference line
|
||||
scale_x_continuous(limits = c(0, 200)) +
|
||||
scale_y_continuous(limits = c(0, 200)) +
|
||||
labs(title = "Model Performance: \nPredicted vs Actual Tonnage/ha",
|
||||
x = "Actual tonnage/ha (Tcha)",
|
||||
y = "Predicted tonnage/ha (Tcha)") +
|
||||
theme_minimal()
|
||||
|
||||
ggplot(pred_rf_current_season, aes(x = Age_days, y = predicted_Tcha)) +
|
||||
geom_point(size = 2, alpha = 0.6) + # Adjust point size and transparency
|
||||
labs(title = "Predicted Yields for Fields Over 300 Days \nOld Yet to Be Harvested",
|
||||
x = "Age (days)",
|
||||
y = "Predicted tonnage/ha (Tcha)") +
|
||||
# facet_wrap(~sub_area) +
|
||||
scale_y_continuous(limits = c(0, 200)) + # Optional: Set limits for y-axis
|
||||
theme_minimal()
|
||||
|
||||
knitr::kable(pred_rf_current_season,
|
||||
digits = 0,
|
||||
caption = "Predicted Tonnage/ha for Fields Over 300 Days Old")
|
||||
```{r plotting_yield_data, echo=FALSE, fig.height=5, fig.width=8, message=FALSE, warning=FALSE}
|
||||
# Display yield prediction visualizations with error handling
|
||||
tryCatch({
|
||||
if (nrow(pred_ffs_rf) > 0) {
|
||||
# Plot model performance (predicted vs actual)
|
||||
ggplot2::ggplot(pred_ffs_rf, ggplot2::aes(y = predicted_Tcha, x = tonnage_ha)) +
|
||||
ggplot2::geom_point(size = 2, alpha = 0.6) +
|
||||
ggplot2::geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "red") +
|
||||
ggplot2::scale_x_continuous(limits = c(0, 200)) +
|
||||
ggplot2::scale_y_continuous(limits = c(0, 200)) +
|
||||
ggplot2::labs(title = "Model Performance: \nPredicted vs Actual Tonnage/ha",
|
||||
x = "Actual tonnage/ha (Tcha)",
|
||||
y = "Predicted tonnage/ha (Tcha)") +
|
||||
ggplot2::theme_minimal()
|
||||
}
|
||||
|
||||
if (nrow(pred_rf_current_season) > 0) {
|
||||
# Plot predicted yields by age
|
||||
ggplot2::ggplot(pred_rf_current_season, ggplot2::aes(x = Age_days, y = predicted_Tcha)) +
|
||||
ggplot2::geom_point(size = 2, alpha = 0.6) +
|
||||
ggplot2::labs(title = "Predicted Yields for Fields Over 300 Days \nOld Yet to Be Harvested",
|
||||
x = "Age (days)",
|
||||
y = "Predicted tonnage/ha (Tcha)") +
|
||||
ggplot2::scale_y_continuous(limits = c(0, 200)) +
|
||||
ggplot2::theme_minimal()
|
||||
|
||||
# Display prediction table
|
||||
knitr::kable(pred_rf_current_season,
|
||||
digits = 0,
|
||||
caption = "Predicted Tonnage/ha for Fields Over 300 Days Old")
|
||||
} else {
|
||||
cat("No fields over 300 days old without harvest data available for yield prediction.")
|
||||
}
|
||||
}, error = function(e) {
|
||||
safe_log(paste("Error in yield prediction visualization:", e$message), "ERROR")
|
||||
cat("Error generating yield prediction visualizations. See log for details.")
|
||||
})
|
||||
```
|
||||
````````````````
|
||||
|
||||
|
|
|
|||
BIN
r_app/Rplots.pdf
BIN
r_app/Rplots.pdf
Binary file not shown.
|
|
@ -1,137 +1,100 @@
|
|||
# nolint start: commented_code_linter, line_length_linter,object_usage_linter.
|
||||
library(sf)
|
||||
library(terra)
|
||||
library(tidyverse)
|
||||
library(lubridate)
|
||||
library(exactextractr)
|
||||
library(readxl)
|
||||
# CI_EXTRACTION.R
|
||||
# ==============
|
||||
# This script processes satellite imagery to extract Canopy Index (CI) values for agricultural fields.
|
||||
# It handles image processing, masking, and extraction of statistics by field/sub-field.
|
||||
#
|
||||
# Usage: Rscript ci_extraction.R [end_date] [offset] [project_dir]
|
||||
# - end_date: End date for processing (YYYY-MM-DD format)
|
||||
# - offset: Number of days to look back from end_date
|
||||
# - project_dir: Project directory name (e.g., "chemba")
|
||||
#
|
||||
|
||||
# Vang alle command line argumenten op
|
||||
args <- commandArgs(trailingOnly = TRUE)
|
||||
# 1. Load required packages
|
||||
# -----------------------
|
||||
suppressPackageStartupMessages({
|
||||
library(sf)
|
||||
library(terra)
|
||||
library(tidyverse)
|
||||
library(lubridate)
|
||||
library(exactextractr)
|
||||
library(readxl)
|
||||
library(here)
|
||||
})
|
||||
|
||||
# Controleer of er ten minste één argument is doorgegeven
|
||||
if (length(args) == 0) {
|
||||
stop("Geen argumenten doorgegeven aan het script")
|
||||
}
|
||||
|
||||
# Converteer het eerste argument naar een numerieke waarde
|
||||
end_date <- as.Date(args[1])
|
||||
if (is.na(end_date)) {
|
||||
end_date <- lubridate::dmy("28-08-2024")
|
||||
}
|
||||
|
||||
offset <- as.numeric(args[2])
|
||||
# Controleer of weeks_ago een geldig getal is
|
||||
if (is.na(offset)) {
|
||||
# stop("Het argument is geen geldig getal")
|
||||
offset <- 7
|
||||
}
|
||||
|
||||
week <- week(end_date)
|
||||
# Converteer het tweede argument naar een string waarde
|
||||
project_dir <- as.character(args[3])
|
||||
|
||||
# Controleer of data_dir een geldige waarde is
|
||||
if (!is.character(project_dir)) {
|
||||
project_dir <- "chemba"
|
||||
}
|
||||
new_project_question = FALSE
|
||||
|
||||
source("parameters_project.R")
|
||||
source("ci_extraction_utils.R")
|
||||
|
||||
dates <- date_list(end_date, offset)
|
||||
print(dates)
|
||||
|
||||
raster_files <- list.files(planet_tif_folder,full.names = T, pattern = ".tif")
|
||||
|
||||
filtered_files <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>%
|
||||
compact() %>%
|
||||
flatten_chr()
|
||||
|
||||
# Remove files that do not exist
|
||||
existing_files <- filtered_files[file.exists(filtered_files)]
|
||||
|
||||
# Check if the list of existing files is empty
|
||||
if (length(existing_files) == 0) {
|
||||
message("No files exist for the given date(s).")
|
||||
stop("Terminating script.")
|
||||
# 2. Process command line arguments
|
||||
# ------------------------------
|
||||
main <- function() {
|
||||
# Capture command line arguments
|
||||
args <- commandArgs(trailingOnly = TRUE)
|
||||
|
||||
# Process end_date argument
|
||||
if (length(args) >= 1 && !is.na(args[1])) {
|
||||
end_date <- as.Date(args[1])
|
||||
if (is.na(end_date)) {
|
||||
warning("Invalid end_date provided. Using default (current date).")
|
||||
#end_date <- Sys.Date()
|
||||
end_date <- "2023-10-01"
|
||||
}
|
||||
} else {
|
||||
#end_date <- Sys.Date()
|
||||
end_date <- "2023-10-01"
|
||||
}
|
||||
|
||||
# Continue with the rest of the script
|
||||
print(existing_files)
|
||||
|
||||
|
||||
vrt_list <- list()
|
||||
|
||||
for (file in existing_files) {
|
||||
v_crop <- create_mask_and_crop(file, field_boundaries, merged_final)
|
||||
emtpy_or_full <- global(v_crop, "notNA")
|
||||
|
||||
vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt"))
|
||||
if(emtpy_or_full[1,] > 100){
|
||||
vrt_list[vrt_file] <- vrt_file
|
||||
|
||||
}else{
|
||||
file.remove(vrt_file)
|
||||
|
||||
# Process offset argument
|
||||
if (length(args) >= 2 && !is.na(args[2])) {
|
||||
offset <- as.numeric(args[2])
|
||||
if (is.na(offset) || offset <= 0) {
|
||||
warning("Invalid offset provided. Using default (7 days).")
|
||||
offset <- 7
|
||||
}
|
||||
} else {
|
||||
offset <- 7
|
||||
}
|
||||
|
||||
message(file, " processed")
|
||||
gc()
|
||||
|
||||
# Process project_dir argument
|
||||
if (length(args) >= 3 && !is.na(args[3])) {
|
||||
project_dir <- as.character(args[3])
|
||||
} else {
|
||||
project_dir <- "chemba"
|
||||
}
|
||||
|
||||
# 3. Initialize project configuration
|
||||
# --------------------------------
|
||||
new_project_question <- FALSE
|
||||
source("r_app/parameters_project.R")
|
||||
source("r_app/ci_extraction_utils.R")
|
||||
|
||||
# 4. Generate date list for processing
|
||||
# ---------------------------------
|
||||
dates <- date_list(end_date, offset)
|
||||
log_message(paste("Processing data for week", dates$week, "of", dates$year))
|
||||
|
||||
# 5. Find and filter raster files by date
|
||||
# -----------------------------------
|
||||
log_message("Searching for raster files")
|
||||
|
||||
tryCatch({
|
||||
# Use the new utility function to find satellite images
|
||||
existing_files <- find_satellite_images(planet_tif_folder, dates$days_filter)
|
||||
log_message(paste("Found", length(existing_files), "raster files for processing"))
|
||||
|
||||
# 6. Process raster files and create VRT
|
||||
# -----------------------------------
|
||||
# Use the new utility function for batch processing
|
||||
vrt_list <- process_satellite_images(existing_files, field_boundaries, merged_final, daily_vrt)
|
||||
|
||||
# 7. Process and combine CI values
|
||||
# ------------------------------
|
||||
# Call the process_ci_values function from utils with all required parameters
|
||||
process_ci_values(dates, field_boundaries, merged_final,
|
||||
field_boundaries_sf, daily_CI_vals_dir, cumulative_CI_vals_dir)
|
||||
|
||||
}, error = function(e) {
|
||||
log_message(paste("Error in main processing:", e$message), level = "ERROR")
|
||||
stop(e$message)
|
||||
})
|
||||
}
|
||||
|
||||
raster_files_NEW <- list.files(merged_final,full.names = T, pattern = ".tif")
|
||||
# Run the main function if the script is executed directly
|
||||
main()
|
||||
|
||||
# Define the path to the file
|
||||
file_path <- here(cumulative_CI_vals_dir, "combined_CI_data.rds")
|
||||
|
||||
# Check if the file exists
|
||||
if (!file.exists(file_path)) {
|
||||
# File does not exist, create it with all available data
|
||||
print("combined_CI_data.rds does not exist. Preparing combined_CI_data.rds file for all available images.")
|
||||
|
||||
# Extract data from all raster files
|
||||
walk(raster_files_NEW, extract_rasters_daily, field_geojson = field_boundaries, quadrants = FALSE, daily_CI_vals_dir)
|
||||
|
||||
# Combine all extracted data
|
||||
extracted_values <- list.files(here(daily_CI_vals_dir), full.names = TRUE)
|
||||
|
||||
pivot_stats <- extracted_values %>%
|
||||
map(readRDS) %>% list_rbind() %>%
|
||||
group_by(sub_field)
|
||||
|
||||
# Save the combined data to the file
|
||||
saveRDS(pivot_stats, file_path)
|
||||
|
||||
print("All CI values extracted from all historic images and saved to combined_CI_data.rds.")
|
||||
|
||||
} else {
|
||||
# File exists, add new data
|
||||
print("combined_CI_data.rds exists, adding the latest image data to the table.")
|
||||
|
||||
# Filter and process the latest data
|
||||
filtered_files <- map(dates$days_filter, ~ raster_files_NEW[grepl(pattern = .x, x = raster_files_NEW)]) %>%
|
||||
compact() %>%
|
||||
flatten_chr()
|
||||
|
||||
walk(filtered_files, extract_rasters_daily, field_geojson = field_boundaries, quadrants = TRUE, daily_CI_vals_dir)
|
||||
|
||||
# Extract new values
|
||||
extracted_values <- list.files(daily_CI_vals_dir, full.names = TRUE)
|
||||
extracted_values <- map(dates$days_filter, ~ extracted_values[grepl(pattern = .x, x = extracted_values)]) %>%
|
||||
compact() %>%
|
||||
flatten_chr()
|
||||
|
||||
pivot_stats <- extracted_values %>%
|
||||
map(readRDS) %>% list_rbind() %>%
|
||||
group_by(sub_field)
|
||||
|
||||
# Load existing data and append new data
|
||||
combined_CI_data <- readRDS(file_path)
|
||||
pivot_stats2 <- bind_rows(pivot_stats, combined_CI_data)
|
||||
|
||||
# Save the updated combined data
|
||||
saveRDS(pivot_stats2, file_path)
|
||||
|
||||
print("All CI values extracted from the latest images and added to combined_CI_data.rds.")
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,53 +1,426 @@
|
|||
# Utils for ci extraction
|
||||
date_list <- function(end_date, offset){
|
||||
offset <- as.numeric(offset) - 1
|
||||
end_date <- as.Date(end_date)
|
||||
# CI_EXTRACTION_UTILS.R
|
||||
# =====================
|
||||
# Utility functions for the SmartCane CI (Chlorophill Index) extraction workflow.
|
||||
# These functions support date handling, raster processing, and data extraction.
|
||||
|
||||
#' Safe logging function that works whether log_message exists or not
|
||||
#'
|
||||
#' @param message The message to log
|
||||
#' @param level The log level (default: "INFO")
|
||||
#' @return NULL (used for side effects)
|
||||
#'
|
||||
safe_log <- function(message, level = "INFO") {
|
||||
if (exists("log_message")) {
|
||||
log_message(message, level)
|
||||
} else {
|
||||
if (level %in% c("ERROR", "WARNING")) {
|
||||
warning(message)
|
||||
} else {
|
||||
message(message)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#' Generate a sequence of dates for processing
|
||||
#'
|
||||
#' @param end_date The end date for the sequence (Date object)
|
||||
#' @param offset Number of days to look back from end_date
|
||||
#' @return A list containing week number, year, and a sequence of dates for filtering
|
||||
#'
|
||||
date_list <- function(end_date, offset) {
|
||||
# Input validation
|
||||
if (!lubridate::is.Date(end_date)) {
|
||||
end_date <- as.Date(end_date)
|
||||
if (is.na(end_date)) {
|
||||
stop("Invalid end_date provided. Expected a Date object or a string convertible to Date.")
|
||||
}
|
||||
}
|
||||
|
||||
offset <- as.numeric(offset)
|
||||
if (is.na(offset) || offset < 1) {
|
||||
stop("Invalid offset provided. Expected a positive number.")
|
||||
}
|
||||
|
||||
# Calculate date range
|
||||
offset <- offset - 1 # Adjust offset to include end_date
|
||||
start_date <- end_date - lubridate::days(offset)
|
||||
|
||||
week <- week(start_date)
|
||||
year <- year(start_date)
|
||||
# Extract week and year information
|
||||
week <- lubridate::week(start_date)
|
||||
year <- lubridate::year(start_date)
|
||||
|
||||
# Generate sequence of dates
|
||||
days_filter <- seq(from = start_date, to = end_date, by = "day")
|
||||
days_filter <- format(days_filter, "%Y-%m-%d") # Format for consistent filtering
|
||||
|
||||
return(list("week" = week, "year" = year, "days_filter" = days_filter))
|
||||
# Log the date range
|
||||
safe_log(paste("Date range generated from", start_date, "to", end_date))
|
||||
|
||||
return(list(
|
||||
"week" = week,
|
||||
"year" = year,
|
||||
"days_filter" = days_filter,
|
||||
"start_date" = start_date,
|
||||
"end_date" = end_date
|
||||
))
|
||||
}
|
||||
|
||||
#' Create a Chlorophill Index (CI) mask from satellite imagery and crop to field boundaries
|
||||
#'
|
||||
#' @param file Path to the satellite image file
|
||||
#' @param field_boundaries Field boundaries vector object
|
||||
#' @param merged_final_dir Directory to save the processed raster
|
||||
#' @return Processed raster object with CI band
|
||||
#'
|
||||
create_mask_and_crop <- function(file, field_boundaries, merged_final_dir) {
|
||||
message("starting ", file)
|
||||
loaded_raster <- rast(file)
|
||||
names(loaded_raster) <- c("Red", "Green", "Blue", "NIR")
|
||||
message("raster loaded")
|
||||
# Validate inputs
|
||||
if (!file.exists(file)) {
|
||||
stop(paste("File not found:", file))
|
||||
}
|
||||
|
||||
CI <- loaded_raster$NIR / loaded_raster$Green - 1
|
||||
if (is.null(field_boundaries)) {
|
||||
stop("Field boundaries are required but were not provided")
|
||||
}
|
||||
|
||||
loaded_raster$CI <- CI
|
||||
message("CI calculated")
|
||||
loaded_raster <- terra::crop(loaded_raster, field_boundaries, mask = TRUE) #%>% CI_func()
|
||||
# Establish file names for output
|
||||
basename_no_ext <- tools::file_path_sans_ext(basename(file))
|
||||
new_file <- here::here(merged_final_dir, paste0(basename_no_ext, ".tif"))
|
||||
vrt_file <- here::here(daily_vrt, paste0(basename_no_ext, ".vrt"))
|
||||
|
||||
loaded_raster[loaded_raster == 0] <- NA
|
||||
new_file <- here(merged_final_dir, 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)
|
||||
|
||||
return(loaded_raster)
|
||||
# Process with error handling
|
||||
tryCatch({
|
||||
# Log processing start
|
||||
safe_log(paste("Processing", basename(file)))
|
||||
|
||||
# Load and prepare raster
|
||||
loaded_raster <- terra::rast(file)
|
||||
|
||||
# Validate raster has necessary bands
|
||||
if (terra::nlyr(loaded_raster) < 4) {
|
||||
stop("Raster must have at least 4 bands (Red, Green, Blue, NIR)")
|
||||
}
|
||||
|
||||
# Name bands for clarity
|
||||
names(loaded_raster) <- c("Red", "Green", "Blue", "NIR")
|
||||
|
||||
# Calculate Canopy Index
|
||||
CI <- loaded_raster$NIR / loaded_raster$Green - 1
|
||||
|
||||
# Add CI to raster and mask invalid values
|
||||
loaded_raster$CI <- CI
|
||||
loaded_raster <- terra::crop(loaded_raster, field_boundaries, mask = TRUE)
|
||||
|
||||
# Replace zeros with NA for better visualization and analysis
|
||||
loaded_raster[loaded_raster == 0] <- NA
|
||||
|
||||
# Write output files
|
||||
terra::writeRaster(loaded_raster, new_file, overwrite = TRUE)
|
||||
terra::vrt(new_file, vrt_file, overwrite = TRUE)
|
||||
|
||||
# Check if the result has enough valid pixels
|
||||
valid_pixels <- terra::global(loaded_raster$CI, "notNA", na.rm=TRUE)
|
||||
|
||||
# Log completion
|
||||
safe_log(paste("Completed processing", basename(file),
|
||||
"- Valid pixels:", valid_pixels[1,]))
|
||||
|
||||
return(loaded_raster)
|
||||
|
||||
}, error = function(e) {
|
||||
err_msg <- paste("Error processing", basename(file), "-", e$message)
|
||||
safe_log(err_msg, "ERROR")
|
||||
return(NULL)
|
||||
}, finally = {
|
||||
# Clean up memory
|
||||
gc()
|
||||
})
|
||||
}
|
||||
|
||||
#' Process a batch of satellite images and create VRT files
|
||||
#'
|
||||
#' @param files Vector of file paths to process
|
||||
#' @param field_boundaries Field boundaries vector object for cropping
|
||||
#' @param merged_final_dir Directory to save processed rasters
|
||||
#' @param daily_vrt_dir Directory to save VRT files
|
||||
#' @param min_valid_pixels Minimum number of valid pixels for a raster to be kept (default: 100)
|
||||
#' @return List of valid VRT files created
|
||||
#'
|
||||
process_satellite_images <- function(files, field_boundaries, merged_final_dir, daily_vrt_dir, min_valid_pixels = 100) {
|
||||
vrt_list <- list()
|
||||
|
||||
safe_log(paste("Starting batch processing of", length(files), "files"))
|
||||
|
||||
# Process each file
|
||||
for (file in files) {
|
||||
# Process each raster file
|
||||
v_crop <- create_mask_and_crop(file, field_boundaries, merged_final_dir)
|
||||
|
||||
# Skip if processing failed
|
||||
if (is.null(v_crop)) {
|
||||
next
|
||||
}
|
||||
|
||||
# Check if the raster has enough valid data
|
||||
valid_data <- terra::global(v_crop, "notNA")
|
||||
vrt_file <- here::here(daily_vrt_dir, paste0(tools::file_path_sans_ext(basename(file)), ".vrt"))
|
||||
|
||||
if (valid_data[1,] > min_valid_pixels) {
|
||||
vrt_list[[vrt_file]] <- vrt_file
|
||||
} else {
|
||||
# Remove VRT files with insufficient data
|
||||
if (file.exists(vrt_file)) {
|
||||
file.remove(vrt_file)
|
||||
}
|
||||
safe_log(paste("Skipping", basename(file), "- insufficient valid data"), "WARNING")
|
||||
}
|
||||
|
||||
# Clean up memory
|
||||
rm(v_crop)
|
||||
gc()
|
||||
}
|
||||
|
||||
safe_log(paste("Completed processing", length(vrt_list), "raster files"))
|
||||
|
||||
return(vrt_list)
|
||||
}
|
||||
|
||||
#' Find satellite image files filtered by date
|
||||
#'
|
||||
#' @param tif_folder Directory containing satellite imagery files
|
||||
#' @param dates_filter Character vector of dates in YYYY-MM-DD format
|
||||
#' @return Vector of file paths matching the date filter
|
||||
#'
|
||||
find_satellite_images <- function(tif_folder, dates_filter) {
|
||||
# Find all raster files
|
||||
raster_files <- list.files(tif_folder, full.names = TRUE, pattern = "\\.tif$")
|
||||
|
||||
if (length(raster_files) == 0) {
|
||||
stop("No raster files found in directory: ", tif_folder)
|
||||
}
|
||||
|
||||
# Filter files by dates
|
||||
filtered_files <- purrr::map(dates_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>%
|
||||
purrr::compact() %>%
|
||||
purrr::flatten_chr()
|
||||
|
||||
# Remove files that do not exist
|
||||
existing_files <- filtered_files[file.exists(filtered_files)]
|
||||
|
||||
# Check if the list of existing files is empty
|
||||
if (length(existing_files) == 0) {
|
||||
stop("No files found matching the date filter: ", paste(dates_filter, collapse = ", "))
|
||||
}
|
||||
|
||||
return(existing_files)
|
||||
}
|
||||
|
||||
#' Extract date from file path
|
||||
#'
|
||||
#' @param file_path Path to the file
|
||||
#' @return Extracted date in YYYY-MM-DD format
|
||||
#'
|
||||
date_extract <- function(file_path) {
|
||||
str_extract(file_path, "\\d{4}-\\d{2}-\\d{2}")
|
||||
date <- stringr::str_extract(file_path, "\\d{4}-\\d{2}-\\d{2}")
|
||||
|
||||
if (is.na(date)) {
|
||||
warning(paste("Could not extract date from file path: ", file_path))
|
||||
}
|
||||
|
||||
return(date)
|
||||
}
|
||||
|
||||
#' Extract CI values from a raster for each field or subfield
|
||||
#'
|
||||
#' @param file Path to the raster file
|
||||
#' @param field_geojson Field boundaries as SF object
|
||||
#' @param quadrants Boolean indicating whether to extract by quadrants
|
||||
#' @param save_dir Directory to save the extracted values
|
||||
#' @return Path to the saved RDS file
|
||||
#'
|
||||
extract_rasters_daily <- function(file, field_geojson, quadrants = TRUE, save_dir) {
|
||||
field_geojson <- sf::st_as_sf(field_geojson)
|
||||
x <- rast(file[1])
|
||||
# Validate inputs
|
||||
if (!file.exists(file)) {
|
||||
stop(paste("File not found: ", file))
|
||||
}
|
||||
|
||||
if (!inherits(field_geojson, "sf") && !inherits(field_geojson, "sfc")) {
|
||||
field_geojson <- sf::st_as_sf(field_geojson)
|
||||
}
|
||||
|
||||
# Extract date from file path
|
||||
date <- date_extract(file)
|
||||
if (is.na(date)) {
|
||||
stop(paste("Could not extract date from file path:", 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)
|
||||
# Log extraction start
|
||||
safe_log(paste("Extracting CI values for", date, "- Using quadrants:", quadrants))
|
||||
|
||||
save_suffix <- if (quadrants){"quadrant"} else {"whole_field"}
|
||||
save_path <- here(save_dir, paste0("extracted_", date, "_", save_suffix, ".rds"))
|
||||
|
||||
saveRDS(pivot_stats, save_path)
|
||||
# Process with error handling
|
||||
tryCatch({
|
||||
# Load raster
|
||||
x <- terra::rast(file)
|
||||
|
||||
# Check if CI band exists
|
||||
if (!"CI" %in% names(x)) {
|
||||
stop("CI band not found in raster")
|
||||
}
|
||||
|
||||
# Extract statistics
|
||||
pivot_stats <- cbind(
|
||||
field_geojson,
|
||||
mean_CI = round(exactextractr::exact_extract(x$CI, field_geojson, fun = "mean"), 2)
|
||||
) %>%
|
||||
sf::st_drop_geometry() %>%
|
||||
dplyr::rename("{date}" := mean_CI)
|
||||
|
||||
# Determine save path
|
||||
save_suffix <- if (quadrants) {"quadrant"} else {"whole_field"}
|
||||
save_path <- here::here(save_dir, paste0("extracted_", date, "_", save_suffix, ".rds"))
|
||||
|
||||
# Save extracted data
|
||||
saveRDS(pivot_stats, save_path)
|
||||
|
||||
# Log success
|
||||
safe_log(paste("Successfully extracted and saved CI values for", date))
|
||||
|
||||
return(save_path)
|
||||
|
||||
}, error = function(e) {
|
||||
err_msg <- paste("Error extracting CI values for", date, "-", e$message)
|
||||
safe_log(err_msg, "ERROR")
|
||||
return(NULL)
|
||||
})
|
||||
}
|
||||
|
||||
#' Combine daily CI values into a single dataset
|
||||
#'
|
||||
#' @param daily_CI_vals_dir Directory containing daily CI values
|
||||
#' @param output_file Path to save the combined dataset
|
||||
#' @return Combined dataset as a tibble
|
||||
#'
|
||||
combine_ci_values <- function(daily_CI_vals_dir, output_file = NULL) {
|
||||
# List all RDS files in the daily CI values directory
|
||||
files <- list.files(path = daily_CI_vals_dir, pattern = "^extracted_.*\\.rds$", full.names = TRUE)
|
||||
|
||||
if (length(files) == 0) {
|
||||
stop("No extracted CI values found in directory:", daily_CI_vals_dir)
|
||||
}
|
||||
|
||||
# Log process start
|
||||
safe_log(paste("Combining", length(files), "CI value files"))
|
||||
|
||||
# Load and combine all files
|
||||
combined_data <- files %>%
|
||||
purrr::map(readRDS) %>%
|
||||
purrr::list_rbind() %>%
|
||||
dplyr::group_by(sub_field)
|
||||
|
||||
# Save if output file is specified
|
||||
if (!is.null(output_file)) {
|
||||
saveRDS(combined_data, output_file)
|
||||
safe_log(paste("Combined CI values saved to", output_file))
|
||||
}
|
||||
|
||||
return(combined_data)
|
||||
}
|
||||
|
||||
#' Update existing CI data with new values
|
||||
#'
|
||||
#' @param new_data New CI data to be added
|
||||
#' @param existing_data_file Path to the existing data file
|
||||
#' @return Updated combined dataset
|
||||
#'
|
||||
update_ci_data <- function(new_data, existing_data_file) {
|
||||
if (!file.exists(existing_data_file)) {
|
||||
warning(paste("Existing data file not found:", existing_data_file))
|
||||
return(new_data)
|
||||
}
|
||||
|
||||
# Load existing data
|
||||
existing_data <- readRDS(existing_data_file)
|
||||
|
||||
# Combine data, handling duplicates by keeping the newer values
|
||||
combined_data <- dplyr::bind_rows(new_data, existing_data) %>%
|
||||
dplyr::distinct() %>%
|
||||
dplyr::group_by(sub_field)
|
||||
|
||||
# Save updated data
|
||||
saveRDS(combined_data, existing_data_file)
|
||||
safe_log(paste("Updated CI data saved to", existing_data_file))
|
||||
|
||||
return(combined_data)
|
||||
}
|
||||
|
||||
#' Process and combine CI values from raster files
|
||||
#'
|
||||
#' @param dates List of dates from date_list()
|
||||
#' @param field_boundaries Field boundaries as vector object
|
||||
#' @param merged_final_dir Directory with processed raster files
|
||||
#' @param field_boundaries_sf Field boundaries as SF object
|
||||
#' @param daily_CI_vals_dir Directory to save daily CI values
|
||||
#' @param cumulative_CI_vals_dir Directory to save cumulative CI values
|
||||
#' @return NULL (used for side effects)
|
||||
#'
|
||||
process_ci_values <- function(dates, field_boundaries, merged_final_dir,
|
||||
field_boundaries_sf, daily_CI_vals_dir,
|
||||
cumulative_CI_vals_dir) {
|
||||
# Find processed raster files
|
||||
raster_files <- list.files(merged_final_dir, full.names = TRUE, pattern = "\\.tif$")
|
||||
|
||||
# Define path for combined CI data
|
||||
combined_ci_path <- here::here(cumulative_CI_vals_dir, "combined_CI_data.rds")
|
||||
|
||||
# Check if the combined CI data file exists
|
||||
if (!file.exists(combined_ci_path)) {
|
||||
# Process all available data if file doesn't exist
|
||||
safe_log("combined_CI_data.rds does not exist. Creating new file with all available data.")
|
||||
|
||||
# Extract data from all raster files
|
||||
purrr::walk(
|
||||
raster_files,
|
||||
extract_rasters_daily,
|
||||
field_geojson = field_boundaries_sf,
|
||||
quadrants = FALSE,
|
||||
save_dir = daily_CI_vals_dir
|
||||
)
|
||||
|
||||
# Combine all extracted data
|
||||
pivot_stats <- combine_ci_values(daily_CI_vals_dir, combined_ci_path)
|
||||
safe_log("All CI values extracted from historic images and saved.")
|
||||
|
||||
} else {
|
||||
# Process only the latest data and add to existing file
|
||||
safe_log("combined_CI_data.rds exists, adding the latest image data.")
|
||||
|
||||
# Filter files by dates
|
||||
filtered_files <- purrr::map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>%
|
||||
purrr::compact() %>%
|
||||
purrr::flatten_chr()
|
||||
|
||||
# Extract data for the new files
|
||||
purrr::walk(
|
||||
filtered_files,
|
||||
extract_rasters_daily,
|
||||
field_geojson = field_boundaries_sf,
|
||||
quadrants = TRUE,
|
||||
save_dir = daily_CI_vals_dir
|
||||
)
|
||||
|
||||
# Filter extracted values files by the current date range
|
||||
extracted_values <- list.files(daily_CI_vals_dir, full.names = TRUE)
|
||||
extracted_values <- purrr::map(dates$days_filter, ~ extracted_values[grepl(pattern = .x, x = extracted_values)]) %>%
|
||||
purrr::compact() %>%
|
||||
purrr::flatten_chr()
|
||||
|
||||
# Combine new values
|
||||
new_pivot_stats <- extracted_values %>%
|
||||
purrr::map(readRDS) %>%
|
||||
purrr::list_rbind() %>%
|
||||
dplyr::group_by(sub_field)
|
||||
|
||||
# Update the combined data file
|
||||
update_ci_data(new_pivot_stats, combined_ci_path)
|
||||
safe_log("CI values from latest images added to combined_CI_data.rds")
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,15 +1,12 @@
|
|||
library(raster)
|
||||
library(rgdal)
|
||||
library(ggplot2)
|
||||
library(dplyr)
|
||||
library(sf)
|
||||
library(terra)
|
||||
|
||||
# Define the file paths
|
||||
raster_dir <- "C:/Users/timon/Resilience BV/4020 SCane ESA DEMO - Documenten/General/4020 SCDEMO Team/4020 TechnicalData/WP3/smartcane/laravel_app/storage/app/chemba/merged_final_tif"
|
||||
polygon_path <- "C:/Users/timon/Resilience BV/4020 SCane ESA DEMO - Documenten/General/4020 SCDEMO Team/4020 TechnicalData/WP3/smartcane/r_app/experiments/pivot.geojson"
|
||||
|
||||
# Load the polygon
|
||||
library(sf)
|
||||
|
||||
# Load the polygon
|
||||
polygon <- st_read(polygon_path)
|
||||
|
||||
|
|
@ -17,9 +14,9 @@ polygon <- st_read(polygon_path)
|
|||
polygon_filtered <- subset(polygon, field == "1.1")
|
||||
|
||||
# List all raster files in the directory
|
||||
library(terra)
|
||||
raster_files <- list.files(raster_dir, pattern = "\\.tif$", full.names = FALSE)
|
||||
raster_file <- raster_files[1]
|
||||
|
||||
# Initialize an empty list to store CI values
|
||||
ci_values <- list()
|
||||
|
||||
|
|
@ -29,19 +26,18 @@ process_raster_file <- function(raster_file) {
|
|||
date <- as.Date(substr(raster_file, 1, 10), format = "%Y-%m-%d")
|
||||
|
||||
# Load the raster
|
||||
raster_data <- rast(file.path(raster_dir, raster_file))
|
||||
raster_data <- terra::rast(file.path(raster_dir, raster_file))
|
||||
|
||||
# Crop the raster using the filtered polygon
|
||||
cropped_raster <- crop(raster_data, polygon_filtered)
|
||||
cropped_raster <- terra::crop(raster_data, polygon_filtered)
|
||||
|
||||
# Extract the CI band (assuming the CI band is the first band)
|
||||
ci_band <- cropped_raster$CI
|
||||
|
||||
# Extract the CI values from the CI band
|
||||
return(data.frame(Date = date, CI = values(ci_band)) %>%
|
||||
return(data.frame(Date = date, CI = terra::values(ci_band)) %>%
|
||||
filter(!is.na(CI)) %>%
|
||||
mutate(Polygon = unique(polygon_filtered$field)))
|
||||
|
||||
}
|
||||
|
||||
# Use walk to apply the function to each raster file
|
||||
|
|
@ -53,7 +49,7 @@ ci_values_df <- do.call(rbind, ci_values)
|
|||
head(ci_values_df)
|
||||
|
||||
# Plot the first raster
|
||||
plot(rast(raster_files[1])[[1]])
|
||||
plot(terra::rast(file.path(raster_dir, raster_files[1]))[[1]])
|
||||
|
||||
# Calculate the mean CI value
|
||||
ci_stats <- ci_values_df %>%
|
||||
|
|
@ -65,23 +61,23 @@ ci_stats <- ci_values_df %>%
|
|||
q50 = quantile(CI, 0.50, na.rm = TRUE),
|
||||
q75 = quantile(CI, 0.75, na.rm = TRUE),
|
||||
q95 = quantile(CI, 0.95, na.rm = TRUE),
|
||||
cv = sd(CI, na.rm = TRUE) / mean(CI, na.rm = TRUE) * 100
|
||||
cv = sd(CI, na.rm = TRUE) / mean(CI, na.rm = TRUE) * 100,
|
||||
.groups = "drop"
|
||||
)
|
||||
|
||||
# Plot the coefficient of variation over time
|
||||
ggplot(ci_stats, aes(x = Date, y = cv, group = Polygon)) +
|
||||
geom_line(color = "green") +
|
||||
labs(title = "Coefficient of Variation (CV) Over Time",
|
||||
x = "Date",
|
||||
y = "Coefficient of Variation (%)") +
|
||||
# Plot the mean CI value over time
|
||||
ggplot(ci_stats, aes(x = Date, y = mean_ci)) +
|
||||
geom_line() +
|
||||
geom_ribbon(aes(ymin = q25, ymax = q75), alpha = 0.2) +
|
||||
labs(title = "Mean CI value over time",
|
||||
x = "Date",
|
||||
y = "Mean CI") +
|
||||
theme_minimal()
|
||||
|
||||
# Plot the CI statistics over time with different quantiles and coefficient of variation
|
||||
ggplot(ci_stats, aes(x = Date, group = Polygon)) +
|
||||
geom_line(aes(y = q50), color = "red") +
|
||||
geom_ribbon(aes(ymin = q25, ymax = q75), fill = "blue", alpha = 0.2) +
|
||||
geom_ribbon(aes(ymin = q5, ymax = q95), fill = "blue", alpha = 0.1) +
|
||||
labs(title = "CI Statistics Over Time",
|
||||
x = "Date",
|
||||
y = "CI Value / Coefficient of Variation (%)") +
|
||||
# Plot the coefficient of variation over time
|
||||
ggplot(ci_stats, aes(x = Date, y = cv)) +
|
||||
geom_line() +
|
||||
labs(title = "Coefficient of variation over time",
|
||||
x = "Date",
|
||||
y = "CV (%)") +
|
||||
theme_minimal()
|
||||
|
|
|
|||
233
r_app/growth_model_utils.R
Normal file
233
r_app/growth_model_utils.R
Normal file
|
|
@ -0,0 +1,233 @@
|
|||
# filepath: c:\Users\timon\Resilience BV\4020 SCane ESA DEMO - Documenten\General\4020 SCDEMO Team\4020 TechnicalData\WP3\smartcane\r_app\growth_model_utils.R
|
||||
#
|
||||
# GROWTH_MODEL_UTILS.R
|
||||
# ===================
|
||||
# Utility functions for growth model interpolation and manipulation.
|
||||
# These functions support the creation of continuous growth models from point measurements.
|
||||
|
||||
#' Safe logging function that works whether log_message exists or not
|
||||
#'
|
||||
#' @param message The message to log
|
||||
#' @param level The log level (default: "INFO")
|
||||
#' @return NULL (used for side effects)
|
||||
#'
|
||||
safe_log <- function(message, level = "INFO") {
|
||||
if (exists("log_message")) {
|
||||
log_message(message, level)
|
||||
} else {
|
||||
if (level %in% c("ERROR", "WARNING")) {
|
||||
warning(message)
|
||||
} else {
|
||||
message(message)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#' Load and prepare the combined CI data
|
||||
#'
|
||||
#' @param data_dir Directory containing the combined CI data
|
||||
#' @return Long-format dataframe with CI values by date
|
||||
#'
|
||||
load_combined_ci_data <- function(data_dir) {
|
||||
file_path <- here::here(data_dir, "combined_CI_data.rds")
|
||||
|
||||
if (!file.exists(file_path)) {
|
||||
stop(paste("Combined CI data file not found:", file_path))
|
||||
}
|
||||
|
||||
safe_log(paste("Loading CI data from:", file_path))
|
||||
|
||||
# Load and transform the data to long format
|
||||
pivot_stats <- readRDS(file_path) %>%
|
||||
dplyr::ungroup() %>%
|
||||
dplyr::group_by(field, sub_field) %>%
|
||||
dplyr::summarise(dplyr::across(everything(), ~ first(stats::na.omit(.))), .groups = "drop")
|
||||
|
||||
pivot_stats_long <- pivot_stats %>%
|
||||
tidyr::gather("Date", value, -field, -sub_field) %>%
|
||||
dplyr::mutate(Date = lubridate::ymd(Date)) %>%
|
||||
tidyr::drop_na(c("value", "Date")) %>%
|
||||
dplyr::mutate(value = as.numeric(value)) %>%
|
||||
dplyr::filter_all(dplyr::all_vars(!is.infinite(.))) %>%
|
||||
dplyr::distinct()
|
||||
|
||||
safe_log(paste("Loaded", nrow(pivot_stats_long), "CI data points"))
|
||||
|
||||
return(pivot_stats_long)
|
||||
}
|
||||
|
||||
#' Extract and interpolate CI data for a specific field and season
|
||||
#'
|
||||
#' @param field_name Name of the field or sub-field
|
||||
#' @param harvesting_data Dataframe with harvesting information
|
||||
#' @param field_CI_data Dataframe with CI measurements
|
||||
#' @param season Year of the growing season
|
||||
#' @return Dataframe with interpolated daily CI values
|
||||
#'
|
||||
extract_CI_data <- function(field_name, harvesting_data, field_CI_data, season) {
|
||||
# Filter harvesting data for the given season and field name
|
||||
filtered_harvesting_data <- harvesting_data %>%
|
||||
dplyr::filter(year == season, sub_field == field_name)
|
||||
|
||||
if (nrow(filtered_harvesting_data) == 0) {
|
||||
safe_log(paste("No harvesting data found for field:", field_name, "in season:", season), "WARNING")
|
||||
return(data.frame())
|
||||
}
|
||||
|
||||
# Filter field CI data for the given field name
|
||||
filtered_field_CI_data <- field_CI_data %>%
|
||||
dplyr::filter(sub_field == field_name)
|
||||
|
||||
# Return an empty data frame if no CI data is found
|
||||
if (nrow(filtered_field_CI_data) == 0) {
|
||||
safe_log(paste("No CI data found for field:", field_name, "in season:", season), "WARNING")
|
||||
return(data.frame())
|
||||
}
|
||||
|
||||
# Log season dates
|
||||
season_start <- filtered_harvesting_data$season_start[1]
|
||||
season_end <- filtered_harvesting_data$season_end[1]
|
||||
ci_date_range <- paste(format(min(filtered_field_CI_data$Date), "%Y-%m-%d"),
|
||||
"to",
|
||||
format(max(filtered_field_CI_data$Date), "%Y-%m-%d"))
|
||||
|
||||
# Create a linear interpolation function for the CI data
|
||||
tryCatch({
|
||||
ApproxFun <- stats::approxfun(x = filtered_field_CI_data$Date, y = filtered_field_CI_data$value)
|
||||
Dates <- seq.Date(min(filtered_field_CI_data$Date), max(filtered_field_CI_data$Date), by = 1)
|
||||
LinearFit <- ApproxFun(Dates)
|
||||
|
||||
# Combine interpolated data with the original CI data
|
||||
CI <- data.frame(Date = Dates, FitData = LinearFit) %>%
|
||||
dplyr::left_join(filtered_field_CI_data, by = "Date") %>%
|
||||
dplyr::filter(Date > filtered_harvesting_data$season_start & Date < filtered_harvesting_data$season_end)
|
||||
|
||||
# If CI is empty after filtering, return an empty dataframe
|
||||
if (nrow(CI) == 0) {
|
||||
safe_log(paste0("No CI data within season dates for field: ", field_name,
|
||||
" (Season: ", season, ", dates: ",
|
||||
format(season_start, "%Y-%m-%d"), " to ",
|
||||
format(season_end, "%Y-%m-%d"),
|
||||
"). Available CI data range: ", ci_date_range),
|
||||
"WARNING")
|
||||
return(data.frame())
|
||||
}
|
||||
|
||||
# Add additional columns
|
||||
CI <- CI %>%
|
||||
dplyr::mutate(
|
||||
DOY = seq(1, n(), 1),
|
||||
model = paste0("Data", season, " : ", field_name),
|
||||
season = season,
|
||||
subField = field_name
|
||||
)
|
||||
|
||||
# Log successful interpolation
|
||||
safe_log(paste0("Successfully interpolated CI data for field: ", field_name,
|
||||
" (Season: ", season, ", dates: ",
|
||||
format(season_start, "%Y-%m-%d"), " to ",
|
||||
format(season_end, "%Y-%m-%d"),
|
||||
"). ", nrow(CI), " data points created."))
|
||||
|
||||
return(CI)
|
||||
}, error = function(e) {
|
||||
safe_log(paste0("Error interpolating CI data for field ", field_name,
|
||||
" in season ", season,
|
||||
" (", format(season_start, "%Y-%m-%d"), " to ",
|
||||
format(season_end, "%Y-%m-%d"),
|
||||
"): ", e$message), "ERROR")
|
||||
return(data.frame())
|
||||
})
|
||||
}
|
||||
|
||||
#' Generate interpolated CI data for all fields and seasons
|
||||
#'
|
||||
#' @param years Vector of years to process
|
||||
#' @param harvesting_data Dataframe with harvesting information
|
||||
#' @param ci_data Long-format dataframe with CI measurements
|
||||
#' @return Dataframe with interpolated daily CI values for all fields/seasons
|
||||
#'
|
||||
generate_interpolated_ci_data <- function(years, harvesting_data, ci_data) {
|
||||
safe_log("Starting CI data interpolation for all fields")
|
||||
|
||||
# Process each year
|
||||
result <- purrr::map_df(years, function(yr) {
|
||||
safe_log(paste("Processing year:", yr))
|
||||
|
||||
# Get the fields harvested in this year with valid season start dates
|
||||
sub_fields <- harvesting_data %>%
|
||||
dplyr::filter(year == yr, !is.na(season_start)) %>%
|
||||
dplyr::pull(sub_field)
|
||||
|
||||
if (length(sub_fields) == 0) {
|
||||
safe_log(paste("No fields with valid season data for year:", yr), "WARNING")
|
||||
return(data.frame())
|
||||
}
|
||||
|
||||
# Filter sub_fields to only include those with value data in ci_data
|
||||
valid_sub_fields <- sub_fields %>%
|
||||
purrr::keep(~ any(ci_data$sub_field == .x))
|
||||
|
||||
if (length(valid_sub_fields) == 0) {
|
||||
safe_log(paste("No fields with CI data for year:", yr), "WARNING")
|
||||
return(data.frame())
|
||||
}
|
||||
|
||||
# Extract and interpolate data for each valid field
|
||||
safe_log(paste("Processing", length(valid_sub_fields), "fields for year:", yr))
|
||||
|
||||
result <- purrr::map(valid_sub_fields, ~ extract_CI_data(.x,
|
||||
harvesting_data = harvesting_data,
|
||||
field_CI_data = ci_data,
|
||||
season = yr)) %>%
|
||||
purrr::list_rbind()
|
||||
|
||||
safe_log(paste("Generated", nrow(result), "interpolated data points for year:", yr))
|
||||
return(result)
|
||||
})
|
||||
|
||||
safe_log(paste("Total interpolated data points:", nrow(result)))
|
||||
return(result)
|
||||
}
|
||||
|
||||
#' Calculate growth metrics for interpolated CI data
|
||||
#'
|
||||
#' @param interpolated_data Dataframe with interpolated CI values
|
||||
#' @return Dataframe with added growth metrics (CI_per_day and cumulative_CI)
|
||||
#'
|
||||
calculate_growth_metrics <- function(interpolated_data) {
|
||||
if (nrow(interpolated_data) == 0) {
|
||||
safe_log("No data provided to calculate growth metrics", "WARNING")
|
||||
return(interpolated_data)
|
||||
}
|
||||
|
||||
result <- interpolated_data %>%
|
||||
dplyr::group_by(model) %>%
|
||||
dplyr::mutate(
|
||||
CI_per_day = FitData - dplyr::lag(FitData),
|
||||
cumulative_CI = cumsum(FitData)
|
||||
)
|
||||
|
||||
return(result)
|
||||
}
|
||||
|
||||
#' Save interpolated growth model data
|
||||
#'
|
||||
#' @param data Dataframe with interpolated growth data
|
||||
#' @param output_dir Directory to save the output
|
||||
#' @param file_name Filename for the output (default: "All_pivots_Cumulative_CI_quadrant_year_v2.rds")
|
||||
#' @return Path to the saved file
|
||||
#'
|
||||
save_growth_model <- function(data, output_dir, file_name = "All_pivots_Cumulative_CI_quadrant_year_v2.rds") {
|
||||
# Create output directory if it doesn't exist
|
||||
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
|
||||
|
||||
# Create full file path
|
||||
file_path <- here::here(output_dir, file_name)
|
||||
|
||||
# Save the data
|
||||
saveRDS(data, file_path)
|
||||
|
||||
safe_log(paste("Interpolated CI data saved to:", file_path))
|
||||
return(file_path)
|
||||
}
|
||||
|
|
@ -1,124 +1,86 @@
|
|||
library(sf)
|
||||
library(terra)
|
||||
library(tidyverse)
|
||||
library(lubridate)
|
||||
library(exactextractr)
|
||||
# filepath: c:\Users\timon\Resilience BV\4020 SCane ESA DEMO - Documenten\General\4020 SCDEMO Team\4020 TechnicalData\WP3\smartcane\r_app\interpolate_growth_model.R
|
||||
#
|
||||
# INTERPOLATE_GROWTH_MODEL.R
|
||||
# =========================
|
||||
# This script interpolates CI (Chlorophyll Index) values between measurement dates
|
||||
# to create a continuous growth model. It generates daily values and cumulative
|
||||
# CI statistics for each field.
|
||||
#
|
||||
# Usage: Rscript interpolate_growth_model.R [project_dir]
|
||||
# - project_dir: Project directory name (e.g., "chemba")
|
||||
#
|
||||
|
||||
# Vang alle command line argumenten op
|
||||
args <- commandArgs(trailingOnly = TRUE)
|
||||
|
||||
# Converteer het tweede argument naar een string waarde
|
||||
project_dir <- as.character(args[1])
|
||||
|
||||
# Controleer of data_dir een geldige waarde is
|
||||
if (!is.character(project_dir)) {
|
||||
project_dir <- "chemba"
|
||||
}
|
||||
|
||||
source("parameters_project.R")
|
||||
source("ci_extraction_utils.R")
|
||||
# Check if the file exists
|
||||
file_path <- here(cumulative_CI_vals_dir, "combined_CI_data.rds")
|
||||
pivot_stats2 <- data.frame()
|
||||
if (file.exists(file_path)) {
|
||||
pivot_stats2 <- readRDS(file_path) %>%
|
||||
ungroup() %>%
|
||||
group_by(field, sub_field) %>%
|
||||
summarise(across(everything(), ~ first(na.omit(.))), .groups = "drop")
|
||||
}
|
||||
head(pivot_stats2)
|
||||
#%>% drop_na(pivot_quadrant)
|
||||
|
||||
# gather data into long format for easier calculation and visualisation
|
||||
pivot_stats_long <- pivot_stats2 %>%
|
||||
tidyr::gather("Date", value, -field, -sub_field ) %>%
|
||||
mutate(#Date = right(Date, 8),
|
||||
Date = lubridate::ymd(Date)
|
||||
) %>%
|
||||
drop_na(c("value","Date")) %>%
|
||||
mutate(value = as.numeric(value))%>%
|
||||
filter_all(all_vars(!is.infinite(.))) %>%
|
||||
# rename(field = pivot_quadrant,
|
||||
# sub_field = field) %>%
|
||||
unique()
|
||||
|
||||
years <- harvesting_data %>%
|
||||
filter(!is.na(season_start)) %>%
|
||||
distinct(year) %>%
|
||||
pull(year)
|
||||
|
||||
extract_CI_data <- function(field_names, harvesting_data, field_CI_data, season) {
|
||||
# Filter harvesting data for the given season and field names
|
||||
filtered_harvesting_data <- harvesting_data %>%
|
||||
filter(year == season, sub_field %in% field_names)
|
||||
|
||||
# Filter field CI data for the given field names
|
||||
filtered_field_CI_data <- field_CI_data %>%
|
||||
filter(sub_field %in% field_names)
|
||||
|
||||
# Return an empty data frame if no CI data is found
|
||||
if (nrow(filtered_field_CI_data) == 0) {
|
||||
return(data.frame())
|
||||
}
|
||||
|
||||
|
||||
# Create a linear interpolation function for the CI data
|
||||
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)
|
||||
# Combine interpolated data with the original CI data
|
||||
|
||||
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)
|
||||
|
||||
# If CI is empty after filtering, return an empty dataframe
|
||||
if (nrow(CI) == 0) {
|
||||
message ('CI empty after filtering')
|
||||
return(data.frame())
|
||||
}
|
||||
|
||||
# Add additional columns if data exists
|
||||
CI <- CI %>%
|
||||
mutate(DOY = seq(1, n(), 1),
|
||||
model = paste0("Data", season, " : ", field_names),
|
||||
season = season,
|
||||
subField = field_names)
|
||||
|
||||
return(CI)
|
||||
}
|
||||
|
||||
message(harvesting_data)
|
||||
|
||||
CI_all <- map_df(years, function(yr) {
|
||||
# yr = 2021
|
||||
message(yr)
|
||||
|
||||
# Get the fields harvested in this year
|
||||
sub_fields <- harvesting_data %>%
|
||||
filter(year == yr) %>%
|
||||
filter(!is.na(season_start)) %>%
|
||||
pull(sub_field)
|
||||
|
||||
|
||||
# Filter sub_fields to only include those with value data in pivot_stats_long
|
||||
valid_sub_fields <- sub_fields %>%
|
||||
keep(~ any(pivot_stats_long$sub_field == .x))
|
||||
|
||||
# Extract data for each valid field
|
||||
map(valid_sub_fields, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = yr)) %>%
|
||||
list_rbind()
|
||||
# 1. Load required packages
|
||||
# -----------------------
|
||||
suppressPackageStartupMessages({
|
||||
library(tidyverse)
|
||||
library(lubridate)
|
||||
library(here)
|
||||
})
|
||||
|
||||
# it will crash here if CI_all is empty and will not overwrite the rds rendering growth_model.R useless
|
||||
# if(nrow(CI_all) > 0){
|
||||
CI_all <- CI_all %>%
|
||||
group_by(model) %>%
|
||||
mutate(CI_per_day = FitData - lag(FitData), cumulative_CI = cumsum(FitData))
|
||||
# }
|
||||
|
||||
|
||||
|
||||
saveRDS(CI_all, here(cumulative_CI_vals_dir,"All_pivots_Cumulative_CI_quadrant_year_v2.rds"))
|
||||
message('rds saved')
|
||||
# 2. Main function to handle interpolation
|
||||
# -------------------------------------
|
||||
main <- function() {
|
||||
# Process command line arguments
|
||||
args <- commandArgs(trailingOnly = TRUE)
|
||||
|
||||
# Get project directory from arguments or use default
|
||||
if (length(args) >= 1 && !is.na(args[1])) {
|
||||
project_dir <- as.character(args[1])
|
||||
} else {
|
||||
project_dir <- "chemba"
|
||||
message("No project_dir provided. Using default:", project_dir)
|
||||
}
|
||||
|
||||
# Initialize project configuration and load utility functions
|
||||
source("r_app/parameters_project.R")
|
||||
source("r_app/growth_model_utils.R")
|
||||
|
||||
log_message("Starting CI growth model interpolation")
|
||||
|
||||
# Load and process the data
|
||||
tryCatch({
|
||||
# Load the combined CI data
|
||||
CI_data <- load_combined_ci_data(cumulative_CI_vals_dir)
|
||||
|
||||
# Validate harvesting data
|
||||
if (is.null(harvesting_data) || nrow(harvesting_data) == 0) {
|
||||
stop("No harvesting data available")
|
||||
}
|
||||
|
||||
# Get the years from harvesting data
|
||||
years <- harvesting_data %>%
|
||||
filter(!is.na(season_start)) %>%
|
||||
distinct(year) %>%
|
||||
pull(year)
|
||||
|
||||
log_message(paste("Processing data for years:", paste(years, collapse = ", ")))
|
||||
|
||||
# Generate interpolated CI data for each year and field
|
||||
CI_all <- generate_interpolated_ci_data(years, harvesting_data, CI_data)
|
||||
|
||||
# Calculate growth metrics and save the results
|
||||
if (nrow(CI_all) > 0) {
|
||||
# Add daily and cumulative metrics
|
||||
CI_all_with_metrics <- calculate_growth_metrics(CI_all)
|
||||
|
||||
# Save the processed data
|
||||
save_growth_model(
|
||||
CI_all_with_metrics,
|
||||
cumulative_CI_vals_dir,
|
||||
"All_pivots_Cumulative_CI_quadrant_year_v2.rds"
|
||||
)
|
||||
} else {
|
||||
log_message("No CI data was generated after interpolation", level = "WARNING")
|
||||
}
|
||||
|
||||
log_message("Growth model interpolation completed successfully")
|
||||
|
||||
}, error = function(e) {
|
||||
log_message(paste("Error in growth model interpolation:", e$message), level = "ERROR")
|
||||
stop(e$message)
|
||||
})
|
||||
}
|
||||
|
||||
# Run the main function if the script is executed directly
|
||||
main()
|
||||
|
|
|
|||
|
|
@ -1,137 +1,101 @@
|
|||
library(sf)
|
||||
library(terra)
|
||||
library(tidyverse)
|
||||
library(lubridate)
|
||||
# filepath: c:\Users\timon\Resilience BV\4020 SCane ESA DEMO - Documenten\General\4020 SCDEMO Team\4020 TechnicalData\WP3\smartcane\r_app\mosaic_creation.R
|
||||
#
|
||||
# MOSAIC_CREATION.R
|
||||
# ===============
|
||||
# This script creates weekly mosaics from daily satellite imagery.
|
||||
# It handles command-line arguments and initiates the mosaic creation process.
|
||||
#
|
||||
# Usage: Rscript mosaic_creation.R [end_date] [offset] [project_dir] [file_name]
|
||||
# - end_date: End date for processing (YYYY-MM-DD format)
|
||||
# - offset: Number of days to look back from end_date
|
||||
# - project_dir: Project directory name (e.g., "chemba")
|
||||
# - file_name: Optional custom output file name
|
||||
#
|
||||
|
||||
# Vang alle command line argumenten op
|
||||
args <- commandArgs(trailingOnly = TRUE)
|
||||
# 1. Load required packages
|
||||
# -----------------------
|
||||
suppressPackageStartupMessages({
|
||||
library(sf)
|
||||
library(terra)
|
||||
library(tidyverse)
|
||||
library(lubridate)
|
||||
library(here)
|
||||
})
|
||||
|
||||
# Controleer of er ten minste één argument is doorgegeven
|
||||
if (length(args) == 0) {
|
||||
stop("Geen argumenten doorgegeven aan het script")
|
||||
}
|
||||
|
||||
# Converteer het eerste argument naar een numerieke waarde
|
||||
end_date <- as.Date(args[1])
|
||||
|
||||
|
||||
offset <- as.numeric(args[2])
|
||||
# Controleer of weeks_ago een geldig getal is
|
||||
if (is.na(offset)) {
|
||||
# stop("Het argument is geen geldig getal")
|
||||
offset <- 7
|
||||
}
|
||||
|
||||
# Converteer het tweede argument naar een string waarde
|
||||
project_dir <- as.character(args[3])
|
||||
|
||||
# Controleer of data_dir een geldige waarde is
|
||||
if (!is.character(project_dir)) {
|
||||
project_dir <- "chemba"
|
||||
}
|
||||
|
||||
source("parameters_project.R")
|
||||
source("mosaic_creation_utils.R")
|
||||
|
||||
week <- week(end_date)
|
||||
|
||||
dates <- date_list(end_date, offset)
|
||||
|
||||
file_name_tif <- as.character(args[4])
|
||||
if (is.na(file_name_tif)) {
|
||||
file_name_tif <- paste0("week_", sprintf("%02d", dates$week), "_", dates$year, ".tif")
|
||||
}
|
||||
|
||||
print(dates)
|
||||
print(file_name_tif)
|
||||
#load pivot geojson
|
||||
# pivot_sf_q <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% vect()
|
||||
|
||||
vrt_files <- list.files(here(daily_vrt),full.names = T)
|
||||
vrt_list <- map(dates$days_filter, ~ vrt_files[grepl(pattern = .x, x = vrt_files)]) %>%
|
||||
compact() %>%
|
||||
flatten_chr()
|
||||
|
||||
raster_files_final <- list.files(merged_final,full.names = T, pattern = ".tif")
|
||||
|
||||
if (length(vrt_list) > 0) {
|
||||
print("vrt list made, preparing mosaic creation by counting cloud cover")
|
||||
# 2. Process command line arguments and run mosaic creation
|
||||
# ------------------------------------------------------
|
||||
main <- function() {
|
||||
# Capture command line arguments
|
||||
args <- commandArgs(trailingOnly = TRUE)
|
||||
|
||||
total_pix_area <-
|
||||
rast(vrt_list[1]) %>% terra::subset(1) %>% setValues(1) %>%
|
||||
crop(field_boundaries, mask = TRUE) %>%
|
||||
global(., fun = "notNA") #%>%
|
||||
|
||||
layer_5_list <- purrr::map(vrt_list, function(vrt_list) {
|
||||
rast(vrt_list[1]) %>% terra::subset(1)
|
||||
}) %>% rast()
|
||||
|
||||
missing_pixels_count <-
|
||||
layer_5_list %>% global(., fun = "notNA") %>%
|
||||
mutate(
|
||||
total_pixels = total_pix_area$notNA,
|
||||
missing_pixels_percentage = round(100 - ((
|
||||
notNA / total_pix_area$notNA
|
||||
) * 100)),
|
||||
thres_5perc = as.integer(missing_pixels_percentage < 5),
|
||||
thres_40perc = as.integer(missing_pixels_percentage < 45)
|
||||
)
|
||||
|
||||
index_5perc <- which(missing_pixels_count$thres_5perc == max(missing_pixels_count$thres_5perc))
|
||||
index_40perc <- which(missing_pixels_count$thres_40perc == max(missing_pixels_count$thres_40perc))
|
||||
|
||||
## Create mosaic
|
||||
|
||||
if (sum(missing_pixels_count$thres_5perc) > 1) {
|
||||
message("More than 1 raster without clouds (<5%), max composite made")
|
||||
|
||||
cloudy_rasters_list <- vrt_list[index_5perc]
|
||||
|
||||
rsrc <- sprc(cloudy_rasters_list)
|
||||
x <- terra::mosaic(rsrc, fun = "max")
|
||||
names(x) <- c("Red", "Green", "Blue", "NIR", "CI")
|
||||
|
||||
} else if (sum(missing_pixels_count$thres_5perc) == 1) {
|
||||
message("Only 1 raster without clouds (<5%)")
|
||||
|
||||
x <- rast(vrt_list[index_5perc[1]])
|
||||
names(x) <- c("Red", "Green", "Blue", "NIR", "CI")
|
||||
|
||||
} else if (sum(missing_pixels_count$thres_40perc) > 1) {
|
||||
message("More than 1 image contains clouds, composite made of <40% cloud cover images")
|
||||
|
||||
cloudy_rasters_list <- vrt_list[index_40perc]
|
||||
|
||||
rsrc <- sprc(cloudy_rasters_list)
|
||||
x <- mosaic(rsrc, fun = "max")
|
||||
names(x) <- c("Red", "Green", "Blue", "NIR", "CI")
|
||||
|
||||
} else if (sum(missing_pixels_count$thres_40perc) == 1) {
|
||||
message("Only 1 image available but contains clouds")
|
||||
|
||||
x <- rast(vrt_list[index_40perc[1]])
|
||||
names(x) <- c("Red", "Green", "Blue", "NIR", "CI")
|
||||
|
||||
} else if (sum(missing_pixels_count$thres_40perc) == 0) {
|
||||
message("No cloud free images available, all images combined")
|
||||
|
||||
rsrc <- sprc(vrt_list)
|
||||
x <- mosaic(rsrc, fun = "max")
|
||||
names(x) <- c("Red", "Green", "Blue", "NIR", "CI")
|
||||
|
||||
# Process end_date argument with default
|
||||
if (length(args) >= 1 && !is.na(args[1])) {
|
||||
end_date <- as.Date(args[1])
|
||||
if (is.na(end_date)) {
|
||||
message("Invalid end_date provided. Using current date.")
|
||||
#end_date <- Sys.Date()
|
||||
end_date <- "2023-10-01" # Default date for testing
|
||||
}
|
||||
} else {
|
||||
# Default to current date if no argument is provided
|
||||
#end_date <- Sys.Date()
|
||||
end_date <- "2023-10-01" # Default date for testing
|
||||
message("No end_date provided. Using current date:", format(end_date))
|
||||
}
|
||||
|
||||
} else{
|
||||
message("No images available this week, empty mosaic created")
|
||||
# Process offset argument with default
|
||||
if (length(args) >= 2 && !is.na(args[2])) {
|
||||
offset <- as.numeric(args[2])
|
||||
if (is.na(offset) || offset <= 0) {
|
||||
message("Invalid offset provided. Using default (7 days).")
|
||||
offset <- 7
|
||||
}
|
||||
} else {
|
||||
# Default to 7 days if no argument is provided
|
||||
offset <- 7
|
||||
message("No offset provided. Using default:", offset, "days")
|
||||
}
|
||||
|
||||
x <- rast(raster_files_final[1]) %>% setValues(0) %>%
|
||||
crop(field_boundaries, mask = TRUE)
|
||||
# Process project_dir argument with default
|
||||
if (length(args) >= 3 && !is.na(args[3])) {
|
||||
project_dir <- as.character(args[3])
|
||||
} else {
|
||||
# Default project directory
|
||||
project_dir <- "chemba"
|
||||
message("No project_dir provided. Using default:", project_dir)
|
||||
}
|
||||
|
||||
names(x) <- c("Red", "Green", "Blue", "NIR", "CI")
|
||||
# 3. Initialize project configuration
|
||||
# --------------------------------
|
||||
source("r_app/parameters_project.R")
|
||||
source("r_app/mosaic_creation_utils.R")
|
||||
|
||||
# 4. Generate date range for processing
|
||||
# ---------------------------------
|
||||
dates <- date_list(end_date, offset)
|
||||
log_message(paste("Processing data for week", dates$week, "of", dates$year))
|
||||
|
||||
# Create output filename
|
||||
file_name_tif <- if (length(args) >= 4 && !is.na(args[4])) {
|
||||
as.character(args[4])
|
||||
} else {
|
||||
paste0("week_", sprintf("%02d", dates$week), "_", dates$year, ".tif")
|
||||
}
|
||||
|
||||
log_message(paste("Output will be saved as:", file_name_tif))
|
||||
|
||||
# 5. Create weekly mosaic using the function from utils
|
||||
# -------------------------------------------------
|
||||
create_weekly_mosaic(
|
||||
dates = dates,
|
||||
field_boundaries = field_boundaries,
|
||||
daily_vrt_dir = daily_vrt,
|
||||
merged_final_dir = merged_final,
|
||||
output_dir = weekly_CI_mosaic,
|
||||
file_name_tif = file_name_tif,
|
||||
create_plots = TRUE
|
||||
)
|
||||
}
|
||||
plot(x$CI, main = paste("CI map ", dates$week))
|
||||
plotRGB(x, main = paste("RGB map ", dates$week))
|
||||
|
||||
file_path_tif <- here(weekly_CI_mosaic ,file_name_tif)
|
||||
writeRaster(x, file_path_tif, overwrite=TRUE)
|
||||
message("Raster written/made at: ", file_path_tif)
|
||||
# Run the main function if the script is executed directly
|
||||
main()
|
||||
|
|
|
|||
|
|
@ -1,13 +1,309 @@
|
|||
# utils for mosaic creation
|
||||
# MOSAIC_CREATION_UTILS.R
|
||||
# ======================
|
||||
# Utility functions for creating weekly mosaics from daily satellite imagery.
|
||||
# These functions support cloud cover assessment, date handling, and mosaic creation.
|
||||
|
||||
date_list <- function(end_date, offset){
|
||||
offset <- as.numeric(offset) - 1
|
||||
end_date <- as.Date(end_date)
|
||||
#' Safe logging function that works whether log_message exists or not
|
||||
#'
|
||||
#' @param message The message to log
|
||||
#' @param level The log level (default: "INFO")
|
||||
#' @return NULL (used for side effects)
|
||||
#'
|
||||
safe_log <- function(message, level = "INFO") {
|
||||
if (exists("log_message")) {
|
||||
log_message(message, level)
|
||||
} else {
|
||||
if (level %in% c("ERROR", "WARNING")) {
|
||||
warning(message)
|
||||
} else {
|
||||
message(message)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#' Generate a sequence of dates for processing
|
||||
#'
|
||||
#' @param end_date The end date for the sequence (Date object)
|
||||
#' @param offset Number of days to look back from end_date
|
||||
#' @return A list containing week number, year, and a sequence of dates for filtering
|
||||
#'
|
||||
date_list <- function(end_date, offset) {
|
||||
# Input validation
|
||||
if (!lubridate::is.Date(end_date)) {
|
||||
end_date <- as.Date(end_date)
|
||||
if (is.na(end_date)) {
|
||||
stop("Invalid end_date provided. Expected a Date object or a string convertible to Date.")
|
||||
}
|
||||
}
|
||||
|
||||
offset <- as.numeric(offset)
|
||||
if (is.na(offset) || offset < 1) {
|
||||
stop("Invalid offset provided. Expected a positive number.")
|
||||
}
|
||||
|
||||
# Calculate date range
|
||||
offset <- offset - 1 # Adjust offset to include 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")
|
||||
# Extract week and year information
|
||||
week <- lubridate::week(start_date)
|
||||
year <- lubridate::year(start_date)
|
||||
|
||||
return(list("week" = week, "year" = year, "days_filter" = days_filter))
|
||||
# Generate sequence of dates
|
||||
days_filter <- seq(from = start_date, to = end_date, by = "day")
|
||||
days_filter <- format(days_filter, "%Y-%m-%d") # Format for consistent filtering
|
||||
|
||||
# Log the date range
|
||||
safe_log(paste("Date range generated from", start_date, "to", end_date))
|
||||
|
||||
return(list(
|
||||
"week" = week,
|
||||
"year" = year,
|
||||
"days_filter" = days_filter,
|
||||
"start_date" = start_date,
|
||||
"end_date" = end_date
|
||||
))
|
||||
}
|
||||
|
||||
#' Create a weekly mosaic from available VRT files
|
||||
#'
|
||||
#' @param dates List from date_list() with date range info
|
||||
#' @param field_boundaries Field boundaries for image cropping
|
||||
#' @param daily_vrt_dir Directory containing VRT files
|
||||
#' @param merged_final_dir Directory with merged final rasters
|
||||
#' @param output_dir Output directory for weekly mosaics
|
||||
#' @param file_name_tif Output filename for the mosaic
|
||||
#' @param create_plots Whether to create visualization plots (default: TRUE)
|
||||
#' @return The file path of the saved mosaic
|
||||
#'
|
||||
create_weekly_mosaic <- function(dates, field_boundaries, daily_vrt_dir,
|
||||
merged_final_dir, output_dir, file_name_tif,
|
||||
create_plots = TRUE) {
|
||||
# Find VRT files for the specified date range
|
||||
vrt_list <- find_vrt_files(daily_vrt_dir, dates)
|
||||
|
||||
# Find final raster files for fallback
|
||||
raster_files_final <- list.files(merged_final_dir, full.names = TRUE, pattern = "\\.tif$")
|
||||
|
||||
# Process the mosaic if VRT files are available
|
||||
if (length(vrt_list) > 0) {
|
||||
safe_log("VRT list created, assessing cloud cover for mosaic creation")
|
||||
|
||||
# Calculate cloud cover statistics
|
||||
missing_pixels_count <- count_cloud_coverage(vrt_list, field_boundaries)
|
||||
|
||||
# Create mosaic based on cloud cover assessment
|
||||
mosaic <- create_mosaic(vrt_list, missing_pixels_count, field_boundaries, raster_files_final)
|
||||
|
||||
} else {
|
||||
safe_log("No VRT files available for the date range, creating empty mosaic", "WARNING")
|
||||
|
||||
# Create empty mosaic if no files are available
|
||||
if (length(raster_files_final) == 0) {
|
||||
stop("No VRT files or final raster files available to create mosaic")
|
||||
}
|
||||
|
||||
mosaic <- terra::rast(raster_files_final[1]) %>%
|
||||
terra::setValues(0) %>%
|
||||
terra::crop(field_boundaries, mask = TRUE)
|
||||
|
||||
names(mosaic) <- c("Red", "Green", "Blue", "NIR", "CI")
|
||||
}
|
||||
|
||||
# Save the mosaic
|
||||
file_path <- save_mosaic(mosaic, output_dir, file_name_tif, create_plots)
|
||||
|
||||
safe_log(paste("Weekly mosaic processing completed for week", dates$week))
|
||||
|
||||
return(file_path)
|
||||
}
|
||||
|
||||
#' Find VRT files within a date range
|
||||
#'
|
||||
#' @param vrt_directory Directory containing VRT files
|
||||
#' @param dates List from date_list() function containing days_filter
|
||||
#' @return Character vector of VRT file paths
|
||||
#'
|
||||
find_vrt_files <- function(vrt_directory, dates) {
|
||||
# Get all VRT files in directory
|
||||
vrt_files <- list.files(here::here(vrt_directory), full.names = TRUE)
|
||||
|
||||
if (length(vrt_files) == 0) {
|
||||
warning("No VRT files found in directory: ", vrt_directory)
|
||||
return(character(0))
|
||||
}
|
||||
|
||||
# Filter files by dates
|
||||
vrt_list <- purrr::map(dates$days_filter, ~ vrt_files[grepl(pattern = .x, x = vrt_files)]) %>%
|
||||
purrr::compact() %>%
|
||||
purrr::flatten_chr()
|
||||
|
||||
# Log results
|
||||
safe_log(paste("Found", length(vrt_list), "VRT files for the date range"))
|
||||
|
||||
return(vrt_list)
|
||||
}
|
||||
|
||||
#' Count missing pixels (clouds) in rasters
|
||||
#'
|
||||
#' @param vrt_list List of VRT files to analyze
|
||||
#' @param field_boundaries Field boundaries vector for masking
|
||||
#' @return Data frame with cloud coverage statistics
|
||||
#'
|
||||
count_cloud_coverage <- function(vrt_list, field_boundaries) {
|
||||
if (length(vrt_list) == 0) {
|
||||
warning("No VRT files provided for cloud coverage calculation")
|
||||
return(NULL)
|
||||
}
|
||||
|
||||
tryCatch({
|
||||
# Calculate total pixel area using the first VRT file
|
||||
total_pix_area <- terra::rast(vrt_list[1]) %>%
|
||||
terra::subset(1) %>%
|
||||
terra::setValues(1) %>%
|
||||
terra::crop(field_boundaries, mask = TRUE) %>%
|
||||
terra::global(fun = "notNA")
|
||||
|
||||
# Extract layer 1 from all rasters (for cloud detection)
|
||||
layer_5_list <- purrr::map(vrt_list, function(file) {
|
||||
terra::rast(file) %>% terra::subset(1)
|
||||
}) %>% terra::rast()
|
||||
|
||||
# Calculate percentage of missing pixels (clouds)
|
||||
missing_pixels_count <- terra::global(layer_5_list, fun = "notNA") %>%
|
||||
dplyr::mutate(
|
||||
total_pixels = total_pix_area$notNA,
|
||||
missing_pixels_percentage = round(100 - ((notNA / total_pix_area$notNA) * 100)),
|
||||
thres_5perc = as.integer(missing_pixels_percentage < 5),
|
||||
thres_40perc = as.integer(missing_pixels_percentage < 45)
|
||||
)
|
||||
|
||||
# Log results
|
||||
safe_log(paste(
|
||||
"Cloud cover assessment completed for", length(vrt_list), "files.",
|
||||
sum(missing_pixels_count$thres_5perc), "files with <5% cloud cover,",
|
||||
sum(missing_pixels_count$thres_40perc), "files with <45% cloud cover"
|
||||
))
|
||||
|
||||
return(missing_pixels_count)
|
||||
}, error = function(e) {
|
||||
warning("Error in cloud coverage calculation: ", e$message)
|
||||
return(NULL)
|
||||
})
|
||||
}
|
||||
|
||||
#' Create a mosaic from VRT files based on cloud coverage
|
||||
#'
|
||||
#' @param vrt_list List of VRT files to create mosaic from
|
||||
#' @param missing_pixels_count Cloud coverage statistics from count_cloud_coverage()
|
||||
#' @param field_boundaries Field boundaries vector for masking (optional)
|
||||
#' @param raster_files_final List of processed raster files to use as fallback
|
||||
#' @return A SpatRaster object with the mosaic
|
||||
#'
|
||||
create_mosaic <- function(vrt_list, missing_pixels_count, field_boundaries = NULL, raster_files_final = NULL) {
|
||||
# If no VRT files, create an empty mosaic
|
||||
if (length(vrt_list) == 0) {
|
||||
if (length(raster_files_final) == 0 || is.null(field_boundaries)) {
|
||||
stop("No VRT files available and no fallback raster files or field boundaries provided")
|
||||
}
|
||||
|
||||
safe_log("No images available for this period, creating empty mosaic", "WARNING")
|
||||
|
||||
x <- terra::rast(raster_files_final[1]) %>%
|
||||
terra::setValues(0) %>%
|
||||
terra::crop(field_boundaries, mask = TRUE)
|
||||
|
||||
names(x) <- c("Red", "Green", "Blue", "NIR", "CI")
|
||||
return(x)
|
||||
}
|
||||
|
||||
# If missing pixel count was not calculated, use all files
|
||||
if (is.null(missing_pixels_count)) {
|
||||
safe_log("No cloud coverage data available, using all images", "WARNING")
|
||||
|
||||
rsrc <- terra::sprc(vrt_list)
|
||||
x <- terra::mosaic(rsrc, fun = "max")
|
||||
names(x) <- c("Red", "Green", "Blue", "NIR", "CI")
|
||||
return(x)
|
||||
}
|
||||
|
||||
# Determine best rasters to use based on cloud coverage
|
||||
index_5perc <- which(missing_pixels_count$thres_5perc == max(missing_pixels_count$thres_5perc))
|
||||
index_40perc <- which(missing_pixels_count$thres_40perc == max(missing_pixels_count$thres_40perc))
|
||||
|
||||
# Create mosaic based on available cloud-free images
|
||||
if (sum(missing_pixels_count$thres_5perc) > 1) {
|
||||
safe_log("Creating max composite from multiple cloud-free images (<5% clouds)")
|
||||
|
||||
cloudy_rasters_list <- vrt_list[index_5perc]
|
||||
rsrc <- terra::sprc(cloudy_rasters_list)
|
||||
x <- terra::mosaic(rsrc, fun = "max")
|
||||
|
||||
} else if (sum(missing_pixels_count$thres_5perc) == 1) {
|
||||
safe_log("Using single cloud-free image (<5% clouds)")
|
||||
|
||||
x <- terra::rast(vrt_list[index_5perc[1]])
|
||||
|
||||
} else if (sum(missing_pixels_count$thres_40perc) > 1) {
|
||||
safe_log("Creating max composite from partially cloudy images (<40% clouds)", "WARNING")
|
||||
|
||||
cloudy_rasters_list <- vrt_list[index_40perc]
|
||||
rsrc <- terra::sprc(cloudy_rasters_list)
|
||||
x <- terra::mosaic(rsrc, fun = "max")
|
||||
|
||||
} else if (sum(missing_pixels_count$thres_40perc) == 1) {
|
||||
safe_log("Using single partially cloudy image (<40% clouds)", "WARNING")
|
||||
|
||||
x <- terra::rast(vrt_list[index_40perc[1]])
|
||||
|
||||
} else {
|
||||
safe_log("No cloud-free images available, using all images", "WARNING")
|
||||
|
||||
rsrc <- terra::sprc(vrt_list)
|
||||
x <- terra::mosaic(rsrc, fun = "max")
|
||||
}
|
||||
|
||||
# Set consistent layer names
|
||||
names(x) <- c("Red", "Green", "Blue", "NIR", "CI")
|
||||
return(x)
|
||||
}
|
||||
|
||||
#' Save a mosaic raster to disk
|
||||
#'
|
||||
#' @param mosaic_raster A SpatRaster object to save
|
||||
#' @param output_dir Directory to save the output
|
||||
#' @param file_name Filename for the output raster
|
||||
#' @param plot_result Whether to create visualizations (default: FALSE)
|
||||
#' @return The file path of the saved raster
|
||||
#'
|
||||
save_mosaic <- function(mosaic_raster, output_dir, file_name, plot_result = FALSE) {
|
||||
# Validate input
|
||||
if (is.null(mosaic_raster)) {
|
||||
stop("No mosaic raster provided to save")
|
||||
}
|
||||
|
||||
# Create output directory if it doesn't exist
|
||||
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
|
||||
|
||||
# Create full file path
|
||||
file_path <- here::here(output_dir, file_name)
|
||||
|
||||
# Save raster
|
||||
terra::writeRaster(mosaic_raster, file_path, overwrite = TRUE)
|
||||
|
||||
# Create plots if requested
|
||||
if (plot_result) {
|
||||
if ("CI" %in% names(mosaic_raster)) {
|
||||
terra::plot(mosaic_raster$CI, main = paste("CI map", file_name))
|
||||
}
|
||||
|
||||
if (all(c("Red", "Green", "Blue") %in% names(mosaic_raster))) {
|
||||
terra::plotRGB(mosaic_raster, main = paste("RGB map", file_name))
|
||||
}
|
||||
}
|
||||
|
||||
# Log save completion
|
||||
safe_log(paste("Mosaic saved to:", file_path))
|
||||
|
||||
return(file_path)
|
||||
}
|
||||
|
|
|
|||
110
r_app/packages.R
Normal file
110
r_app/packages.R
Normal file
|
|
@ -0,0 +1,110 @@
|
|||
# packages.R
|
||||
#
|
||||
# PACKAGE MANAGEMENT FOR SMARTCANE
|
||||
# ===============================
|
||||
# This script centralizes all package dependencies for the SmartCane project.
|
||||
# It installs missing packages and loads all required libraries.
|
||||
#
|
||||
|
||||
#' Check and install packages if needed
|
||||
#'
|
||||
#' @param pkg_list List of packages to check and install
|
||||
#' @param install_missing Whether to install missing packages
|
||||
#' @return Vector of packages that couldn't be installed (if any)
|
||||
#'
|
||||
check_and_install_packages <- function(pkg_list, install_missing = TRUE) {
|
||||
# Check which packages are already installed
|
||||
is_installed <- pkg_list %in% rownames(installed.packages())
|
||||
missing_pkgs <- pkg_list[!is_installed]
|
||||
|
||||
# Install missing packages if requested
|
||||
failed_pkgs <- character(0)
|
||||
if (length(missing_pkgs) > 0) {
|
||||
if (install_missing) {
|
||||
message("Installing ", length(missing_pkgs), " missing packages...")
|
||||
for (pkg in missing_pkgs) {
|
||||
tryCatch({
|
||||
install.packages(pkg, repos = "https://cran.rstudio.com/", dependencies = TRUE)
|
||||
message(" Installed: ", pkg)
|
||||
}, error = function(e) {
|
||||
warning("Failed to install package: ", pkg)
|
||||
warning("Error: ", e$message)
|
||||
failed_pkgs <<- c(failed_pkgs, pkg)
|
||||
})
|
||||
}
|
||||
} else {
|
||||
message("The following packages are required but not installed:")
|
||||
message(paste(missing_pkgs, collapse = ", "))
|
||||
failed_pkgs <- missing_pkgs
|
||||
}
|
||||
} else {
|
||||
message("All required packages are already installed.")
|
||||
}
|
||||
|
||||
return(failed_pkgs)
|
||||
}
|
||||
|
||||
#' Load all required packages for SmartCane project
|
||||
#'
|
||||
#' @param verbose Whether to show messages during loading
|
||||
#' @return Logical indicating success (TRUE if all packages loaded)
|
||||
#'
|
||||
load_smartcane_packages <- function(verbose = FALSE) {
|
||||
# Define all required packages
|
||||
required_packages <- c(
|
||||
# Geospatial packages
|
||||
"sf", # Simple Features for spatial vector data
|
||||
"terra", # Raster data processing
|
||||
"exactextractr", # Fast extraction from rasters
|
||||
"raster", # Legacy raster package (for compatibility)
|
||||
|
||||
# Data manipulation
|
||||
"tidyverse", # Collection of data manipulation packages
|
||||
"lubridate", # Date manipulation
|
||||
"readxl", # Excel file reading
|
||||
"stringr", # String manipulation
|
||||
"purrr", # Functional programming tools
|
||||
|
||||
# Visualization
|
||||
"ggplot2", # Advanced plotting
|
||||
"leaflet", # Interactive maps
|
||||
"plotly", # Interactive plots
|
||||
|
||||
# Project management
|
||||
"here", # Path handling
|
||||
|
||||
# Document generation
|
||||
"knitr", # Dynamic report generation
|
||||
"rmarkdown" # R Markdown processing
|
||||
)
|
||||
|
||||
# Check and install missing packages
|
||||
failed_pkgs <- check_and_install_packages(required_packages)
|
||||
|
||||
# Load all installed packages
|
||||
success <- TRUE
|
||||
for (pkg in setdiff(required_packages, failed_pkgs)) {
|
||||
tryCatch({
|
||||
if (verbose) message("Loading package: ", pkg)
|
||||
suppressPackageStartupMessages(library(pkg, character.only = TRUE))
|
||||
}, error = function(e) {
|
||||
warning("Failed to load package: ", pkg)
|
||||
warning("Error: ", e$message)
|
||||
success <- FALSE
|
||||
})
|
||||
}
|
||||
|
||||
# Report any issues
|
||||
if (length(failed_pkgs) > 0) {
|
||||
warning("The following packages could not be installed: ",
|
||||
paste(failed_pkgs, collapse = ", "))
|
||||
success <- FALSE
|
||||
}
|
||||
|
||||
return(success)
|
||||
}
|
||||
|
||||
# Run the loading function if the script is sourced directly
|
||||
if (!exists("skip_package_loading") || !skip_package_loading) {
|
||||
load_smartcane_packages()
|
||||
}
|
||||
|
|
@ -1,82 +1,225 @@
|
|||
library(here)
|
||||
library('readxl')
|
||||
#chemba
|
||||
# filepath: c:\Users\timon\Resilience BV\4020 SCane ESA DEMO - Documenten\General\4020 SCDEMO Team\4020 TechnicalData\WP3\smartcane\r_app\parameters_project.R
|
||||
#
|
||||
# PARAMETERS_PROJECT.R
|
||||
# ====================
|
||||
# This script defines project parameters, directory structures, and loads field boundaries.
|
||||
# It establishes all the necessary paths and creates required directories for the SmartCane project.
|
||||
|
||||
laravel_storage_dir <- here("laravel_app/storage/app", project_dir)
|
||||
reports_dir <- here(laravel_storage_dir, "reports")
|
||||
log_dir <- here(laravel_storage_dir, "logs")
|
||||
planet_tif_folder <- here(laravel_storage_dir, "merged_tif")
|
||||
merged_final <- here(laravel_storage_dir, "merged_final_tif")
|
||||
planet_tif_folder <- here(laravel_storage_dir, "merged_tif")
|
||||
merged_final <- here(laravel_storage_dir, "merged_final_tif")
|
||||
data_dir <- here(laravel_storage_dir, "Data")
|
||||
extracted_CI_dir <- here(data_dir, "extracted_ci")
|
||||
daily_CI_vals_dir <- here(extracted_CI_dir, "daily_vals")
|
||||
cumulative_CI_vals_dir <- here(extracted_CI_dir, "cumulative_vals")
|
||||
weekly_CI_mosaic <- here(laravel_storage_dir, "weekly_mosaic")
|
||||
daily_vrt <- here(data_dir, "vrt")
|
||||
harvest_dir <- here(data_dir, "HarvestData")
|
||||
# 1. Load required libraries
|
||||
# -------------------------
|
||||
suppressPackageStartupMessages({
|
||||
library(here)
|
||||
library(readxl)
|
||||
library(sf)
|
||||
library(dplyr)
|
||||
library(tidyr)
|
||||
})
|
||||
|
||||
dir.create(here(laravel_storage_dir), showWarnings = FALSE)
|
||||
dir.create(here(reports_dir), showWarnings = FALSE)
|
||||
dir.create(here(data_dir), showWarnings = FALSE)
|
||||
dir.create(here(log_dir), showWarnings = FALSE)
|
||||
dir.create(here(extracted_CI_dir), showWarnings = FALSE)
|
||||
dir.create(here(daily_CI_vals_dir), showWarnings = FALSE)
|
||||
dir.create(here(cumulative_CI_vals_dir), showWarnings = FALSE)
|
||||
dir.create(here(weekly_CI_mosaic),showWarnings = FALSE)
|
||||
dir.create(here(daily_vrt), showWarnings = FALSE)
|
||||
dir.create(merged_final,showWarnings = FALSE)
|
||||
dir.create(harvest_dir,showWarnings = FALSE)
|
||||
|
||||
field_boundaries_sf <- st_read(here(data_dir, "pivot.geojson"), crs = 4326)
|
||||
|
||||
names(field_boundaries_sf) <- c("field", "sub_field", "geometry")
|
||||
|
||||
field_boundaries <- field_boundaries_sf %>% terra::vect()
|
||||
|
||||
harvesting_data <- read_excel(here(data_dir, "harvest.xlsx")) %>%
|
||||
dplyr::select(
|
||||
c(
|
||||
"field",
|
||||
"sub_field",
|
||||
"year",
|
||||
"season_start",
|
||||
"season_end",
|
||||
"age",
|
||||
"sub_area",
|
||||
"tonnage_ha"
|
||||
)
|
||||
) %>%
|
||||
mutate(
|
||||
field = as.character(field),
|
||||
sub_field = as.character(sub_field),
|
||||
year = as.numeric(year),
|
||||
season_start = as.Date(season_start, format="%d/%m/%Y"),
|
||||
season_end = as.Date(season_end, format="%d/%m/%Y"),
|
||||
age = as.numeric(age),
|
||||
sub_area = as.character(sub_area),
|
||||
tonnage_ha = as.numeric(tonnage_ha)
|
||||
) %>%
|
||||
mutate(
|
||||
season_end = case_when(season_end > Sys.Date() ~ Sys.Date(),
|
||||
is.na(season_end) ~ Sys.Date(),
|
||||
TRUE ~ season_end),
|
||||
age = round(as.numeric(season_end - season_start) / 7, 0)
|
||||
# 2. Define project directory structure
|
||||
# -----------------------------------
|
||||
setup_project_directories <- function(project_dir) {
|
||||
# Base directories
|
||||
laravel_storage_dir <- here("laravel_app/storage/app", project_dir)
|
||||
|
||||
# Main subdirectories
|
||||
dirs <- list(
|
||||
reports = here(laravel_storage_dir, "reports"),
|
||||
logs = here(laravel_storage_dir, "logs"),
|
||||
data = here(laravel_storage_dir, "Data"),
|
||||
tif = list(
|
||||
merged = here(laravel_storage_dir, "merged_tif"),
|
||||
final = here(laravel_storage_dir, "merged_final_tif")
|
||||
),
|
||||
weekly_mosaic = here(laravel_storage_dir, "weekly_mosaic"),
|
||||
extracted_ci = list(
|
||||
base = here(laravel_storage_dir, "Data/extracted_ci"),
|
||||
daily = here(laravel_storage_dir, "Data/extracted_ci/daily_vals"),
|
||||
cumulative = here(laravel_storage_dir, "Data/extracted_ci/cumulative_vals")
|
||||
),
|
||||
vrt = here(laravel_storage_dir, "Data/vrt"),
|
||||
harvest = here(laravel_storage_dir, "Data/HarvestData")
|
||||
)
|
||||
|
||||
log_file <- here(log_dir, paste0(format(Sys.Date(), "%Y%m%d"), ".log"))
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# Create a logging function
|
||||
log_message <- function(message) {
|
||||
cat(message, "\n", file = log_file, append = TRUE)
|
||||
|
||||
# Create all directories
|
||||
for (dir_path in unlist(dirs)) {
|
||||
dir.create(dir_path, showWarnings = FALSE, recursive = TRUE)
|
||||
}
|
||||
|
||||
# Return directory structure for use in other functions
|
||||
return(list(
|
||||
laravel_storage_dir = laravel_storage_dir,
|
||||
reports_dir = dirs$reports,
|
||||
log_dir = dirs$logs,
|
||||
data_dir = dirs$data,
|
||||
planet_tif_folder = dirs$tif$merged,
|
||||
merged_final = dirs$tif$final,
|
||||
daily_CI_vals_dir = dirs$extracted_ci$daily,
|
||||
cumulative_CI_vals_dir = dirs$extracted_ci$cumulative,
|
||||
weekly_CI_mosaic = dirs$weekly_mosaic,
|
||||
daily_vrt = dirs$vrt,
|
||||
harvest_dir = dirs$harvest,
|
||||
extracted_CI_dir = dirs$extracted_ci$base
|
||||
))
|
||||
}
|
||||
|
||||
|
||||
log_head <- function(list) {
|
||||
log_message(paste(capture.output(str(head(list))), collapse = "\n"))
|
||||
# 3. Load field boundaries
|
||||
# ----------------------
|
||||
load_field_boundaries <- function(data_dir) {
|
||||
field_boundaries_path <- here(data_dir, "pivot.geojson")
|
||||
|
||||
if (!file.exists(field_boundaries_path)) {
|
||||
stop(paste("Field boundaries file not found at path:", field_boundaries_path))
|
||||
}
|
||||
|
||||
tryCatch({
|
||||
field_boundaries_sf <- st_read(field_boundaries_path, crs = 4326, quiet = TRUE)
|
||||
names(field_boundaries_sf) <- c("field", "sub_field", "geometry")
|
||||
field_boundaries <- terra::vect(field_boundaries_sf)
|
||||
|
||||
return(list(
|
||||
field_boundaries_sf = field_boundaries_sf,
|
||||
field_boundaries = field_boundaries
|
||||
))
|
||||
}, error = function(e) {
|
||||
stop(paste("Error loading field boundaries:", e$message))
|
||||
})
|
||||
}
|
||||
|
||||
# 4. Load harvesting data
|
||||
# ---------------------
|
||||
load_harvesting_data <- function(data_dir) {
|
||||
harvest_file <- here(data_dir, "harvest.xlsx")
|
||||
|
||||
if (!file.exists(harvest_file)) {
|
||||
warning(paste("Harvest data file not found at path:", harvest_file))
|
||||
return(NULL)
|
||||
}
|
||||
|
||||
tryCatch({
|
||||
harvesting_data <- read_excel(harvest_file) %>%
|
||||
dplyr::select(
|
||||
c(
|
||||
"field",
|
||||
"sub_field",
|
||||
"year",
|
||||
"season_start",
|
||||
"season_end",
|
||||
"age",
|
||||
"sub_area",
|
||||
"tonnage_ha"
|
||||
)
|
||||
) %>%
|
||||
mutate(
|
||||
field = as.character(field),
|
||||
sub_field = as.character(sub_field),
|
||||
year = as.numeric(year),
|
||||
season_start = as.Date(season_start, format="%d/%m/%Y"),
|
||||
season_end = as.Date(season_end, format="%d/%m/%Y"),
|
||||
age = as.numeric(age),
|
||||
sub_area = as.character(sub_area),
|
||||
tonnage_ha = as.numeric(tonnage_ha)
|
||||
) %>%
|
||||
mutate(
|
||||
season_end = case_when(
|
||||
season_end > Sys.Date() ~ Sys.Date(),
|
||||
is.na(season_end) ~ Sys.Date(),
|
||||
TRUE ~ season_end
|
||||
),
|
||||
age = round(as.numeric(season_end - season_start) / 7, 0)
|
||||
)
|
||||
|
||||
return(harvesting_data)
|
||||
}, error = function(e) {
|
||||
warning(paste("Error loading harvesting data:", e$message))
|
||||
return(NULL)
|
||||
})
|
||||
}
|
||||
|
||||
# 5. Define logging functions globally first
|
||||
# ---------------------------------------
|
||||
# Create a simple default log function in case setup_logging hasn't been called yet
|
||||
log_message <- function(message, level = "INFO") {
|
||||
timestamp <- format(Sys.time(), "%Y-%m-%d %H:%M:%S")
|
||||
formatted_message <- paste0("[", level, "] ", timestamp, " - ", message)
|
||||
cat(formatted_message, "\n")
|
||||
}
|
||||
|
||||
log_head <- function(list, level = "INFO") {
|
||||
log_message(paste(capture.output(str(head(list))), collapse = "\n"), level)
|
||||
}
|
||||
|
||||
# 6. Set up full logging system with file output
|
||||
# -------------------------------------------
|
||||
setup_logging <- function(log_dir) {
|
||||
log_file <- here(log_dir, paste0(format(Sys.Date(), "%Y%m%d"), ".log"))
|
||||
|
||||
# Create enhanced log functions
|
||||
log_message <- function(message, level = "INFO") {
|
||||
timestamp <- format(Sys.time(), "%Y-%m-%d %H:%M:%S")
|
||||
formatted_message <- paste0("[", level, "] ", timestamp, " - ", message)
|
||||
cat(formatted_message, "\n", file = log_file, append = TRUE)
|
||||
|
||||
# Also print to console for debugging
|
||||
if (level %in% c("ERROR", "WARNING")) {
|
||||
cat(formatted_message, "\n")
|
||||
}
|
||||
}
|
||||
|
||||
log_head <- function(list, level = "INFO") {
|
||||
log_message(paste(capture.output(str(head(list))), collapse = "\n"), level)
|
||||
}
|
||||
|
||||
# Update the global functions with the enhanced versions
|
||||
assign("log_message", log_message, envir = .GlobalEnv)
|
||||
assign("log_head", log_head, envir = .GlobalEnv)
|
||||
|
||||
return(list(
|
||||
log_file = log_file,
|
||||
log_message = log_message,
|
||||
log_head = log_head
|
||||
))
|
||||
}
|
||||
|
||||
# 7. Initialize the project
|
||||
# ----------------------
|
||||
# Export project directories and settings
|
||||
initialize_project <- function(project_dir) {
|
||||
# Set up directory structure
|
||||
dirs <- setup_project_directories(project_dir)
|
||||
|
||||
# Set up logging
|
||||
logging <- setup_logging(dirs$log_dir)
|
||||
|
||||
# Load field boundaries
|
||||
boundaries <- load_field_boundaries(dirs$data_dir)
|
||||
|
||||
# Load harvesting data
|
||||
harvesting_data <- load_harvesting_data(dirs$data_dir)
|
||||
|
||||
# Return all initialized components
|
||||
return(c(
|
||||
dirs,
|
||||
list(
|
||||
logging = logging,
|
||||
field_boundaries = boundaries$field_boundaries,
|
||||
field_boundaries_sf = boundaries$field_boundaries_sf,
|
||||
harvesting_data = harvesting_data
|
||||
)
|
||||
))
|
||||
}
|
||||
|
||||
# When script is sourced, initialize with the global project_dir variable if it exists
|
||||
if (exists("project_dir")) {
|
||||
# Now we can safely log before initialization
|
||||
log_message(paste("Initializing project with directory:", project_dir))
|
||||
|
||||
project_config <- initialize_project(project_dir)
|
||||
|
||||
# Expose all variables to the global environment
|
||||
list2env(project_config, envir = .GlobalEnv)
|
||||
|
||||
# Log project initialization completion
|
||||
log_message(paste("Project initialized with directory:", project_dir))
|
||||
} else {
|
||||
warning("project_dir variable not found. Please set project_dir before sourcing parameters_project.R")
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,3 +1,34 @@
|
|||
# REPORT_UTILS.R
|
||||
# =============
|
||||
# Utility functions for generating SmartCane reports with visualizations.
|
||||
# These functions support the creation of maps, charts and report elements
|
||||
# for the CI_report_dashboard_planet.Rmd document.
|
||||
|
||||
#' Safe logging function that works whether log_message exists or not
|
||||
#'
|
||||
#' @param message The message to log
|
||||
#' @param level The log level (default: "INFO")
|
||||
#' @return NULL (used for side effects)
|
||||
#'
|
||||
safe_log <- function(message, level = "INFO") {
|
||||
if (exists("log_message")) {
|
||||
log_message(message, level)
|
||||
} else {
|
||||
if (level %in% c("ERROR", "WARNING")) {
|
||||
warning(message)
|
||||
} else {
|
||||
message(message)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#' Creates a sub-chunk for use within RMarkdown documents
|
||||
#'
|
||||
#' @param g A ggplot object to render in the sub-chunk
|
||||
#' @param fig_height Height of the figure in inches
|
||||
#' @param fig_width Width of the figure in inches
|
||||
#' @return NULL (writes chunk directly to output)
|
||||
#'
|
||||
subchunkify <- function(g, fig_height=7, fig_width=5) {
|
||||
g_deparsed <- paste0(deparse(
|
||||
function() {g}
|
||||
|
|
@ -12,14 +43,56 @@ subchunkify <- function(g, fig_height=7, fig_width=5) {
|
|||
")
|
||||
|
||||
cat(knitr::knit(text = knitr::knit_expand(text = sub_chunk), quiet = TRUE))
|
||||
}
|
||||
}
|
||||
|
||||
#' Creates a Chlorophyll Index map for a pivot
|
||||
#'
|
||||
#' @param pivot_raster The raster data containing CI values
|
||||
#' @param pivot_shape The shape of the pivot field
|
||||
#' @param pivot_spans Additional boundary data for the field
|
||||
#' @param show_legend Whether to show the legend (default: FALSE)
|
||||
#' @param legend_is_portrait Whether to show the legend in portrait orientation (default: FALSE)
|
||||
#' @param week Week number to display in the title
|
||||
#' @param age Age of the crop in weeks
|
||||
#' @param borders Whether to display field borders (default: FALSE)
|
||||
#' @return A tmap object with the CI map
|
||||
#'
|
||||
create_CI_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = F, legend_is_portrait = F, week, age, borders = FALSE){
|
||||
# Input validation
|
||||
if (missing(pivot_raster) || is.null(pivot_raster)) {
|
||||
stop("pivot_raster is required")
|
||||
}
|
||||
if (missing(pivot_shape) || is.null(pivot_shape)) {
|
||||
stop("pivot_shape is required")
|
||||
}
|
||||
if (missing(pivot_spans) || is.null(pivot_spans)) {
|
||||
stop("pivot_spans is required")
|
||||
}
|
||||
if (missing(week) || is.null(week)) {
|
||||
stop("week parameter is required")
|
||||
}
|
||||
if (missing(age) || is.null(age)) {
|
||||
stop("age parameter is required")
|
||||
}
|
||||
|
||||
# Create the base map
|
||||
map <- tm_shape(pivot_raster, unit = "m") +
|
||||
tm_raster(breaks = c(0,0.5,1,2,3,4,5,6,7,Inf), palette = "RdYlGn",legend.is.portrait = legend_is_portrait ,midpoint = NA) +
|
||||
tm_raster(breaks = c(0,0.5,1,2,3,4,5,6,7,Inf),
|
||||
palette = "RdYlGn",
|
||||
legend.is.portrait = legend_is_portrait,
|
||||
midpoint = NA,
|
||||
title = "CI") +
|
||||
tm_layout(main.title = paste0("\nMax CI week ", week,"\n", age, " weeks old"),
|
||||
main.title.size = 0.7, legend.show = show_legend)
|
||||
main.title.size = 0.7,
|
||||
legend.show = show_legend,
|
||||
legend.position = c("left", "bottom"),
|
||||
legend.width = 0.3,
|
||||
legend.height = 0.3,
|
||||
legend.text.size = 0.6,
|
||||
legend.title.size = 0.7,
|
||||
legend.outside = FALSE)
|
||||
|
||||
# Add borders if requested
|
||||
if (borders) {
|
||||
map <- map +
|
||||
tm_shape(pivot_shape) +
|
||||
|
|
@ -32,12 +105,55 @@ create_CI_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend =
|
|||
return(map)
|
||||
}
|
||||
|
||||
#' Creates a Chlorophyll Index difference map between two weeks
|
||||
#'
|
||||
#' @param pivot_raster The raster data containing CI difference values
|
||||
#' @param pivot_shape The shape of the pivot field
|
||||
#' @param pivot_spans Additional boundary data for the field
|
||||
#' @param show_legend Whether to show the legend (default: FALSE)
|
||||
#' @param legend_is_portrait Whether to show the legend in portrait orientation (default: FALSE)
|
||||
#' @param week_1 First week number for comparison
|
||||
#' @param week_2 Second week number for comparison
|
||||
#' @param age Age of the crop in weeks
|
||||
#' @param borders Whether to display field borders (default: TRUE)
|
||||
#' @return A tmap object with the CI difference map
|
||||
#'
|
||||
create_CI_diff_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = F, legend_is_portrait = F, week_1, week_2, age, borders = TRUE){
|
||||
map <- tm_shape(pivot_raster, unit = "m") +
|
||||
tm_raster(breaks = c(-3,-2,-1,0,1,2,3), palette = "RdYlGn",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 = 0.7, legend.show = show_legend)
|
||||
# Input validation
|
||||
if (missing(pivot_raster) || is.null(pivot_raster)) {
|
||||
stop("pivot_raster is required")
|
||||
}
|
||||
if (missing(pivot_shape) || is.null(pivot_shape)) {
|
||||
stop("pivot_shape is required")
|
||||
}
|
||||
if (missing(pivot_spans) || is.null(pivot_spans)) {
|
||||
stop("pivot_spans is required")
|
||||
}
|
||||
if (missing(week_1) || is.null(week_1) || missing(week_2) || is.null(week_2)) {
|
||||
stop("week_1 and week_2 parameters are required")
|
||||
}
|
||||
if (missing(age) || is.null(age)) {
|
||||
stop("age parameter is required")
|
||||
}
|
||||
|
||||
# Create the difference map
|
||||
map <- tm_shape(pivot_raster, unit = "m") +
|
||||
tm_raster(breaks = c(-3,-2,-1,0,1,2,3),
|
||||
palette = "RdYlGn",
|
||||
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 = 0.7,
|
||||
legend.show = show_legend,
|
||||
legend.position = c("left", "bottom"),
|
||||
legend.width = 0.3,
|
||||
legend.height = 0.3,
|
||||
legend.text.size = 0.6,
|
||||
legend.title.size = 0.7,
|
||||
legend.outside = FALSE)
|
||||
|
||||
# Add borders if requested
|
||||
if (borders) {
|
||||
map <- map +
|
||||
tm_shape(pivot_shape) +
|
||||
|
|
@ -50,165 +166,282 @@ create_CI_diff_map <- function(pivot_raster, pivot_shape, pivot_spans, show_lege
|
|||
return(map)
|
||||
}
|
||||
|
||||
#' Creates a visualization of CI data for a specific pivot field
|
||||
#'
|
||||
#' @param pivotName The name or ID of the pivot field to visualize
|
||||
#' @return NULL (adds output directly to R Markdown document)
|
||||
#'
|
||||
ci_plot <- function(pivotName){
|
||||
# pivotName = "1.1"
|
||||
pivotShape <- AllPivots0 %>% terra::subset(field %in% pivotName) %>% st_transform(crs(CI))
|
||||
age <- harvesting_data %>% dplyr::filter(field %in% pivotName) %>% sort("year") %>% tail(., 1) %>% dplyr::select(age) %>% unique() %>% pull() %>% round()
|
||||
|
||||
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)
|
||||
|
||||
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(season_start) %>% unique()
|
||||
|
||||
joined_spans2 <- AllPivots0 %>% st_transform(crs(pivotShape)) %>% dplyr::filter(field %in% pivotName) #%>% unique() %>% 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, borders = borders)
|
||||
CImap_m1 <- create_CI_map(singlePivot_m1, AllPivots2, joined_spans2, show_legend= F, legend_is_portrait = F, week = week_minus_1, age = age -1, borders = borders)
|
||||
CImap <- create_CI_map(singlePivot, AllPivots2, joined_spans2, show_legend= F, legend_is_portrait = F, week = week, age = age, borders = borders)
|
||||
|
||||
|
||||
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, borders = borders)
|
||||
CI_max_abs_three_week <- create_CI_diff_map(abs_CI_three_week, AllPivots2, joined_spans2, show_legend = T, legend_is_portrait = T, week_1 = week, week_2 = week_minus_3, age = age, borders = borders)
|
||||
|
||||
tst <- tmap_arrange(CImap_m2, CImap_m1, CImap,CI_max_abs_last_week, CI_max_abs_three_week, nrow = 1)
|
||||
|
||||
cat(paste("## Field", pivotName, "-", age, "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)
|
||||
|
||||
}
|
||||
|
||||
cum_ci_plot <- function(pivotName){
|
||||
|
||||
# pivotName = "3a13"
|
||||
data_ci <- CI_quadrant %>% filter(field == pivotName)
|
||||
|
||||
if (nrow(data_ci) == 0) {
|
||||
return(cum_ci_plot2(pivotName)) # Return an empty data frame if no data is found
|
||||
# Input validation
|
||||
if (missing(pivotName) || is.null(pivotName) || pivotName == "") {
|
||||
stop("pivotName is required")
|
||||
}
|
||||
if (!exists("AllPivots0") || !exists("CI") || !exists("CI_m1") || !exists("CI_m2")) {
|
||||
stop("Required global variables (AllPivots0, CI, CI_m1, CI_m2) not found")
|
||||
}
|
||||
if (!exists("harvesting_data")) {
|
||||
stop("harvesting_data not found")
|
||||
}
|
||||
|
||||
data_ci2 <- data_ci %>% mutate(CI_rate = cumulative_CI/DOY,
|
||||
week = week(Date))%>% group_by(field) %>%
|
||||
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)
|
||||
|
||||
|
||||
unique_seasons <- unique(date_preperation_perfect_pivot$season)
|
||||
|
||||
g <- ggplot(data= data_ci2 %>% filter(season %in% unique_seasons)) +
|
||||
facet_wrap(~season, scales = "free_x") +
|
||||
geom_line( aes(Date, mean_rolling10, col = sub_field, group = sub_field)) +
|
||||
labs(title = paste("14 day rolling MEAN CI rate - Pivot ", pivotName),
|
||||
color = "Field name")+
|
||||
scale_x_date(date_breaks = "1 month", date_labels = "%m-%Y") +
|
||||
theme_minimal() +
|
||||
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, 3.2, 10)
|
||||
|
||||
# Extract pivot shape and age data
|
||||
tryCatch({
|
||||
pivotShape <- AllPivots0 %>% terra::subset(field %in% pivotName) %>% sf::st_transform(terra::crs(CI))
|
||||
age <- harvesting_data %>%
|
||||
dplyr::filter(field %in% pivotName) %>%
|
||||
sort("year") %>%
|
||||
tail(., 1) %>%
|
||||
dplyr::select(age) %>%
|
||||
unique() %>%
|
||||
pull() %>%
|
||||
round()
|
||||
|
||||
# Filter for the specific pivot
|
||||
AllPivots2 <- AllPivots0 %>% dplyr::filter(field %in% pivotName)
|
||||
|
||||
# Create crop masks for different timepoints using terra functions
|
||||
singlePivot <- terra::crop(CI, pivotShape) %>% terra::mask(., pivotShape)
|
||||
singlePivot_m1 <- terra::crop(CI_m1, pivotShape) %>% terra::mask(., pivotShape)
|
||||
singlePivot_m2 <- terra::crop(CI_m2, pivotShape) %>% terra::mask(., pivotShape)
|
||||
|
||||
# Create difference maps
|
||||
abs_CI_last_week <- terra::crop(last_week_dif_raster_abs, pivotShape) %>% terra::mask(., pivotShape)
|
||||
abs_CI_three_week <- terra::crop(three_week_dif_raster_abs, pivotShape) %>% terra::mask(., pivotShape)
|
||||
|
||||
# Get planting date
|
||||
planting_date <- harvesting_data %>%
|
||||
dplyr::filter(field %in% pivotName) %>%
|
||||
ungroup() %>%
|
||||
dplyr::select(season_start) %>%
|
||||
unique()
|
||||
|
||||
# Create spans for borders
|
||||
joined_spans2 <- AllPivots0 %>%
|
||||
sf::st_transform(sf::st_crs(pivotShape)) %>%
|
||||
dplyr::filter(field %in% pivotName)
|
||||
|
||||
# Create the maps for different timepoints
|
||||
CImap_m2 <- create_CI_map(singlePivot_m2, AllPivots2, joined_spans2,
|
||||
show_legend = TRUE, legend_is_portrait = TRUE,
|
||||
week = week_minus_2, age = age - 2, borders = borders)
|
||||
|
||||
CImap_m1 <- create_CI_map(singlePivot_m1, AllPivots2, joined_spans2,
|
||||
show_legend = FALSE, legend_is_portrait = FALSE,
|
||||
week = week_minus_1, age = age - 1, borders = borders)
|
||||
|
||||
CImap <- create_CI_map(singlePivot, AllPivots2, joined_spans2,
|
||||
show_legend = FALSE, legend_is_portrait = FALSE,
|
||||
week = week, age = age, borders = borders)
|
||||
|
||||
# Create difference maps - only show legend on the second one to avoid redundancy
|
||||
CI_max_abs_last_week <- create_CI_diff_map(abs_CI_last_week, AllPivots2, joined_spans2,
|
||||
show_legend = FALSE, legend_is_portrait = TRUE,
|
||||
week_1 = week, week_2 = week_minus_1, age = age, borders = borders)
|
||||
|
||||
CI_max_abs_three_week <- create_CI_diff_map(abs_CI_three_week, AllPivots2, joined_spans2,
|
||||
show_legend = TRUE, legend_is_portrait = TRUE,
|
||||
week_1 = week, week_2 = week_minus_3, age = age, borders = borders)
|
||||
|
||||
# Arrange the maps
|
||||
tst <- tmap_arrange(CImap_m2, CImap_m1, CImap, CI_max_abs_last_week, CI_max_abs_three_week, nrow = 1)
|
||||
|
||||
# Output heading and map to R Markdown
|
||||
cat(paste("## Field", pivotName, "-", age, "weeks after planting/harvest", "\n"))
|
||||
print(tst)
|
||||
|
||||
}, error = function(e) {
|
||||
safe_log(paste("Error creating CI plot for pivot", pivotName, ":", e$message), "ERROR")
|
||||
cat(paste("## Field", pivotName, "- Error creating visualization", "\n"))
|
||||
cat(paste("Error:", e$message, "\n"))
|
||||
})
|
||||
}
|
||||
|
||||
cum_ci_plot <- function(pivotName, plot_type = "value", facet_on = TRUE) {
|
||||
|
||||
# pivotName = "3a13"
|
||||
data_ci <- CI_quadrant %>% filter(field == pivotName)
|
||||
|
||||
if (nrow(data_ci) == 0) {
|
||||
return(cum_ci_plot2(pivotName)) # Return an empty data frame if no data is found
|
||||
#' Creates a plot showing Chlorophyll Index data over time for a pivot field
|
||||
#'
|
||||
#' @param pivotName The name or ID of the pivot field to visualize
|
||||
#' @param plot_type Type of plot to generate: "value", "CI_rate", or "cumulative_CI"
|
||||
#' @param facet_on Whether to facet the plot by season (TRUE) or overlay all seasons (FALSE)
|
||||
#' @return NULL (adds output directly to R Markdown document)
|
||||
#'
|
||||
cum_ci_plot <- function(pivotName, plot_type = "value", facet_on = FALSE) {
|
||||
# Input validation
|
||||
if (missing(pivotName) || is.null(pivotName) || pivotName == "") {
|
||||
stop("pivotName is required")
|
||||
}
|
||||
if (!exists("CI_quadrant")) {
|
||||
stop("Required global variable CI_quadrant not found")
|
||||
}
|
||||
if (!plot_type %in% c("value", "CI_rate", "cumulative_CI")) {
|
||||
stop("plot_type must be one of: 'value', 'CI_rate', or 'cumulative_CI'")
|
||||
}
|
||||
|
||||
data_ci2 <- data_ci %>%
|
||||
mutate(CI_rate = cumulative_CI / DOY,
|
||||
week = week(Date)) %>%
|
||||
group_by(field) %>%
|
||||
mutate(mean_CIrate_rolling_10_days = rollapplyr(CI_rate, width = 10, FUN = mean, partial = TRUE),
|
||||
mean_rolling_10_days = rollapplyr(value, width = 10, FUN = mean, partial = TRUE))
|
||||
|
||||
data_ci2 <- data_ci2 %>% mutate(season = as.factor(season))
|
||||
|
||||
date_preperation_perfect_pivot <- data_ci2 %>%
|
||||
group_by(season) %>%
|
||||
summarise(min_date = min(Date),
|
||||
max_date = max(Date),
|
||||
days = max_date - min_date)
|
||||
|
||||
unique_seasons <- sort(unique(date_preperation_perfect_pivot$season), decreasing = TRUE)[1:3]
|
||||
|
||||
# Determine the y aesthetic based on the plot type
|
||||
y_aesthetic <- switch(plot_type,
|
||||
"CI_rate" = "mean_CIrate_rolling_10_days",
|
||||
"cumulative_CI" = "cumulative_CI",
|
||||
"value" = "mean_rolling_10_days")
|
||||
|
||||
y_label <- switch(plot_type,
|
||||
"CI_rate" = "10-Day Rolling Mean CI Rate (cumulative CI / age)",
|
||||
"cumulative_CI" = "Cumulative CI",
|
||||
"value" = "10-Day Rolling Mean CI")
|
||||
|
||||
if (facet_on) {
|
||||
g <- ggplot(data = data_ci2 %>% filter(season %in% unique_seasons)) +
|
||||
facet_wrap(~season, scales = "free_x") +
|
||||
geom_line(aes_string(x = "Date", y = y_aesthetic, col = "sub_field", group = "sub_field")) +
|
||||
labs(title = paste("Plot of", y_label, "for Pivot", pivotName),
|
||||
color = "Field Name",
|
||||
y = y_label) +
|
||||
scale_x_date(date_breaks = "1 month", date_labels = "%m-%Y") +
|
||||
theme_minimal() +
|
||||
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))
|
||||
} else {
|
||||
g <- ggplot(data = data_ci2 %>% filter(season %in% unique_seasons)) +
|
||||
geom_line(aes_string(x = "DOY", y = y_aesthetic, col = "season", group = "season")) +
|
||||
labs(title = paste("Plot of", y_label, "for Pivot", pivotName),
|
||||
color = "Season",
|
||||
y = y_label,
|
||||
x = "Age of Crop (Days)") +
|
||||
theme_minimal() +
|
||||
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, 3.2, 10)
|
||||
# Filter data for the specified pivot
|
||||
tryCatch({
|
||||
data_ci <- CI_quadrant %>% filter(field == pivotName)
|
||||
|
||||
if (nrow(data_ci) == 0) {
|
||||
safe_log(paste("No CI data found for pivot", pivotName), "WARNING")
|
||||
return(cum_ci_plot2(pivotName)) # Use fallback function when no data is available
|
||||
}
|
||||
|
||||
# Process data
|
||||
data_ci2 <- data_ci %>%
|
||||
mutate(CI_rate = cumulative_CI / DOY,
|
||||
week = week(Date)) %>%
|
||||
group_by(field) %>%
|
||||
mutate(mean_CIrate_rolling_10_days = rollapplyr(CI_rate, width = 10, FUN = mean, partial = TRUE),
|
||||
mean_rolling_10_days = rollapplyr(value, width = 10, FUN = mean, partial = TRUE))
|
||||
|
||||
data_ci2 <- data_ci2 %>% mutate(season = as.factor(season))
|
||||
|
||||
# Prepare date information by season
|
||||
date_preperation_perfect_pivot <- data_ci2 %>%
|
||||
group_by(season) %>%
|
||||
summarise(min_date = min(Date),
|
||||
max_date = max(Date),
|
||||
days = max_date - min_date)
|
||||
|
||||
# Get the 3 most recent seasons
|
||||
unique_seasons <- sort(unique(date_preperation_perfect_pivot$season), decreasing = TRUE)[1:3]
|
||||
|
||||
# Determine the y aesthetic based on the plot type
|
||||
y_aesthetic <- switch(plot_type,
|
||||
"CI_rate" = "mean_CIrate_rolling_10_days",
|
||||
"cumulative_CI" = "cumulative_CI",
|
||||
"value" = "mean_rolling_10_days")
|
||||
|
||||
y_label <- switch(plot_type,
|
||||
"CI_rate" = "10-Day Rolling Mean CI Rate (cumulative CI / age)",
|
||||
"cumulative_CI" = "Cumulative CI",
|
||||
"value" = "10-Day Rolling Mean CI")
|
||||
|
||||
# Create plot with either facets by season or overlay by DOY
|
||||
if (facet_on) {
|
||||
g <- ggplot(data = data_ci2 %>% filter(season %in% unique_seasons)) +
|
||||
facet_wrap(~season, scales = "free_x") +
|
||||
geom_line(aes_string(x = "Date", y = y_aesthetic, col = "sub_field", group = "sub_field")) +
|
||||
labs(title = paste("Plot of", y_label, "for Pivot", pivotName),
|
||||
color = "Field Name",
|
||||
y = y_label) +
|
||||
scale_x_date(date_breaks = "1 month", date_labels = "%m-%Y") +
|
||||
theme_minimal() +
|
||||
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))
|
||||
} else {
|
||||
g <- ggplot(data = data_ci2 %>% filter(season %in% unique_seasons)) +
|
||||
geom_line(aes_string(x = "DOY", y = y_aesthetic, col = "season", group = "season")) +
|
||||
labs(title = paste("Plot of", y_label, "for Pivot", pivotName),
|
||||
color = "Season",
|
||||
y = y_label,
|
||||
x = "Age of Crop (Days)") +
|
||||
theme_minimal() +
|
||||
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))
|
||||
}
|
||||
|
||||
# Output plot to R Markdown
|
||||
subchunkify(g, 3.2, 10)
|
||||
|
||||
}, error = function(e) {
|
||||
safe_log(paste("Error creating CI trend plot for pivot", pivotName, ":", e$message), "ERROR")
|
||||
cum_ci_plot2(pivotName) # Use fallback function in case of error
|
||||
})
|
||||
}
|
||||
|
||||
#' Fallback function for creating CI visualization when data is missing
|
||||
#'
|
||||
#' @param pivotName The name or ID of the pivot field to visualize
|
||||
#' @return NULL (adds output directly to R Markdown document)
|
||||
#'
|
||||
cum_ci_plot2 <- function(pivotName){
|
||||
end_date <- Sys.Date()
|
||||
start_date <- end_date %m-% months(11) # 11 months ago from end_date
|
||||
date_seq <- seq.Date(from = start_date, to = end_date, by = "month")
|
||||
midpoint_date <- start_date + (end_date - start_date) / 2
|
||||
g <- ggplot() +
|
||||
scale_x_date(limits = c(start_date, end_date), date_breaks = "1 month", date_labels = "%m-%Y") +
|
||||
scale_y_continuous(limits = c(0, 4)) +
|
||||
labs(title = paste("14 day rolling MEAN CI rate - Field ", pivotName),
|
||||
x = "Date", y = "CI Rate") +
|
||||
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)) +
|
||||
annotate("text", x = midpoint_date, y = 2, label = "No data available", size = 6, hjust = 0.5)
|
||||
# Input validation
|
||||
if (missing(pivotName) || is.null(pivotName) || pivotName == "") {
|
||||
stop("pivotName is required")
|
||||
}
|
||||
|
||||
subchunkify(g, 3.2, 10)
|
||||
}
|
||||
# Create a simple plot showing "No data available"
|
||||
tryCatch({
|
||||
end_date <- Sys.Date()
|
||||
start_date <- end_date %m-% months(11) # 11 months ago from end_date
|
||||
date_seq <- seq.Date(from = start_date, to = end_date, by = "month")
|
||||
midpoint_date <- start_date + (end_date - start_date) / 2
|
||||
|
||||
g <- ggplot() +
|
||||
scale_x_date(limits = c(start_date, end_date), date_breaks = "1 month", date_labels = "%m-%Y") +
|
||||
scale_y_continuous(limits = c(0, 4)) +
|
||||
labs(title = paste("14 day rolling MEAN CI rate - Field ", pivotName),
|
||||
x = "Date", y = "CI Rate") +
|
||||
theme_minimal() +
|
||||
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)) +
|
||||
annotate("text", x = midpoint_date, y = 2, label = "No data available", size = 6, hjust = 0.5)
|
||||
|
||||
subchunkify(g, 3.2, 10)
|
||||
|
||||
}, error = function(e) {
|
||||
safe_log(paste("Error creating fallback CI plot for pivot", pivotName, ":", e$message), "ERROR")
|
||||
cat(paste("No data available for field", pivotName, "\n"))
|
||||
})
|
||||
}
|
||||
|
||||
#' Gets the file path for a specific week's mosaic
|
||||
#'
|
||||
#' @param mosaic_path Base directory containing mosaic files
|
||||
#' @param input_date Reference date to calculate from
|
||||
#' @param week_offset Number of weeks to offset from input date (positive or negative)
|
||||
#' @return File path to the requested week's mosaic TIF file
|
||||
#'
|
||||
get_week_path <- function(mosaic_path, input_date, week_offset) {
|
||||
# Input validation
|
||||
if (missing(mosaic_path) || is.null(mosaic_path) || mosaic_path == "") {
|
||||
stop("mosaic_path is required")
|
||||
}
|
||||
if (missing(input_date)) {
|
||||
stop("input_date is required")
|
||||
}
|
||||
|
||||
tryCatch({
|
||||
# Convert input_date to Date object (in case it's a string)
|
||||
input_date <- as.Date(input_date)
|
||||
if (is.na(input_date)) {
|
||||
stop("Invalid input_date. Expected a Date object or a string convertible to Date.")
|
||||
}
|
||||
|
||||
# Validate week_offset
|
||||
week_offset <- as.integer(week_offset)
|
||||
if (is.na(week_offset)) {
|
||||
stop("Invalid week_offset. Expected an integer value.")
|
||||
}
|
||||
|
||||
# Get the start of the week for the input date (adjust to Monday as the start of the week)
|
||||
start_of_week <- lubridate::floor_date(input_date, unit = "week", week_start = 1)
|
||||
|
||||
# Calculate the new date after applying the week offset
|
||||
target_date <- start_of_week + lubridate::weeks(week_offset)
|
||||
|
||||
# Get the week number and year of the target date
|
||||
target_week <- sprintf("%02d", lubridate::isoweek(target_date)) # Left-pad week number with a zero if needed
|
||||
target_year <- lubridate::isoyear(target_date)
|
||||
|
||||
# Generate the file path for the target week
|
||||
path_to_week <- here::here(mosaic_path, paste0("week_", target_week, "_", target_year, ".tif"))
|
||||
|
||||
# Log the path calculation
|
||||
safe_log(paste("Calculated path for week", target_week, "of year", target_year, ":", path_to_week), "INFO")
|
||||
|
||||
# Return the path
|
||||
return(path_to_week)
|
||||
|
||||
}, error = function(e) {
|
||||
safe_log(paste("Error calculating week path:", e$message), "ERROR")
|
||||
stop(e$message)
|
||||
})
|
||||
}
|
||||
|
||||
|
|
|
|||
38
r_app/run_tests.R
Normal file
38
r_app/run_tests.R
Normal file
|
|
@ -0,0 +1,38 @@
|
|||
# run_tests.R
|
||||
#
|
||||
# TEST RUNNER FOR SMARTCANE
|
||||
# =======================
|
||||
# This script runs all tests for the SmartCane project.
|
||||
# Usage: Rscript run_tests.R [test_pattern]
|
||||
# - test_pattern: Optional regex pattern to match test files (default: all test_*.R files)
|
||||
#
|
||||
|
||||
# Process command line arguments
|
||||
args <- commandArgs(trailingOnly = TRUE)
|
||||
test_pattern <- if (length(args) > 0) args[1] else "^test_.+\\.R$"
|
||||
|
||||
# Set working directory to script location
|
||||
script_dir <- dirname(normalizePath(sys.frame(1)$ofile))
|
||||
setwd(script_dir)
|
||||
|
||||
# Source test framework
|
||||
source("tests/test_framework.R")
|
||||
|
||||
# Set up test environment
|
||||
env <- setup_test_env()
|
||||
|
||||
# Run tests
|
||||
cat("Running tests with pattern:", test_pattern, "\n\n")
|
||||
success <- run_tests(test_pattern)
|
||||
|
||||
# Clean up
|
||||
teardown_test_env()
|
||||
|
||||
# Exit with appropriate status code
|
||||
if (success) {
|
||||
cat("\n✓ All tests passed successfully!\n")
|
||||
quit(save = "no", status = 0)
|
||||
} else {
|
||||
cat("\n✗ Some tests failed.\n")
|
||||
quit(save = "no", status = 1)
|
||||
}
|
||||
68
r_app/tests/test_date_functions.R
Normal file
68
r_app/tests/test_date_functions.R
Normal file
|
|
@ -0,0 +1,68 @@
|
|||
# test_date_functions.R
|
||||
#
|
||||
# Tests for date-related functions in ci_extraction_utils.R
|
||||
#
|
||||
|
||||
# Load the test framework
|
||||
source("tests/test_framework.R")
|
||||
|
||||
# Set up test environment
|
||||
env <- setup_test_env()
|
||||
|
||||
# Load the functions to test
|
||||
source("../ci_extraction_utils.R")
|
||||
|
||||
# Test the date_list function
|
||||
test_that("date_list creates correct date sequences", {
|
||||
# Test with a specific date and offset
|
||||
dates <- date_list(as.Date("2023-01-15"), 7)
|
||||
|
||||
# Check the structure
|
||||
expect_type(dates, "list")
|
||||
expect_equal(names(dates), c("week", "year", "days_filter", "start_date", "end_date"))
|
||||
|
||||
# Check the values
|
||||
expect_equal(dates$week, lubridate::week(as.Date("2023-01-09")))
|
||||
expect_equal(dates$year, 2023)
|
||||
expect_equal(dates$start_date, as.Date("2023-01-09"))
|
||||
expect_equal(dates$end_date, as.Date("2023-01-15"))
|
||||
expect_equal(length(dates$days_filter), 7)
|
||||
expect_equal(dates$days_filter[1], "2023-01-09")
|
||||
expect_equal(dates$days_filter[7], "2023-01-15")
|
||||
|
||||
# Test with a different offset
|
||||
dates_short <- date_list(as.Date("2023-01-15"), 3)
|
||||
expect_equal(length(dates_short$days_filter), 3)
|
||||
expect_equal(dates_short$days_filter, c("2023-01-13", "2023-01-14", "2023-01-15"))
|
||||
|
||||
# Test with string date
|
||||
dates_string <- date_list("2023-01-15", 5)
|
||||
expect_equal(dates_string$days_filter,
|
||||
c("2023-01-11", "2023-01-12", "2023-01-13", "2023-01-14", "2023-01-15"))
|
||||
|
||||
# Test error handling
|
||||
expect_error(date_list("invalid-date", 7),
|
||||
"Invalid end_date provided")
|
||||
expect_error(date_list("2023-01-15", -1),
|
||||
"Invalid offset provided")
|
||||
})
|
||||
|
||||
# Test the date_extract function
|
||||
test_that("date_extract correctly extracts dates from file paths", {
|
||||
# Test with various file path formats
|
||||
expect_equal(date_extract("/some/path/2023-01-15_image.tif"), "2023-01-15")
|
||||
expect_equal(date_extract("path/to/planet_2023-01-15.tif"), "2023-01-15")
|
||||
expect_equal(date_extract("c:\\path\\with\\windows\\2023-01-15_file.tif"), "2023-01-15")
|
||||
expect_equal(date_extract("2023-01-15.tif"), "2023-01-15")
|
||||
expect_equal(date_extract("prefix-2023-01-15-suffix.tif"), "2023-01-15")
|
||||
|
||||
# Test with invalid file paths
|
||||
expect_warning(result <- date_extract("no-date-here.tif"), "Could not extract date")
|
||||
expect_true(is.na(result))
|
||||
})
|
||||
|
||||
# Clean up
|
||||
teardown_test_env()
|
||||
|
||||
# Print success message
|
||||
cat("Date function tests completed successfully\n")
|
||||
120
r_app/tests/test_framework.R
Normal file
120
r_app/tests/test_framework.R
Normal file
|
|
@ -0,0 +1,120 @@
|
|||
# test_framework.R
|
||||
#
|
||||
# TEST FRAMEWORK FOR SMARTCANE
|
||||
# ===========================
|
||||
# This script provides a simple testing framework for the SmartCane project.
|
||||
# It includes utilities for setting up test environments and running tests.
|
||||
#
|
||||
|
||||
# Install required packages if not available
|
||||
if (!require("testthat", quietly = TRUE)) {
|
||||
install.packages("testthat", repos = "https://cran.rstudio.com/")
|
||||
}
|
||||
library(testthat)
|
||||
|
||||
# Define paths for testing
|
||||
test_root <- file.path(normalizePath(".."), "tests")
|
||||
test_data_dir <- file.path(test_root, "test_data")
|
||||
|
||||
# Create test directories if they don't exist
|
||||
dir.create(test_data_dir, recursive = TRUE, showWarnings = FALSE)
|
||||
|
||||
# Set up a test environment with all necessary data
|
||||
setup_test_env <- function() {
|
||||
# Add working directory to the path
|
||||
.libPaths(c(.libPaths(), normalizePath("..")))
|
||||
|
||||
# Source required files with minimal dependencies
|
||||
tryCatch({
|
||||
source(file.path(normalizePath(".."), "packages.R"))
|
||||
skip_package_loading <- TRUE
|
||||
|
||||
# Load minimal dependencies for tests
|
||||
required_packages <- c("lubridate", "stringr", "purrr", "dplyr", "testthat")
|
||||
for (pkg in required_packages) {
|
||||
if (!require(pkg, character.only = TRUE, quietly = TRUE)) {
|
||||
warning(paste("Package", pkg, "not available, some tests may fail"))
|
||||
}
|
||||
}
|
||||
}, error = function(e) {
|
||||
warning("Error loading dependencies: ", e$message)
|
||||
})
|
||||
|
||||
# Set up test logging
|
||||
assign("log_message", function(message, level = "INFO") {
|
||||
cat(paste0("[", level, "] ", message, "\n"))
|
||||
}, envir = .GlobalEnv)
|
||||
|
||||
# Create a mock project structure
|
||||
test_project <- list(
|
||||
project_dir = "test_project",
|
||||
data_dir = test_data_dir,
|
||||
daily_CI_vals_dir = file.path(test_data_dir, "extracted_ci", "daily_vals"),
|
||||
cumulative_CI_vals_dir = file.path(test_data_dir, "extracted_ci", "cumulative_vals"),
|
||||
merged_final = file.path(test_data_dir, "merged_final"),
|
||||
daily_vrt = file.path(test_data_dir, "daily_vrt")
|
||||
)
|
||||
|
||||
# Create the directories
|
||||
for (dir in test_project) {
|
||||
if (is.character(dir)) {
|
||||
dir.create(dir, recursive = TRUE, showWarnings = FALSE)
|
||||
}
|
||||
}
|
||||
|
||||
return(test_project)
|
||||
}
|
||||
|
||||
# Clean up test environment
|
||||
teardown_test_env <- function() {
|
||||
# Clean up only test-created files if needed
|
||||
# We'll leave the main directories for inspection
|
||||
}
|
||||
|
||||
# Run all tests in a directory
|
||||
run_tests <- function(pattern = "^test_.+\\.R$") {
|
||||
test_files <- list.files(
|
||||
path = test_root,
|
||||
pattern = pattern,
|
||||
full.names = TRUE
|
||||
)
|
||||
|
||||
# Exclude this file
|
||||
test_files <- test_files[!grepl("test_framework\\.R$", test_files)]
|
||||
|
||||
if (length(test_files) == 0) {
|
||||
cat("No test files found matching pattern:", pattern, "\n")
|
||||
return(FALSE)
|
||||
}
|
||||
|
||||
cat("Found", length(test_files), "test files:\n")
|
||||
cat(paste(" -", basename(test_files)), sep = "\n")
|
||||
cat("\n")
|
||||
|
||||
# Run each test file
|
||||
results <- lapply(test_files, function(file) {
|
||||
cat("Running tests in:", basename(file), "\n")
|
||||
tryCatch({
|
||||
source(file, local = TRUE)
|
||||
cat("✓ Tests completed\n\n")
|
||||
TRUE
|
||||
}, error = function(e) {
|
||||
cat("✗ Error:", e$message, "\n\n")
|
||||
FALSE
|
||||
})
|
||||
})
|
||||
|
||||
# Summary
|
||||
success_count <- sum(unlist(results))
|
||||
cat("\nTest Summary:", success_count, "of", length(test_files),
|
||||
"test files completed successfully\n")
|
||||
|
||||
return(all(unlist(results)))
|
||||
}
|
||||
|
||||
# If this script is run directly, run all tests
|
||||
if (!interactive() && (basename(sys.frame(1)$ofile) == "test_framework.R")) {
|
||||
setup_test_env()
|
||||
run_tests()
|
||||
teardown_test_env()
|
||||
}
|
||||
Loading…
Reference in a new issue