Detect temporal dimension based on "refsys" in st_extract()

[loreabad6]

Hi @edzer

I created this PR to fix the following line:

stars/R/extract.R

Line 92 in 06e04a1

tm = match("time", names(st_dimensions(x))) # FIXME: select based on refsys in time classes

In the PR, the temporal dimension is automatically detected, no matter its name, based on refsys ("POSIXct", "POSIXt", "Date", "PCICt"). In case there is more than one, the first one will be taken. If there is no temporal dimension, there is an error message. A reprex of the changes:

library(stars)
library(dplyr)
# spatiotemporal data from the oceans,
# adapted from a stars vignette
x = c(
  "avhrr-only-v2.19810901.nc",
  "avhrr-only-v2.19810902.nc",
  "avhrr-only-v2.19810903.nc",
  "avhrr-only-v2.19810904.nc",
  "avhrr-only-v2.19810905.nc",
  "avhrr-only-v2.19810906.nc",
  "avhrr-only-v2.19810907.nc",
  "avhrr-only-v2.19810908.nc",
  "avhrr-only-v2.19810909.nc"
)
# see the second vignette:
# options(timeout = max(1000, getOption("timeout")))
# install.packages("starsdata", repos = "http://pebesma.staff.ifgi.de", type = "source")
file_list = system.file(paste0("netcdf/", x), package = "starsdata")
(y = read_stars(file_list, quiet = TRUE))
#> stars object with 4 dimensions and 4 attributes
#> attribute(s), summary of first 1e+05 cells:
#>                Min. 1st Qu. Median       Mean 3rd Qu. Max.  NA's
#> sst [°*C]     -1.80   -1.19  -1.05 -0.3201670   -0.20 9.36 13360
#> anom [°*C]    -4.69   -0.06   0.52  0.2299385    0.71 3.70 13360
#> err [°*C]      0.11    0.30   0.30  0.2949421    0.30 0.48 13360
#> ice [percent]  0.01    0.73   0.83  0.7657695    0.87 1.00 27377
#> dimension(s):
#>      from   to         offset  delta  refsys x/y
#> x       1 1440              0   0.25      NA [x]
#> y       1  720             90  -0.25      NA [y]
#> zlev    1    1          0 [m]     NA      NA    
#> time    1    9 1981-09-01 UTC 1 days POSIXct

# Create sf object with time column
t = y |> 
  st_bbox() |> 
  st_sample(9) |> 
  st_as_sf() |> 
  mutate(ts = st_get_dimension_values(y, 'time'))
z = y |> adrop()
st_extract(z, t, time_column = "ts")
#> Simple feature collection with 9 features and 6 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 46.93849 ymin: -46.62712 xmax: 345.7044 ymax: 87.64539
#> CRS:           NA
#>           sst       anom        err            ice       time         ts
#> 1  2.73 [°*C] 0.84 [°*C] 0.30 [°*C]   NA [percent] 1981-09-01 1981-09-01
#> 2  6.65 [°*C] 0.97 [°*C] 0.31 [°*C] 0.18 [percent] 1981-09-02 1981-09-02
#> 3 15.52 [°*C] 0.21 [°*C] 0.23 [°*C]   NA [percent] 1981-09-03 1981-09-03
#> 4    NA [°*C]   NA [°*C]   NA [°*C]   NA [percent] 1981-09-04 1981-09-04
#> 5  5.13 [°*C] 0.76 [°*C] 0.20 [°*C]   NA [percent] 1981-09-05 1981-09-05
#> 6    NA [°*C]   NA [°*C]   NA [°*C]   NA [percent] 1981-09-06 1981-09-06
#> 7 -1.18 [°*C] 0.57 [°*C] 0.30 [°*C] 0.84 [percent] 1981-09-07 1981-09-07
#> 8    NA [°*C]   NA [°*C]   NA [°*C]   NA [percent] 1981-09-08 1981-09-08
#> 9 23.92 [°*C] 0.93 [°*C] 0.11 [°*C]   NA [percent] 1981-09-09 1981-09-09
#>                            x
#> 1  POINT (150.5463 73.36827)
#> 2  POINT (59.94191 68.90588)
#> 3  POINT (321.4299 49.91673)
#> 4  POINT (263.9572 33.92841)
#> 5 POINT (345.7044 -46.62712)
#> 6  POINT (46.93849 20.81668)
#> 7  POINT (126.5627 87.64539)
#> 8  POINT (299.7781 1.075033)
#> 9  POINT (214.7992 32.88862)

