Focal filter with stars

[michaeldorman]

This is a general question - hope it's appropriate to ask here.

As someone who just started to work with stars (as replacement for raster), I find the raster::focal function the most significant thing missing from stars at the moment.

According to the migration table (#122), raster::focal is marked as "Not present, could be implemented".

I was wondering whether there are any more concrete development paths in mind -
Is it likely that that a focal filter function will be implemented in stars? If not, do you think it makes sense to use a wrapper that transforms stars<->raster and uses raster::focal? Perhaps you are aware of any general focal filtering function in R, which is aimed at array object and thus may work on stars objects as is?

Thanks for the great work - so far enjoyed using stars very much!

[edzer]

After applying this patch, you can try this:

library(stars)
tif = system.file("tif/L7_ETMs.tif", package = "stars")
x1 = read_stars(tif)
foc = function(x, w) {
	raster::as.matrix(raster::focal(raster::raster(x), w))
}
f = st_apply(x1, 3, foc, w = matrix(1, 3, 3))
plot(f)

[michaeldorman]

Nice! I will use this approach, thank you very much for the help and quick reply!

[florisvdh]

I think that the raster::raster(x) part of the above solution should be as(x, "Raster"), i.e. the function provided by stars. It seems that raster::raster() does not yet accept stars objects.

[edzer]

Did you try? Functions passed as argument to st_apply receive sub-arrays, in this case matrices, not (subsetted) stars objects.

[florisvdh]

You are right; my mistake! Thanks for the clarification.

[jan-abel-inwt]

Hi there, i hope this the right place. The Discussion above was helping me a lot. Thanks!

I did the same as described above using terra::focal instead of raster::focal and it's working fine so far.

apply_focal_filter <- function(stars_slice, w) {
  # Convert the stars slice to a terra raster
  raster <- rast(stars_slice)
  # Apply the filter
  filtered_raster <- focal(
    raster,
    w,
    fun = "sum",
    expand = TRUE
  )
  mat <- as.matrix(filtered_raster)
  # mat <- transpose_flat_mat(mat)
  return(mat)
}

sf <- st_apply(
    s,
    MARGIN = "date_time", # Specifies the dimension over which to apply function
    FUN = apply_focal_filter, # The filter function wrapper for terra::focal
    w = gauss_kernel, # The Weight matrix for the filter
    CLUSTER = my.cluster, # Cluster for parallel processing
    PROGRESS = TRUE # Display progress
  )
  # Rename Attribute / Feature
  sf <- setNames(sf, "kfz_per_hour_neigh")
  # Save the filtered result into the results list
  s_filtered <- sf
}

However, in my particular case, my data was transposed (mirrored at the diagonal) by the as.matrix() function. I believe this issue might be related to the orientation of the delta of my stars object. I am generating the stars object from a dataframe with a x, y, date_time and a attribute column:

dims <- c("x", "y", "date_time")
crs <- 25833

s <- st_as_stars(input, dims = dims)
st_crs(s) <- st_crs(crs)
s

stars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
              Min. 1st Qu. Median     Mean 3rd Qu. Max.
kfz_per_hour     0       0      0 201.7703     287 2987
dimension(s):
          from   to         offset   delta                refsys x/y
x            1   21         387575      50 ETRS89 / UTM zone 33N [x]
y            1   21        5823175     -50 ETRS89 / UTM zone 33N [y]
date_time    1 8760 2022-01-01 CET 1 hours               POSIXct    

I found a workaround by transposing the matrix after applying the filter with this function:

transpose_flat_mat <- function(mat) {
  # get dimension & name
  dim <- dim(mat)
  sq_dim <- sqrt(dim[1])
  name <- colnames(mat)

  # transform to square matrix & transpose
  dim(mat) <- c(sq_dim, sq_dim)
  mat <- t(mat)

  # set back to orig dimensions
  dim(mat) <- dim
  colnames(mat) <- name
  return(matrix)
}

Is there a more elegant solution to this issue?

[edzer]

Reading the as.matrix() help of package terra: I'm guessing you may want to try

  mat <- as.matrix(filtered_raster, wide = TRUE)

in your apply_focal_filter() function.

[jan-abel-inwt]

Thanks a lot!