0

first post :)

I've been transitioning my R code from sp() to sf()/stars(), and one thing I'm still trying to grasp is accounting for the area in my grids.

Here's an example code to explain what I mean.

library(stars)
library(tidyverse)

# Reading in an example tif file, from stars() vignette
tif = system.file("tif/L7_ETMs.tif", package = "stars")
x = read_stars(tif)
x

# Get areas for each grid of the x object. Returns stars object with "area" in units of [m^2]
x_area <- st_area(x)
x_area

I tried loosely adopting code from this vignette (https://github.com/r-spatial/stars/blob/master/vignettes/stars5.Rmd) to divide each value in x by it's grid area, and it's not working as expected (perhaps because my objects are stars and not sf?)

x$test1 = x$L7_ETMs.tif / x_area  # Some computationally intensive calculation seems to happen, but doesn't produce the results I expect?

x$test1 = x$L7_ETMs.tif / x_area$area # Throws error, "non-conformable arrays"

What does seem to work is the following.

x %>%
  mutate(test1 = L7_ETMs.tif / units::set_units(as.numeric(x_area$area), m^2))

Here are the concerns I have with this code.

  1. I worry that as I turn the x_area$area (a matrix, areas in lat/lon) into a numeric vector, I may mess up the lat/lon matching between the grid and it's area. I did some rough testing to see if the areas match up the way I expect them to, but can't escape the worry that this could lead to errors that are difficult to catch.

  2. It just doesn't seem clean that I start with "x_area" in the correct units, only to remove then set the units again during the computation.

Can someone suggest a "cleaner" implementation for what I'm trying to do, i.e. multiplying or dividing grids by its area while maintaining units throughout? Or convince me that the code I have is fine?

Thanks!

JJK
  • 3
  • 1

1 Answers1

0

I do not know how to improve the stars code, but you can compare the results you get with this

tif <- system.file("tif/L7_ETMs.tif", package = "stars")
library(terra)
r <- rast(tif)

a <- cellSize(r, sum=FALSE)
x <- r / a

With planar data you could do this when it is safe to assume there is no distortion (generally not the case, but it can be the case)

y <- r / prod(res(r))
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • Thanks for your comments! The "planar cells" aspect of the code in my post was an unintended consequence of my (lazy) attempt to create a reproducible example, my real data is lon/lats with variations in the grid areas. Also, thanks for pointing out the "terra" package (yours?), I'm reading through the documentation and find it very interesting (including perhaps the very best, most straight-forward explanation of how rasters are different from spatial vectors that I had ever read!), will definitely explore what I might be able to do with it! – JJK Feb 24 '21 at 19:29
  • I found `terra::cellSize()` to be the function I was looking for to cross-check my `st_area()` calculations. Also note that `st_area()` requires the sf/stars object to have a CRS to return areas in m^2 units. The "non-conforming array" error in my original post was actually due to a matrix dimension mismatch problem, as the areas just have 2 dimensions (x, y), while `x` in my example above had a third dimension (`band`). – JJK Jun 16 '21 at 23:46