# Different temporal dimension name
z1 = z |> st_set_dimensions(names = c("x", "y", "datetime"))
st_extract(z1, t, time_column = "ts")
#> Simple feature collection with 9 features and 6 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 46.93849 ymin: -46.62712 xmax: 345.7044 ymax: 87.64539
#> CRS:           NA
#>           sst       anom        err            ice   datetime         ts
#> 1  2.73 [°*C] 0.84 [°*C] 0.30 [°*C]   NA [percent] 1981-09-01 1981-09-01
#> 2  6.65 [°*C] 0.97 [°*C] 0.31 [°*C] 0.18 [percent] 1981-09-02 1981-09-02
#> 3 15.52 [°*C] 0.21 [°*C] 0.23 [°*C]   NA [percent] 1981-09-03 1981-09-03
#> 4    NA [°*C]   NA [°*C]   NA [°*C]   NA [percent] 1981-09-04 1981-09-04
#> 5  5.13 [°*C] 0.76 [°*C] 0.20 [°*C]   NA [percent] 1981-09-05 1981-09-05
#> 6    NA [°*C]   NA [°*C]   NA [°*C]   NA [percent] 1981-09-06 1981-09-06
#> 7 -1.18 [°*C] 0.57 [°*C] 0.30 [°*C] 0.84 [percent] 1981-09-07 1981-09-07
#> 8    NA [°*C]   NA [°*C]   NA [°*C]   NA [percent] 1981-09-08 1981-09-08
#> 9 23.92 [°*C] 0.93 [°*C] 0.11 [°*C]   NA [percent] 1981-09-09 1981-09-09
#>                            x
#> 1  POINT (150.5463 73.36827)
#> 2  POINT (59.94191 68.90588)
#> 3  POINT (321.4299 49.91673)
#> 4  POINT (263.9572 33.92841)
#> 5 POINT (345.7044 -46.62712)
#> 6  POINT (46.93849 20.81668)
#> 7  POINT (126.5627 87.64539)
#> 8  POINT (299.7781 1.075033)
#> 9  POINT (214.7992 32.88862)

# No temporal dimension
z2 = z1[,,,1, drop = T]
st_extract(z2, t, time_column = "ts")
#> Error in st_extract.stars(z2, t, time_column = "ts"): cannot match times: x does not have a temporal dimension

# More than one temporal dimension
z3 = st_set_dimensions(y, 3, values = st_get_dimension_values(z1, "datetime")[1])
st_extract(z3, t, time_column = "ts")
#> Simple feature collection with 9 features and 6 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 46.93849 ymin: -46.62712 xmax: 345.7044 ymax: 87.64539
#> CRS:           NA
#>          sst       anom       err          ice       zlev         ts
#> 1 2.73 [°*C] 0.84 [°*C] 0.3 [°*C] NA [percent] 1981-09-01 1981-09-01
#> 2   NA [°*C]   NA [°*C]  NA [°*C] NA [percent]       <NA> 1981-09-02
#> 3   NA [°*C]   NA [°*C]  NA [°*C] NA [percent]       <NA> 1981-09-03
#> 4   NA [°*C]   NA [°*C]  NA [°*C] NA [percent]       <NA> 1981-09-04
#> 5   NA [°*C]   NA [°*C]  NA [°*C] NA [percent]       <NA> 1981-09-05
#> 6   NA [°*C]   NA [°*C]  NA [°*C] NA [percent]       <NA> 1981-09-06
#> 7   NA [°*C]   NA [°*C]  NA [°*C] NA [percent]       <NA> 1981-09-07
#> 8   NA [°*C]   NA [°*C]  NA [°*C] NA [percent]       <NA> 1981-09-08
#> 9   NA [°*C]   NA [°*C]  NA [°*C] NA [percent]       <NA> 1981-09-09
#>                            x
#> 1  POINT (150.5463 73.36827)
#> 2  POINT (59.94191 68.90588)
#> 3  POINT (321.4299 49.91673)
#> 4  POINT (263.9572 33.92841)
#> 5 POINT (345.7044 -46.62712)
#> 6  POINT (46.93849 20.81668)
#> 7  POINT (126.5627 87.64539)
#> 8  POINT (299.7781 1.075033)
#> 9  POINT (214.7992 32.88862)

[edzer]

Great, thanks!

[edzer]

The defaults for time_column btw stem from using sftime, as in

library(sftime)
tt = st_as_sftime(t)
st_extract(z3, tt)