Helper function to get raster resolution

[kadyb]

As far as I know, the only way to find the raster resolution is to use the st_dimensions(). In {raster} and {terra} there is dedicated res() function. Could you consider adding dedicated st_res() function to {stars} too?

library("stars")
tif = system.file("tif/L7_ETMs.tif", package = "stars")
r = read_stars(tif)

st_dimensions(r)

st_res = function(x) {
  res = sapply(st_dimensions(x), "[[", "delta")
  res[!is.na(res)]
}

st_res(r)
#>     x     y
#>  28.5 -28.5

[edzer]

Yes, great idea. Maybe it is easier to let the NA in; the current implementation also reports the time resolution.

[kadyb]

Here the second attempt:

# spatial - if TRUE, it returns only the spatial resolution of the raster,
# otherwise it returns the resolution of all dimensions (can be NA); logical
st_res = function(x, spatial = TRUE) {
  res = sapply(st_dimensions(x), "[[", "delta")
  res[1:2] = abs(res[1:2])
  if (spatial) {
    res[1:2]
  } else res
}

tif = system.file("nc/reduced.nc", package = "stars")
r = read_stars(tif)
st_res(r, spatial = TRUE)
#> x y
#> 2 2
st_res(r, spatial = FALSE)
#>   x    y zlev time
#>   2    2   NA   NA

I'm not sure if 1:2 indexing always refers to XY coordinates, so maybe it would be better to use select by name?

Edit: I also added abs() to the results, but I'm not sure. {terra} returns positive values, but QGIS positive and negative.

[edzer]

Will spatial also be the vertical dimension? I like to have the option for only x/y raster; getting the x/y raster dimensions is done by

> attr(st_dimensions(r), "raster")$dimensions
[1] "x" "y"

If we want to do resolutions of all dimensions, we may have to return a list rather than a numeric vector, and provide a class & print method, as resolutions will have units.

[kadyb]

Will spatial also be the vertical dimension?

I used "spatial resolution" in terms of remote sensing, but from a GIS perspective, I think the vertical dimension should be included as "spatial" too. I rather misnamed this argument, because my intention was to return only XY resolution.

I like to have the option for only x/y raster

+1

[kadyb]

Thanks!