Stars: Extracting cell values by point gives 'NA' for duplicated points

[michaeldorman]

Extracting cell values from stars object by sf points seems to produce NA for all duplicated (same x-y) points.

To extract cell values, I'm using the method mentioned in the wiki:

aggregate(stars_object, sf_object, function(x) x[1], as_points = FALSE)

Here is a reproducible example.

First, creating a sample raster and points, with some of the points being duplicated:

# Create 'stars' raster
library(stars)
d = st_dimensions(x = 1:ncol(volcano), y = 1:nrow(volcano))
r = st_as_stars(t(volcano))
r = st_set_dimensions(r, 1, offset = 0, delta = 1)
r = st_set_dimensions(r, 2, offset = nrow(volcano), delta = -1)

# Create 'sf' points
pnt = data.frame(
  x = c(30, 30, 48, 48, 48, 12)+0.5, 
  y = c(23, 23, 41, 41, 41, 71)+0.5
)
pnt = st_as_sf(pnt, coords = c("x", "y"))

Here is a plot of the stars and sf layers:

plot(r, reset = FALSE)
plot(st_geometry(pnt), add = TRUE, col = "red")

And this is the extraction result, which unexpectedly gives NA for the duplicated points:

# Extract - 'stars'
x = aggregate(r, pnt, function(x) x[1], as_points = FALSE)
x[[1]]

## [1] 135  NA 123  NA  NA 130

For comparison, extracting with raster::extract, or with sf::st_join, works as expected, i.e., returning duplicated cell values rather than NA:

# Extract - 'raster'
library(raster)
extract(as(r, "Raster"), pnt)

## [1] 135 135 123 123 123 130

# Extract - 'sf'
st_join(pnt, st_as_sf(r))$A1

## [1] 135 135 123 123 123 130

Thanks!

[dvictori]

Just bumped into the same problem here. Was about to report it when I saw your issue. I was not aware that we could use st_join. So will do that for now. Thanks

[edzer]

I need to rethink whether this is good, or bad. I guess the thinking now behind it is that by aggregating, you assign elements of x to maximally one feature in by, so that if you'd sum things, that you don't get double counts.

Using aggregate to extract at points was a hack anyway!

[michaeldorman]

I understand, thanks!

[edzer]

I'm trying to develop a robust st_join.stars method, in branch st_join, that would do a

> st_join(r, pnt)$A1
[1] 135 135 123 123 123 130

so without the need for st_as_sf(r) (which is expensive); it works for 2D, but may still change quite a bit; will bring it up as a new issue.