st_extract() with time matching is not generalizable to stars objects with n-ary dimensions (n>4) #739
Description
Extracting values from a stars object with x, y, date and other dimensions, gives incomplete/incorrect results.
Finding the temporal indices for time matching should work differently, considering the remainder dimensions and possibly identifying which one is the temporal dimension.
library(sf)
library(stars)
library(sftime)
date = seq(
from = as.Date("2020-10-01"),
to = as.Date("2020-10-03"),
length.out = 2
)
set.seed(1285)
x = seq(-0.5, 1, length.out = 20)
y = seq(0, 1.5, length.out = 20)
s1 = st_as_stars(data.frame(expand.grid(x = x, y = y), z = 1:400/400))
s2 = st_as_stars(data.frame(expand.grid(x = x, y = y), z = 400:1/400))
rdc = c(s1,s2, along = "date") |>
st_set_dimensions(
"date",
values = date,
) |>
st_set_crs(4326)
(rmulti1 = c(rdc, rdc*100, along = "foo"))
#> stars object with 4 dimensions and 1 attribute
#> attribute(s):
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> z 0.0025 0.499375 0.99375 25.31312 50.0625 100
#> dimension(s):
#> from to offset delta refsys x/y
#> x 1 20 -0.5395 0.07895 WGS 84 [x]
#> y 1 20 1.539 -0.07895 WGS 84 [y]
#> date 1 2 2020-10-01 2 days Date
#> foo 1 2 NA NA NA
(rmulti2 = aperm(rmulti1, c("x", "y", "foo", "date")))
#> stars object with 4 dimensions and 1 attribute
#> attribute(s):
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> z 0.0025 0.499375 0.99375 25.31312 50.0625 100
#> dimension(s):
#> from to offset delta refsys x/y
#> x 1 20 -0.5395 0.07895 WGS 84 [x]
#> y 1 20 1.539 -0.07895 WGS 84 [y]
#> foo 1 2 NA NA NA
#> date 1 2 2020-10-01 2 days Date
set.seed(82192)
pnt = st_sample(st_as_sfc(st_bbox(rdc)), 2)
(pntsf = st_sf(date = date, geometry = pnt) |> st_as_sftime())
#> Spatiotemporal feature collection with 2 features and 0 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 0.5551259 ymin: 1.321098 xmax: 0.6256081 ymax: 1.404827
#> Geodetic CRS: WGS 84
#> Time column with class: 'Date'.
#> Ranging from 2020-10-01 to 2020-10-03.
#> geometry date
#> 1 POINT (0.5551259 1.404827) 2020-10-01
#> 2 POINT (0.6256081 1.321098) 2020-10-03
## works fine
st_extract(rdc, pnt) |> st_as_sf(long = TRUE)
#> Simple feature collection with 4 features and 2 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 0.5551259 ymin: 1.321098 xmax: 0.6256081 ymax: 1.404827
#> Geodetic CRS: WGS 84
#> date z geometry
#> 1 2020-10-01 0.9350 POINT (0.5551259 1.404827)
#> 2 2020-10-01 0.8875 POINT (0.6256081 1.321098)
#> 3 2020-10-03 0.0675 POINT (0.5551259 1.404827)
#> 4 2020-10-03 0.1150 POINT (0.6256081 1.321098)
st_extract(rdc, pntsf)
#> Simple feature collection with 2 features and 2 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 0.5551259 ymin: 1.321098 xmax: 0.6256081 ymax: 1.404827
#> Geodetic CRS: WGS 84
#> z date geometry
#> 1 0.935 2020-10-01 POINT (0.5551259 1.404827)
#> 2 0.115 2020-10-03 POINT (0.6256081 1.321098)
## multidimensional - date as third dimension
st_extract(rmulti1, pnt) |> st_as_sf(long = TRUE)
#> Simple feature collection with 8 features and 3 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 0.5551259 ymin: 1.321098 xmax: 0.6256081 ymax: 1.404827
#> Geodetic CRS: WGS 84
#> date foo z geometry
#> 1 2020-10-01 1 0.9350 POINT (0.5551259 1.404827) *
#> 2 2020-10-01 1 0.8875 POINT (0.6256081 1.321098)
#> 3 2020-10-03 1 0.0675 POINT (0.5551259 1.404827)
#> 4 2020-10-03 1 0.1150 POINT (0.6256081 1.321098) *
#> 5 2020-10-01 2 93.5000 POINT (0.5551259 1.404827) **
#> 6 2020-10-01 2 88.7500 POINT (0.6256081 1.321098)
#> 7 2020-10-03 2 6.7500 POINT (0.5551259 1.404827)
#> 8 2020-10-03 2 11.5000 POINT (0.6256081 1.321098) **
# Incomplete result (missing lines marked with ** above)
st_extract(rmulti1, pntsf)
#> Simple feature collection with 2 features and 2 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 0.5551259 ymin: 1.321098 xmax: 0.6256081 ymax: 1.404827
#> Geodetic CRS: WGS 84
#> z date geometry
#> 1 0.935 2020-10-01 POINT (0.5551259 1.404827)
#> 2 0.115 2020-10-03 POINT (0.6256081 1.321098)
## multidimensional - date as fourth dimension
st_extract(rmulti2, pnt) |> st_as_sf(long = TRUE)
#> Simple feature collection with 8 features and 3 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 0.5551259 ymin: 1.321098 xmax: 0.6256081 ymax: 1.404827
#> Geodetic CRS: WGS 84
#> foo date z geometry
#> 1 1 2020-10-01 0.9350 POINT (0.5551259 1.404827) *
#> 2 1 2020-10-01 0.8875 POINT (0.6256081 1.321098)
#> 3 2 2020-10-01 93.5000 POINT (0.5551259 1.404827) **
#> 4 2 2020-10-01 88.7500 POINT (0.6256081 1.321098) ?
#> 5 1 2020-10-03 0.0675 POINT (0.5551259 1.404827)
#> 6 1 2020-10-03 0.1150 POINT (0.6256081 1.321098) **
#> 7 2 2020-10-03 6.7500 POINT (0.5551259 1.404827)
#> 8 2 2020-10-03 11.5000 POINT (0.6256081 1.321098) **
# Incorrect result (adds a non-matching value (?) and skips others (**))
st_extract(rmulti2, pntsf)
#> Simple feature collection with 2 features and 2 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 0.5551259 ymin: 1.321098 xmax: 0.6256081 ymax: 1.404827
#> Geodetic CRS: WGS 84
#> z date geometry
#> 1 0.935 2020-10-01 POINT (0.5551259 1.404827)
#> 2 88.750 2020-10-03 POINT (0.6256081 1.321098)Quick fix
Restrict time matching to only spatiotemporal dimensions.
sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#> setting value
#> version R version 4.4.2 (2024-10-31)
#> os Ubuntu 24.04.1 LTS
#> system x86_64, linux-gnu
#> ui X11
#> language (EN)
#> collate en_US.UTF-8
#> ctype en_US.UTF-8
#> tz Etc/UTC
#> date 2025-02-15
#> pandoc 3.5 @ /usr/bin/ (via rmarkdown)
#> quarto 1.5.57 @ /usr/local/bin/quarto
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#> package * version date (UTC) lib source
#> abind * 1.4-8 2024-09-12 [1] RSPM (R 4.4.0)
#> class 7.3-22 2023-05-03 [2] CRAN (R 4.4.2)
#> classInt 0.4-10 2023-09-05 [1] RSPM (R 4.4.0)
#> cli 3.6.3 2024-06-21 [1] RSPM (R 4.4.0)
#> DBI 1.2.3 2024-06-02 [1] RSPM (R 4.4.0)
#> digest 0.6.37 2024-08-19 [1] RSPM (R 4.4.0)
#> e1071 1.7-16 2024-09-16 [1] RSPM (R 4.4.0)
#> evaluate 1.0.1 2024-10-10 [1] RSPM (R 4.4.0)
#> fastmap 1.2.0 2024-05-15 [1] RSPM (R 4.4.0)
#> fs 1.6.5 2024-10-30 [1] RSPM (R 4.4.0)
#> glue 1.8.0 2024-09-30 [1] RSPM (R 4.4.0)
#> htmltools 0.5.8.1 2024-04-04 [1] RSPM (R 4.4.0)
#> KernSmooth 2.23-24 2024-05-17 [2] CRAN (R 4.4.2)
#> knitr 1.48 2024-07-07 [1] RSPM (R 4.4.0)
#> lifecycle 1.0.4 2023-11-07 [1] RSPM (R 4.4.0)
#> lwgeom 0.2-15 2024-11-03 [1] Github (r-spatial/lwgeom@4569f09)
#> magrittr 2.0.3 2022-03-30 [1] RSPM (R 4.4.0)
#> proxy 0.4-27 2022-06-09 [1] RSPM (R 4.4.0)
#> Rcpp 1.0.13 2024-07-17 [1] RSPM (R 4.4.0)
#> reprex 2.1.1 2024-07-06 [1] RSPM (R 4.4.0)
#> rlang 1.1.4 2024-06-04 [1] RSPM (R 4.4.0)
#> rmarkdown 2.28 2024-08-17 [1] RSPM (R 4.4.0)
#> rstudioapi 0.17.1 2024-10-22 [1] RSPM (R 4.4.0)
#> s2 1.1.7 2024-07-17 [1] RSPM (R 4.4.0)
#> sessioninfo 1.2.3 2025-02-05 [1] RSPM (R 4.4.0)
#> sf * 1.0-19 2024-11-03 [1] Github (r-spatial/sf@dc02f71)
#> sftime * 0.3.0 2024-09-11 [1] RSPM (R 4.4.0)
#> stars * 0.6-9 2025-02-15 [1] Github (loreabad6/stars@main)
#> units 0.8-5 2023-11-28 [1] RSPM (R 4.4.0)
#> withr 3.0.2 2024-10-28 [1] RSPM (R 4.4.0)
#> wk 0.9.4 2024-10-11 [1] RSPM (R 4.4.0)
#> xfun 0.48 2024-10-03 [1] RSPM (R 4.4.0)
#> yaml 2.3.10 2024-07-26 [1] RSPM (R 4.4.0)
#>
#> [1] /usr/local/lib/R/site-library
#> [2] /usr/local/lib/R/library
#> * ── Packages attached to the search path.
#>
#> ──────────────────────────────────────────────────────────────────────────────