Can't calculate EVI2 from Sentinel-2 raster

[evhersh]

Hello,

Thanks for the awesome package, I've been using it a ton for my work and it's helped exponentially increase my productivity!

I have an issue relating to calculating indices that require additional parameters, in this case EVI2 from Sentinel-2. I think the problem is that I can't seem to specify the 'g' and 'L' parameters required for the function - perhaps I am not understanding how the 'extra_objects' parameter is supposed to work? But it also gives the same error when I replace the 'g' and 'L' parameters and change the formula itself. I haven't had issues with the other indices I've tried to calculate (NDVI, BI, OSAVI).

I've pasted some code that I tried below:

> calculate_indices(
+   raster = xg.s2.rast,
+   indices = asi[asi$short_name == "EVI2", ],
+   output_filename = here("data", "rasters", "bigrast", "layers", "EVI2.tif"),
+   extra_objects = list(g=2.5, L=1)
+ )

Error in `calculate_indices()`:
! Some indices cannot be calculated using the available image bands.
Run `rlang::last_trace()` to see where the error occurred.

> names(xg.s2.rast)
 [1] "A"   "B"   "G"   "R"   "RE1" "RE2" "RE3" "N"   "N2"  "WV"  "S1"  "S2" 

> asi[asi$short_name == "EVI2", ]
# A tibble: 1 × 9
  application_domain bands     contributor                date_of_addition formula                         long_name                          platforms reference short_name
  <chr>              <list>    <chr>                      <chr>            <chr>                           <chr>                              <list>    <chr>     <chr>     
1 vegetation         <chr [4]> https://github.com/davemlz 2021-04-07       g * (N - R) / (N + 2.4 * R + L) Two-Band Enhanced Vegetation Index <chr [6]> https://… EVI2      
> asi[asi$short_name == "EVI2", ]$platforms
[[1]]
[1] "Sentinel-2"    "Landsat-OLI"   "Landsat-TM"    "Landsat-ETM+"  "MODIS"         "Planet-Fusion"

# this results in the same error
asi[asi$short_name == "EVI2", ]$formula <- "2.5 * (N - R) / (N + 2.4 * R + 1)"

calculate_indices(
  raster = xg.s2.rast,
  indices = asi[asi$short_name == "EVI2", ],
  output_filename = here("data", "rasters", "bigrast", "layers", "EVI2.tif")
)

Error in `calculate_indices()`:
! Some indices cannot be calculated using the available image bands.

[mikemahoney218]

Thanks for the bug report! I took a few notes while scoping this out, using the self-contained example:

library(sf)

aoi <- spatialsample::boston_canopy |>
  sf::st_union()

sentinel2 <- get_sentinel2_imagery(
  aoi, 
  "2024-06-01", 
  "2024-09-01"
)

evi <- spectral_indices()
evi <- evi[evi$short_name == "EVI2", ]

try(
  calculate_indices(
    raster = sentinel2,
    indices = evi,
    output_filename = tempfile(fileext = ".tif")
  )
)

evi$formula <- "2.5 * (N - R) / (N + 2.4 * R + 1)"
try(
  calculate_indices(
    raster = sentinel2,
    indices = evi,
    output_filename = tempfile(fileext = ".tif")
  )
)

This bug shows up because this function naively checks the bands column. That's not documented anywhere and makes it challenging to provide custom functions to this...

rsi/R/calculate_indices.R

Lines 201 to 222 in 2707b7e

	check_indices <- function(
	remap_band_names,
	indices,
	call = rlang::caller_env()
	) {
	good <- TRUE
	if (!is.null(remap_band_names)) {
	good <- vapply(
	indices$bands,
	function(bands) all(bands %in% unlist(remap_band_names)),
	logical(1)
	)
	}
	if (!any(good)) {
	rlang::abort(
	"Some indices cannot be calculated using the available image bands.",
	call = call,
	class = "rsi_missing_indices"
	)
	}
	}

I'll look into that. For now, as a workaround, this works:

evi$bands <- list("R")
try(
  calculate_indices(
    raster = sentinel2,
    indices = evi,
    output_filename = tempfile(fileext = ".tif")
  )
)

That's definitely not ideal, though...

[mikemahoney218]

I guess this is kinda documented:

indices
A data frame of indices to compute. The intent is for this function to work with subsets of spectral_indices, but any data frame with columns formula (containing a string representation of the equation used to calculate the index), bands (a list column containing character vectors of the necessary bands) and short_name (which will be used as the band name) will work.

But I don't mention why you need the bands column. The fix might be to document that better...

Thinking out loud: The reason to throw this error is (1) so you don't waste time calculating 1000s of indices before getting a failure (2) because the failure that gets thrown is along the lines of:

Error in eval(calc) : object 'g' not found

Which isn't super obvious. So probably we want to keep this check intact...

[mikemahoney218]

My thinking is that I'm going to document this better (#95) but leave the behavior as-is. The documentation now says:

#' @section Custom indices:
#' While this function is designed to work with subsets of [spectral_indices()],
#' you can also provide your own custom indices by providing a data frame with
#' at least three columns:
#'
#' * `formula` should contain a string representation of the equation used to
#' calculate the index. This string will be interpreted as R code, and should use
#' normal R arithmetic operators. Any additional functions used in calculating the
#' index will need to be passed through `extra_objects`.
#' * `bands` should be a list column containing character vectors identifying the
#' bands required to calculate the index. These vectors will be checked against the
#' band names of `raster` to make sure that all necessary bands are available
#' before indices are calculated.
#' * `short_name` should contain a string which will be used as the band name in
#' the final output file.
#'
#' Sometimes you might want to tweak the formula used for an existing index, for
#' instance to set non-data parameters to constant values. In that case, make sure
#' you also update `bands` to reflect just the bands required by your updated formula,
#' so that the function won't error unnecessarily.

Hopefully that explains exactly what happened here -- the function still thought g and L were required, because the bands column still advertised them -- and helps others avoid it in the future.

Thanks again for the report!

[github-actions]

This issue has been automatically locked. If you believe you have found a related problem, please file a new issue (with a reprex: https://reprex.tidyverse.org) and link to this issue.