# What is the canonical way to perform raster calculations using R stars?

I have a data cube that I am handling in the R `stars` package. It has 3 dimensions (latitude, longitude and time) and 2 attributes (northwards and eastwards components of wind). A dummy equivalent dataset can be generated as follows:

``````library(units); library(stars)

(w <- system.file("nc/bcsd_obs_1999.nc", package = "stars") %>%
read_stars() %>% drop_units() %>% setNames(., c("u10", "v10")))
``````

I want to calculate the resultant magnitude of wind for every pixel using simple trigonometry. My data is considerably larger (87m pixels), so efficiency of computation is very important

What is the canonical way to perform such raster calculations?

Ideally I'm looking for answers that are generalisable to any raster calculation problem using R `stars`, not just the specific example I give here - mostly because I know I'll need to be doing more complicated calculations in the future.

I've found three options:

Method X

``````system.time(
x <- sqrt(w\$u10^2 + w\$v10^2)
)
# 0.001s
class(x)
# [1] array
x <- st_as_stars(x, dimensions = attr(w, "dimensions")) # Hack it back to a stars object
``````

This is the most intuitive method (at least to me) and is very fast - but the resulting object has lost the dimensions and is no longer a `stars` object, so presumably this isn't the canonical method.

Method Y

``````system.time(
y <- sqrt(w[1,,,]^2 + w[2,,,]^2)
)
# 0.006s
``````

This generates a `stars` object but takes 6 times as long as the previous method and is less intuitive.

Method Z

``````wm <- merge(w)
system.time(
z <- st_apply(wm, 1:3, function(ws) sqrt(sum(ws^2)) )
)
# 0.167s
``````

This seems the most canonically "`stars`" method, but requires reorganisation of the data and fitting a function that could get complicated and prone to errors if computing something involving many attribute layers. Plus it's super slow even discounting the required merge operation to convert the attributes to a dimension.

A variation on method Z:

``````system.time(
z <- st_apply(wm, 1:3, function(u,v) sqrt(u^2+v^2) )
)
``````

The speed-up comes from the fact that in method Z, the function is called once for each pixel, in this variation it is called only once with the sub-matrices as argument. See also here.

Revisiting this code a month later I thought of an obvious alternative method that is both intuitive and retains the `stars` object:

``````x <- sqrt(w["u10"]^2 + w["v10"]^2)
``````

This is still slightly slower than method X, even when including converting the resulting object back to a `stars` object (twice as long in the latter case). But it's faster and more elegant than any of the other methods.