Getting-started

Setup and Requirements

library(cercospoRa)
library(data.table)
library(terra)
#> terra 1.8.42
#> 
#> Attaching package: 'terra'
#> The following object is masked from 'package:data.table':
#> 
#>     shift
library(sf)
#> Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE
library(ggplot2)

The cercospoRa package estimates the daily infection values which tracks the progress towards a Cercospora beticola epidemic on sugar beet. The package requires formatted hourly weather data in a format which the model recognises.

Format weather station data

Included in the package is a weather data set from a station adjacent to a sugar beet field trial.

The weather data is first cleaned into to a format that the model can use without returning errors due to inconsistencies and NAs, which are common in weather data. This will require the user to think about how to replace the NAs or state default values to complete the data set.

# classify to data.table
wthr <- data.table(weathr)

# Use POSIXct formatted time.
wthr[,Time := as.POSIXct(paste0(Datum, " ",Stunde,":00"),tz = "UTC")]

# Nominate Latitude and Longitude location of the weather station. 
# While not needed in cercospoRa some plant disease models will use location to 
#  decide the closest weather station to pull weather from
wthr[, c("lon","lat") := list(9.916,51.41866)]

# weather is hourly and will error if we don't specify a standard deviation of 
#  weather direction. This is intentional to force the user to decide how variable
#  the wind direction data could be.
wthr[, wd_std := 20]

# remove all data after September as it contains missing data
wthr <- wthr[Datum < as.POSIXct("2022-10-01")]

# set NA wind speed values to zero
wthr[is.na(WG200), WG200 := 0]

Now we can use the format_weather() function.

wthr <- format_weather(wthr,
                         POSIXct_time = "Time",
                         time_zone = "UTC",
                         temp = "T200",
                         rain = "N100",
                         rh = "F200",
                         wd = "WR200",
                         ws = "WG200",
                         station = "Station",
                         lon = "lon",
                         lat = "lat",
                         wd_sd = "wd_std",
                         data_check = FALSE # this stops the function from checking for faults
                         )
# As the data is formatted closely enough for what is expected for the model. 
# We can elect to turn the data_check off so

Calculate epidemic onset

The following function calculates the earliest date from which a cercospora leaf spot epidemic could commence. start and end provide the time window for which the model should run. Canopy closure (c_closure) indicates when the model should start point for the model. This is due to an observation by Wolf canopy closure pr

This date is determined by the negative prognosis models to be the last day for which a epidemic is likely not to start.

cercospoRa::calc_epidemic_onset(start = as.POSIXct("2022-04-25",tz = "UTC"),
                    end = as.POSIXct("2022-09-30",tz = "UTC"),
                    c_closure = as.POSIXct("2022-07-30",tz = "UTC"),
                    weather = wthr,
                    cultivar_sus = 6)
#> [1] "2022-08-17 UTC"

Estimate canopy closure

Canopy closure is traditionally observed manually. The following methods permits canopy closure date estimation from remotely sensed data. UAV or satellite data can be used to produce georeferenced leaf area index (LAI) rasters.

# Get file location of example rasters with LAI values
image_files <- list.files(system.file("extdata", "uav_img",package = "cercospoRa"),
                          pattern = "tif",
                          full.names = TRUE)

# Read in data and check for consistency  
epidemic_onset_param <-
  read_sb_growth_parameter(img_files = image_files,
                           img_dates = as.POSIXct(c("2022-06-14","2022-06-28"),
                                                  tz = "UTC"),
                           target_res = 10)

terra::plot(epidemic_onset_param)

Calculate growth rate

calc_r_x0 calculates the sugar beet growth rate based on a minimum of two images.

param_rxt <- calc_r_x0(epidemic_onset_param,
                      min_r = 0.02,
                      max_r = 0.05,
                      k = 6)

Calculate canopy closure dates

Calculate canopy closure uses a logistic regression to estimate the canopy closure date for each raster cell.

canopy_closure <- calc_c_closure(param_rxt,
                                 x1 = 1.3,
                                 k=6)
terra::plot(canopy_closure)

Finally we can use a specialised wrapper function to calculate the earliest estimated cercospora leaf spot epidemic onset dates when chemical intervention might be needed.

epidemic_onset_map <- 
  calc_epidemic_onset_from_image(start =as.POSIXct("2022-04-25",tz = "UTC"),
                                 end = as.POSIXct("2022-09-30",tz = "UTC"),
                                 cc_r = canopy_closure,
                                 weather = wthr,
                                 cultivar_sus = 6)
epidemic_onset_map
#> class       : SpatRaster 
#> dimensions  : 17, 29, 1  (nrow, ncol, nlyr)
#> resolution  : 0.0001103147, 6.885404e-05  (x, y)
#> extent      : 9.915826, 9.919025, 51.49783, 51.499  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory
#> name        :      lyr.1 
#> min value   : 1653264000 
#> max value   : 1659571200

Plot the heterogeneity of onset date

terra::plot(epidemic_onset_map)

Convert numeric date back to POSIXct format

onset_dates <- as.POSIXct(terra::values(epidemic_onset_map),
                          tz = "UTC",
                          origin = "1970-01-01")

Plot date values

# convert to a data.table
dtdate <- data.table(pixel = seq_along(onset_dates),
                     ODates = as.IDate(onset_dates))

# Summarise the range of epidemic onset dates
summary(onset_dates[is.na(onset_dates)== FALSE])
#>                  Min.               1st Qu.                Median 
#> "2022-05-23 00:00:00" "2022-06-19 12:00:00" "2022-07-07 00:00:00" 
#>                  Mean               3rd Qu.                  Max. 
#> "2022-07-02 16:03:58" "2022-07-17 00:00:00" "2022-08-04 00:00:00"

# Plot dates on boxplot
dtdate[is.na(ODates) == FALSE,] |>
ggplot(aes(y=ODates))+
  geom_boxplot()