OtsuFire (version 0.1.4)

calculate_doy_flags: Calculate DOY-Based Flags and Summary Statistics from Raster and Fire Polygons

Description

This function extracts Day-of-Year (DOY) values from a raster (e.g., from a Landsat composite), masked by burned area polygons, and computes per-polygon statistics: mode, median, and optionally percentiles (e.g., 5th, 25th, 95th). These DOY values are converted to calendar dates, and their components (day, month, year) are also included.

Optionally, the function can threshold the DOY band around the median and/or mode of each polygon using user-defined DOY windows (e.g., \pm10 days), create binary raster masks, and vectorize them using GDAL.

Use this function to assign a representative burning date (DOY) to each fire event and to discard events with unusually high internal variation in DOY among their pixels.

The function supports either a single shapefile (original burned areas) or a named list of shapefiles or `sf` objects generated with different Otsu thresholds.

Value

A named list with two elements:

sf_objects

Named list of `sf` outputs containing DOY statistics and date components for each input polygon set.

shp_paths

Named list of output shapefile paths, including DOY summary shapefiles and any polygonized threshold layers.

Arguments

raster

A multi-band `SpatRaster` object. One of the bands must contain DOY values.

doy_band

Integer. The index of the DOY band in the raster (default is 2).

polygons_sf

A single shapefile path, an `sf` object, or a named list of shapefiles or `sf` objects. Names are used to label outputs (e.g., `"ge50"`, `"ge100"`).

output_dir

Directory where output raster and shapefile files will be saved.

year

Calendar year (numeric or character) used for DOY-to-date conversion and filenames.

doy_thresholds

Numeric vector of DOY windows (e.g., `c(10, 15)`) to apply around the selected statistic(s).

stats

Statistic(s) to use for thresholding. One of `"median"`, `"mode"`, or `"both"` (default `"both"`).

calc_percentiles

Logical. If `TRUE`, compute additional DOY percentiles. Default is `TRUE`.

percentiles

Numeric vector of percentiles to calculate (default: `c(0.05, 0.25, 0.95)`).

polygonize

Logical. If `TRUE`, vectorize the DOY binary masks using GDAL (`gdal_polygonize.py`). Default is `TRUE`.

python_exe

Path to the Python executable. Required if `polygonize = TRUE`.

gdal_polygonize_script

Path to `gdal_polygonize.py`. Required if `polygonize = TRUE`.

get_mode

Optional custom function to calculate the mode. Default function provided.

get_median

Optional custom function to calculate the median. Default function provided.

get_quantile

Optional custom function to calculate percentiles. Default function provided.

keep_all_polygons

Logical. If `TRUE`, all polygons are kept and a flag column is added indicating which polygons meet the threshold condition. If `FALSE`, only polygons meeting the condition are saved. Default: `TRUE`.

References

Quintero, N., Viedma, O., Veraverbeke, S., & Moreno, J. M. (2025). *Optimising Regional Fire Severity Mapping Using Pixel-Based Image Compositing*. SSRN 4929831.

Examples

if (FALSE) {
# Case 1: Multiple Otsu-based burned area shapefiles
rast <- terra::rast("rbr_doy_stack_1987.tif")
polys <- list(
  ge50 = "burned_areas_1987_otsu_ge50.shp",
  ge100 = "burned_areas_1987_otsu_ge100.shp"
)
calculate_doy_flags(
  raster = rast,
  doy_band = 2,
  polygons_sf = polys,
  output_dir = "output/",
  year = 1987,
  stats = "both",
  doy_thresholds = c(10, 15),
  calc_percentiles = TRUE,
  percentiles = c(0.05, 0.25, 0.95),
  polygonize = TRUE,
  python_exe = "/usr/bin/python",
  gdal_polygonize_script = "/usr/bin/gdal_polygonize.py"
)

# Case 2: Single shapefile without polygonization
calculate_doy_flags(
  raster = rast,
  polygons_sf = "burned_areas_1987_origval.shp",
  output_dir = "output/",
  year = 1987,
  stats = "median",
  doy_thresholds = 10,
  calc_percentiles = FALSE,
  polygonize = FALSE
)
}