moving window with flexible function

[spono]

Ciao JR, I was thinking that a possible idea to enhance find_trees would be giving it the possibility to accept more functions (or a set of parameters) according to a second variable provided along with the height from CHM.

Basically, if I provide a stack with "height" + "forest_type" (conifer/broadleaf, high forest/coppice, etc.), the moving window applies a function according to a rule based on the ancillary layer. This would allow to be a bit more specific already during the process instead of -e.g.- running twice the algorithm for the different areas.

I guess it could be something like this:

params <- data.frame(a = c(0.05, 0.07), b = c(1.5, 2))

f <- function(x,y,z) { x * y + z }

find_trees(chm, layer2, lmf(f(a,b), params, hmin=2, shape = "circular"))

but I don't know how complex this would be to implement on the actual tools.

[Jean-Romain]

Could you be more specific? It sounds like a interesting request but I did not understand the details. Why f has 3 parameters but only 2 are used in f(a,b)? How values a and b are selected? How does the "forest type" is evaluated at run time? Using merge_spatial and a shapefile (layer2)? In this case I could design something more generic like:

f <- function(x,y,z) { x * y + z }
find_trees(chm, lmf(f(Z, Intensity, ScanAngleRank)))

Obviously Z * Intensity + ScanAngle cannot give a good windows but something like that becomes possible

f <- function(x,y) { ifelse(y == 1, 0.05*x + 1.5, 0.07*x + 2) }
las <- merge_spatial(las, layer, "foresttype")
find_trees(chm, lmf(f(Z, foresttype)))

[spono]

Why f has 3 parameters but only 2 are used in f(a,b)?

it was just a draft, I didn't really mean to write a working example, but yeah: you definitely got the point in your second example! ...and passing through merge_spatial could be an idea. With a function designed in that way, I guess it might only need a method to operate on a stack with named layers and should be fine

How does the "forest type" is evaluated at run time?

not at run time: right now I see it only as a data fusion solution, i.e. the ancillary layer could come from an aerial/satellite image previously classified or from polygons of a management plan, etc etc.

[Jean-Romain]

it was just a draft, I didn't really mean to write a working example

Certainly! I understand that. But a syntactically valid example helps for the understanding :wink:

Below a quick and dirty proof of concept tweaking the code of lmf. It does what I described above. It works with v4 because I used the function plugin_itd() which does not exist in previous version

lmf_extended = function(ws, hmin = 2, shape = c("circular", "square"), ws_args = NULL)
{
  shape <- match.arg(shape)
  circ  <- shape == "circular"
  ws    <- lazyeval::uq(ws)
  hmin  <- lazyeval::uq(hmin)
  ws_args  <- lazyeval::uq(ws_args)

  f = function(las)
  {
    if (is.function(ws))
    {
      if (!is.null(ws_args))
      {
        args <- lapply(ws_args, function(x) if (x %in% names(las)) las@data[[x]] else x)
        ws <- do.call(ws, args)
        b <- las$Z < hmin
        ws[b] <- min(ws)
      }
      else
      {
        ws    <- ws(las@data$Z)
        b <- las$Z < hmin
        ws[b] <- ws(hmin)
      }

      n <- npoints(las)
      if (!is.numeric(ws)) stop("The function 'ws' did not return a correct output. ", call. = FALSE)
      if (any(ws <= 0))    stop("The function 'ws' returned negative or null values.", call. = FALSE)
      if (anyNA(ws))       stop("The function 'ws' returned NA values.",               call. = FALSE)
      if (length(ws) != n) stop("The function 'ws' did not return a correct output.",  call. = FALSE)

    }
    else if (!is.numeric(ws))
    {
      stop("'ws' must be a number or a function", call. = FALSE)
    }

    lidR:::force_autoindex(las) <- lidR:::LIDRGRIDPARTITION
    return(lidR:::C_lmf(las, ws, hmin, circ, lidR:::getThread()))
  }

  f <- plugin_itd(f , TRUE)
  return(f)
}

LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile)
ws <- function(x,y) { 0.1*x + 0.01*y }

algo <- lmf_extended(ws,  ws_args = list(x = "Z", y = "Intensity"))
u = locate_trees(las, algo)               
plot(sf::st_geometry(u))

algo <- lmf_extended(ws,  ws_args = list(x = "Z", y = 100))
u = locate_trees(las, algo)               
plot(sf::st_geometry(u))

algo
#> Object of class lidR algorithm
#> Algorithm for: individual tree detection 
#> Designed to be used with: locate_trees 
#> Native C++ parallelization: yes 
#> Parameters: 
#>  - circ = TRUE <logical>
#>  - hmin = 2 <numeric>
#>  - shape = circular <character>
#>  - ws <function>
#>  - ws_args <list>

[spono]

whoa, amazing! unfortunately I won't be able to test it soon but I crave for it. I'll give you some feedback asap. In the meantime, thank you very much, as usual

[Jean-Romain]

Added in v4.0.0 with minimal modification + backward compatibility