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

  x <- sqrt(w$u10^2 + w$v10^2)
# 0.001s
# [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

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

2 Answers 2


A variation on method Z:

  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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.