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

In OtsuFire: Fire Scars, Severity and Regeneration Mapping Using 'Otsu' Thresholding

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.

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'.

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.

Note

Examples require large external raster files (hosted on Zenodo) and depend on external software (Python, GDAL). Therefore, they are wrapped in dontrun to avoid errors during R CMD check and to ensure portability.

The DOY raster band should be derived from a Landsat composite generated in Google Earth Engine (GEE), following the approach described by Quintero et al. (2025).

References

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

Examples

## Not run: 
# 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
)

## End(Not run)

OtsuFire documentation built on June 13, 2025, 3:01 a.m.