how to annotate trajectories with stars in space and time

r-spatial / stars

Spatiotemporal Arrays, Raster and Vector Data Cubes

https://r-spatial.github.io/stars/

Apache License 2.0

554 stars 94 forks source link

how to annotate trajectories with stars in space and time #352

Closed bniebuhr closed 3 years ago

bniebuhr commented 3 years ago

Hi,

I have animal movement trajectories (with defined by (x, y, t); could also be any other trajectory) and a NetCDF file with one environmental variable (e.g. snow depth) over space and time. I want to be able to annotate my trajectory with a new column, snow_depth, in which matches this variable by the location and by the time (day) of the record. I imagine something like raster::extract that can cross both the spatial location and the "temporal location" of each animal position record. Can I do that with stars? I tried testing this with a single position with aggregate(start_object, sf_object, function(x) x[1], as_point = F), it was quite time consuming and it kept records of all times for a given position. I would also like to keep my points as a tibble or sf object, and just add a new column with the variable extracted.

I looked at that using other packages such as tidync and ncdfgeom but I could not find an easy solution for that yet, in vognettes or forums.

My data looks like:

> amt::amt_fisher
# A tibble: 14,230 x 6
         x_       y_ t_                  sex   id    name 
 *    <dbl>    <dbl> <dttm>              <chr> <chr> <chr>
 1 1782673. 2402297. 2009-02-11 12:16:45 M     M1    Leroy
 2 1782680. 2402297. 2009-02-11 12:31:38 M     M1    Leroy
 3 1782683. 2402292. 2009-02-11 12:45:48 M     M1    Leroy

and my netcdf file looks like:

> test
stars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
 sd_2007.nc [mm]   
 Min.   :  -1.000  
 1st Qu.:  -1.000  
 Median :  -1.000  
 Mean   :   6.793  
 3rd Qu.:  -1.000  
 Max.   :1464.000  
dimension(s):
     from   to         offset  delta      refsys point
x       1 1195             NA     NA SWEREF99 TM    NA
y       1 1550             NA     NA SWEREF99 TM    NA
time    1  365 2007-01-01 UTC 1 days     POSIXct    NA
                              values    
x     [1195x1550] -75000,...,1119000 [x]
y    [1195x1550] 6450000,...,7999000 [y]
time                            NULL    
curvilinear grid

I do not have a reprex but I can provide one, if necessary. It would be nice to have something like this in the package datasets to be able to test.

edzer commented 3 years ago

I agree. Would be great if you could provide a reprex, along with a script with what you did.

bniebuhr commented 3 years ago

Hi,

Maybe this could work, tell me if there is something missing for a reprex.

First I take some datasets from starsdata of ocean variables following the steps of one of the stars vignette. They vary in space and time, with daily temporal resolution.

library(stars)
library(dplyr)
library(abind)
library(lubridate)

# 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:
# 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))
z <- y %>% adrop()

# this data set look like
stars object with 3 dimensions and 4 attributes
attribute(s), summary of first 1e+05 cells:
   sst [°*C]       anom [°*C]      err [°*C]     ice [percent]  
 Min.   :-1.80   Min.   :-4.69   Min.   :0.110   Min.   :0.010  
 1st Qu.:-1.19   1st Qu.:-0.06   1st Qu.:0.300   1st Qu.:0.730  
 Median :-1.05   Median : 0.52   Median :0.300   Median :0.830  
 Mean   :-0.32   Mean   : 0.23   Mean   :0.295   Mean   :0.766  
 3rd Qu.:-0.20   3rd Qu.: 0.71   3rd Qu.:0.300   3rd Qu.:0.870  
 Max.   : 9.36   Max.   : 3.70   Max.   :0.480   Max.   :1.000  
 NA's   :13360   NA's   :13360   NA's   :13360   NA's   :27377  
dimension(s):
     from   to                   offset  delta  refsys point values    
x       1 1440                        0   0.25      NA    NA   NULL [x]
y       1  720                       90  -0.25      NA    NA   NULL [y]
time    1    9 1981-09-01 02:00:00 CEST 1 days POSIXct    NA   NULL

Then I create a random trajectory across the Atlantic Ocean with only a few positions, just to play with.

# create trajectory
# points somewhere in Atlantic Ocean
traj <- dplyr::tibble(x = runif(5, -7.814687, 10.249546),
                      y = runif(5, -34.707068, -19.200096),
                      time = lubridate::ymd("19810901") + days(c(1,2,3,5,7))) %>% 
  sf::st_as_sf(coords = c("x", "y"), remove = F)
traj

# this look like
Simple feature collection with 5 features and 3 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -6.652136 ymin: -33.8658 xmax: 8.353633 ymax: -19.87982
CRS:            NA
# A tibble: 5 x 4
       x     y time                    geometry
   <dbl> <dbl> <date>                   <POINT>
1 -4.00  -19.9 1981-09-02 (-3.997794 -19.87982)
2 -6.65  -22.0 1981-09-03 (-6.652136 -21.98427)
3  8.35  -29.1 1981-09-04  (8.353633 -29.05652)
4 -0.636 -33.9 1981-09-06 (-0.6364778 -33.8658)
5  5.52  -29.6 1981-09-08  (5.522454 -29.61568)

Let's say we want to create two new columns in this sf object, with the sst and anom variables from the z stars object, at the corresponding positions and time. What would be the best solution for that?

Sorry, it is not a good realistic example, but something I came up with since I guess I cannot share the data I am using at this moment.

edzer commented 3 years ago

Good enough! This is a first shot, e.g.

> st_extract(z, traj, time_column = "time")
Simple feature collection with 5 features and 5 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -5.29408 ymin: -32.72663 xmax: 7.563879 ymax: -20.8987
CRS:            NA
          sst        anom        err          ice       time
1    NA [°*C]    NA [°*C]   NA [°*C] NA [percent] 1981-09-02
2 17.41 [°*C]  0.32 [°*C] 0.30 [°*C] NA [percent] 1981-09-03
3 15.24 [°*C] -1.16 [°*C] 0.13 [°*C] NA [percent] 1981-09-04
4    NA [°*C]    NA [°*C]   NA [°*C] NA [percent] 1981-09-06
5    NA [°*C]    NA [°*C]   NA [°*C] NA [percent] 1981-09-08
                     geometry
1 POINT (-2.513204 -32.72663)
2  POINT (7.563879 -26.32871)
3  POINT (6.400503 -32.01026)
4 POINT (-0.4375339 -20.8987)
5  POINT (-5.29408 -31.20643)
>

time matching is still incomplete, e.g. for cases where a cube has POSIXct time stamps but indicates intervals (point is TRUE). Also the returned time is that of the data cube, not that of the pts argument.

edzer commented 3 years ago

Handles interval, point=TRUE, and returns both cube and pts time fields (as long as they don't have the same name!)

bniebuhr commented 3 years ago

Dear @edzer , thank you for that!

However, I could not test that. I tried installing stars from Github but I got the same error. I tried in different computers. And from the error, I don't have a clue on the reason...

devtools::install_github("https://github.com/r-spatial/stars.git")

Using github PAT from envvar GITHUB_PAT
Downloading GitHub repo r-spatial/stars@HEAD

The downloaded binary packages are in ...
    C:\Users\AppData\Local\Temp\RtmpaG98db\downloaded_packages
√  checking for file 'C:\Users\AppData\Local\Temp\RtmpaG98db\remotes3a643b4e7ba6\r-spatial-stars-87ad17d/DESCRIPTION' (792ms)
-  preparing 'stars': (2.3s)
√  checking DESCRIPTION meta-information ... 
-  checking for LF line-endings in source and make files and shell scripts (474ms)
-  checking for empty or unneeded directories
-  building 'stars_0.4-4.tar.gz'

* installing *source* package 'stars' ...
** using staged installation
** R
** demo
** inst
** byte-compile and prepare package for lazy loading
Error: (converted from warning) package 'sf' was built under R version 3.6.3
Execution halted
ERROR: lazy loading failed for package 'stars'
* removing 'C:/Users/Documents/R/R-3.6.1/library/stars'
* restoring previous 'C:/Users/Documents/R/R-3.6.1/library/stars'
Error: Failed to install 'stars' from GitHub:
  (converted from warning) installation of package ‘C:/Users/AppData/Local/Temp/RtmpaG98db/file3a64357f54df/stars_0.4-4.tar.gz’ had non-zero exit status

Nowosad commented 3 years ago

@bniebuhr I seems you have an old R version (R-3.6.1). Have you tried to update it to the more recent one?

bniebuhr commented 3 years ago

@Nowosad thanks for the hint, and sorry about such a dumb error. It works now.

@edzer that's fantastic, thanks! I'll try it with a real example now, and I come back here if there are issues (or even if it works).

edzer commented 3 years ago

That would be good; the implementation now does a two-step: first extract all time slices at all points, as st_extract already did, and then match the time slice to the time point. This may be both demanding in terms of resources and slow. Also, I'm not sure it can handle more than three dimensions data cubes (spatial raster + time).

bniebuhr commented 3 years ago

Ok, now I tried with a real example.

First, I noticed the extraction does not work for all projections. When I used EPSG:3006 (SWEREF99 TM), I got the error below. The description of the code is further down.

Error in colrow_from_xy(st_coordinates(pts), x, NA_outside = TRUE) : 
  colrow_from_xy not supported for curvilinear objects

# packages
library(remotes)
# remotes::install_github("r-spatial/sf") # sf version 0.97 required
# remotes::install_github("https://github.com/r-spatial/stars.git")

library(sf)
library(stars)

# parameters
crs_to_use <- 3006

# study area
# already in epsg:3006
study_area <- sf::read_sf("availability_general_mittadalen.shp")

# read snow depth layer
# originally in epsg:32633
sd_2007 <- stars::read_stars("sd_2007.nc") %>% 
  sf::st_set_crs(32633) %>%
  sf::st_transform(crs_to_use)

# cut spatiotemporal cube for the study area
# it was not necessary, but I did
sd_2007_study_area <- sd_2007[study_area]

# test for one point
pt <-  dplyr::tibble(x = 425879, y = 6891469) %>% 
  sf::st_as_sf(coords = c(1,2), crs = st_crs(study_area), remove = F) %>% 
  dplyr::mutate(day = lubridate::make_date(2007, 01, 10))

(pt_ann <- stars::st_extract(sd_2007_study_area, pt, time_column = "day", point = TRUE))

# Error in colrow_from_xy(st_coordinates(pts), x, NA_outside = TRUE) : 
#   colrow_from_xy not supported for curvilinear objects

If I then change to EPSG:32633 (WGS84, UTM zone 33N), everything works fine.

Now, when I run the real example in epsg:32633, I get this other error:


# having loaded the rest in the code above

# real example, data in epsg:3006
load("data_cleaned_winter_SK.RData")
gps_data <- gps_clean10
inds <- unique(gps_clean10$id)

# for the first individual
ind1 <- gps_clean10 %>% 
  dplyr::filter(id == inds[1]) %>% 
  sf::st_as_sf(coords = c("x", "y"), remove = F, crs = 3006) %>% 
  sf::st_transform(crs_to_use)

# extract
ind1_ann <- stars:st_extract(sd_2007_study_area, ind1, time_column = "timestamp")

#Error in stars:st_extract(sd_2007_study_area, ind1, time_column = "timestamp") : 
#  NA/NaN argument
#In addition: Warning messages:
#1: In CPL_read_gdal(as.character(x), as.character(options), as.character(driver),  :
#  GDAL Message 1: Recode from UTF-8 to CP_ACP failed with the error: "Invalid argument".
#2: In CPL_extract(f, pts, as.logical(bilinear)) :
#  GDAL Message 1: Recode from UTF-8 to CP_ACP failed with the error: "Invalid argument".
#3: In stars:st_extract(sd_2007_study_area, ind1, time_column = "timestamp") :
#  numerical expression has 4 elements: only the first used

anyNA(ind1$timestamp) # FALSE
any(is.nan(ind1$timestamp)) # FALSE

Any idea of what this can be? I cannot attach the data here but I can send it by email or share it in a private repo.

edzer commented 3 years ago

Would be great if you could post the output of traceback() after the second error.

The first one is not because of the CRS but because you used st_transform() on a raster, which creates a curvilinear grid, which (so far) has no support in st_extract. Use st_warp() instead to warp the grid to a new grid in your target CRS.

bniebuhr commented 3 years ago

Dear @edzer

Sorry, I don't know exactly what went wrong, now it is working! I put the code here again with new questions.

library(sf)
library(stars)

# parameters
crs_to_use <- 3006

# study area
study_area <- sf::read_sf("data/availability_general_mittadalen.shp") 

# read daily snow depth layer
sd_2007 <- stars::read_stars("data/sd_2007.nc") %>% 
  sf::st_set_crs(32633) %>%
  stars::st_warp(crs = crs_to_use)

# cut spatiotemporal cube for the study area
# it was not necessary, but I did it
sd_2007_study_area <- sd_2007[study_area]

plot(sd_2007_study_area[,,,1], reset = F)
plot(st_geometry(study_area), add = T)

# test for one point
pt <-  dplyr::tibble(x = 425879, y = 6891469) %>% 
  sf::st_as_sf(coords = c(1,2), crs = st_crs(study_area), remove = F) %>% 
  dplyr::mutate(day = lubridate::make_date(2007, 01, 10))

(pt_ann <- stars::st_extract(sd_2007_study_area, pt, time_column = "day", point = TRUE))

# real example
load("data/data_cleaned_winter_SK.RData")
gps_data <- gps_clean10
inds <- unique(gps_clean10$id)

# for the first individual
ind1 <- gps_clean10 %>% 
  dplyr::filter(id == inds[1]) %>% 
  dplyr::mutate(day = lubridate::date(timestamp)) %>% 
  sf::st_as_sf(coords = c("x", "y"), remove = F, crs = 3006)

# extract

# this does not work
(ind1_ann <- stars::st_extract(sd_2007_study_area, ind1, time_column = "timestamp"))

Simple feature collection with 220 features and 3 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 403276.5 ymin: 6880588 xmax: 436997.6 ymax: 6906352
projected CRS:  SWEREF99 TM
First 10 features:
   sd_2007.nc time           timestamp                 geometry
1     NA [mm] <NA> 2007-11-17 17:06:00   POINT (403312 6906345)
2     NA [mm] <NA> 2007-11-17 20:06:00 POINT (403276.5 6906352)
3     NA [mm] <NA> 2007-11-17 23:06:00 POINT (404331.1 6905805)
4     NA [mm] <NA> 2007-11-18 02:06:00 POINT (404822.7 6905665)
5     NA [mm] <NA> 2007-11-18 05:06:00   POINT (405077 6905559)
6     NA [mm] <NA> 2007-11-18 08:07:00 POINT (405142.4 6905241)
7     NA [mm] <NA> 2007-11-18 11:07:00 POINT (406259.3 6904447)
8     NA [mm] <NA> 2007-11-18 14:06:00 POINT (407219.6 6904131)
9     NA [mm] <NA> 2007-11-18 17:06:00   POINT (408226 6903713)
10    NA [mm] <NA> 2007-11-18 20:06:00 POINT (408242.8 6903706)

# this works
(ind1_ann <- stars::st_extract(sd_2007_study_area, ind1, time_column = "day"))

Simple feature collection with 220 features and 3 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 403276.5 ymin: 6880588 xmax: 436997.6 ymax: 6906352
projected CRS:  SWEREF99 TM
First 10 features:
   sd_2007.nc       time        day                 geometry
1    211 [mm] 2007-11-17 2007-11-17   POINT (403312 6906345)
2    211 [mm] 2007-11-17 2007-11-17 POINT (403276.5 6906352)
3    214 [mm] 2007-11-17 2007-11-17 POINT (404331.1 6905805)
4    182 [mm] 2007-11-18 2007-11-18 POINT (404822.7 6905665)
5    182 [mm] 2007-11-18 2007-11-18   POINT (405077 6905559)
6    182 [mm] 2007-11-18 2007-11-18 POINT (405142.4 6905241)
7    167 [mm] 2007-11-18 2007-11-18 POINT (406259.3 6904447)
8    164 [mm] 2007-11-18 2007-11-18 POINT (407219.6 6904131)
9    156 [mm] 2007-11-18 2007-11-18   POINT (408226 6903713)
10   156 [mm] 2007-11-18 2007-11-18 POINT (408242.8 6903706)

So, the function is working now, when I pick only the day of the records and use the time_column = day to match the time in the netcdf stars object. This stars object also has daily data, so this is why I guess it works.

When I use the original timestamp column (time_column = timestamp, which is a POSIXct object with a record every third hour), however, I get NA for all snow depth and stars time rows in the output. For some reason it is not finding the right day for each record. When I look at the match_time function in the source code it seems to me it should work, since it has as.Date for both objects, but for some reason it doesn't.

bniebuhr commented 3 years ago

Anyway, the way it is solves my current problem, since I can just have the day as a new column and use it.

Yet, I think expanding this would be great and there could be several options. I'll give some examples.

1) stars object with gaps: my point dataset has positions every day, and the stars object has a record everyday, but there are some gaps of a few days. If a given point is on a date for which I don't have spatial background information, it would be nice if I could choose whether to (i) just keep NA, (ii) get the closes value for which there is a value (such as the previous or next day when there is information), or even (iii) interpolating the last and next records in time, in the stars object, and getting an average value for my point.

2) stars and point objects with different temporal resolution: my points were sampled 4 times a day but in my stars object I have records only every 3 days, or at irregular times. Maybe this is the same case as number 1. If I just use st_extract as it is, I'll have 4 positions annotated with background data (for a given day), then have NA for the two following days, and so on. It would be good if we could get the closest value on time or an interpolation as options. I guess this kind of case is common when one creates spatiotemporal layers of weather variables based on interpolations of weather station data, which often has gaps.

3) objects with higher temporal resolution: my points are collected every 10 minutes and my background data (e.g. radar images) are collected every 5 minutes. Here it will be probably very hard to match exactly the same time (e.g. 2020-11-25 04:20:03) and it does not make sense to think in the interval of days, so maybe it would be necessary to find the closest timestamp to a given point time record.

Maybe you have already thought about it or it is already implemented elsewhere, but it would be great to have all these things in this same function.

bniebuhr commented 3 years ago

Another example that came to me later:

objects with different temporal resolution: my points are collected two times a day, but the stars object has a single map for each year, and the datasets encompass several years. in this case, there will not be variation across time in the variables within starts, only across space; however, there will be variation across space and time if we consider the multiple years. In this case, it would be nice to set what is the temporal resolution of the stars object and the temporal resolution of the sf (or sf_time) point object.

edzer commented 3 years ago

Could you please show the output of st_dimensions(sd_2007)?

bniebuhr commented 3 years ago

Yes.

st_dimensions(sd_2007)

     from   to         offset  delta      refsys point values x/y
x       1 1196         -75500   1000 SWEREF99 TM    NA   NULL [x]
y       1 1550        7999500  -1000 SWEREF99 TM    NA   NULL [y]
time    1  365 2007-01-01 UTC 1 days     POSIXct    NA   NULL

edzer commented 3 years ago

That is what I thought. For time, the field point is NA, indicating that (as default) the time dimension does not relate to time intervals, but to time instances. Trajectory points that do not coincide exactly with these time instances will find no match. If you'd set

attr(sd_2007, "dimensions")$time$point = FALSE

then the match will work.

The following commit will also support interpolation in time; demonstrated in the last example below:

```r library(stars) # Loading required package: abind # Loading required package: sf # Linking to GEOS 3.8.1, GDAL 3.1.3, PROJ 7.1.1 library(stars) library(dplyr) # Attaching package: ‘dplyr’ # The following objects are masked from ‘package:stats’: # filter, lag # The following objects are masked from ‘package:base’: # intersect, setdiff, setequal, union library(abind) library(lubridate) # Attaching package: ‘lubridate’ # The following objects are masked from ‘package:base’: # date, intersect, setdiff, union # 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" ) set.seed(131) # see the second vignette: # 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: # sst [°*C] anom [°*C] err [°*C] ice [percent] # Min. :-1.80 Min. :-4.69 Min. :0.110 Min. :0.010 # 1st Qu.:-1.19 1st Qu.:-0.06 1st Qu.:0.300 1st Qu.:0.730 # Median :-1.05 Median : 0.52 Median :0.300 Median :0.830 # Mean :-0.32 Mean : 0.23 Mean :0.295 Mean :0.766 # 3rd Qu.:-0.20 3rd Qu.: 0.71 3rd Qu.:0.300 3rd Qu.:0.870 # Max. : 9.36 Max. : 3.70 Max. :0.480 Max. :1.000 # NA's :13360 NA's :13360 NA's :13360 NA's :27377 # dimension(s): # from to offset delta refsys point values x/y # x 1 1440 0 0.25 NA NA NULL [x] # y 1 720 90 -0.25 NA NA NULL [y] # zlev 1 1 0 [m] NA NA NA NULL # time 1 9 1981-09-01 UTC 1 days POSIXct NA NULL z <- y %>% adrop() # create trajectory # points somewhere in Atlantic Ocean traj <- dplyr::tibble(x = runif(6, 0, 360), y = runif(6, -90, 90), time_pts = lubridate::ymd("19810901") + days(c(1,2,3,5,7,12))) %>% sf::st_as_sf(coords = c("x", "y"), remove = F) st_extract(z, traj, time_column = "time_pts") # Simple feature collection with 6 features and 6 fields # geometry type: POINT # dimension: XY # bbox: xmin: 44.97918 ymin: -47.18619 xmax: 304.6849 ymax: 74.69014 # CRS: NA # sst anom err ice time time_pts # 1 28.45 [°*C] -0.05 [°*C] 0.35 [°*C] NA [percent] 1981-09-02 1981-09-02 # 2 5.07 [°*C] 0.44 [°*C] 0.38 [°*C] NA [percent] 1981-09-03 1981-09-03 # 3 15.60 [°*C] -0.97 [°*C] 0.18 [°*C] NA [percent] 1981-09-04 1981-09-04 # 4 1.44 [°*C] -0.90 [°*C] 0.31 [°*C] NA [percent] 1981-09-06 1981-09-06 # 5 NA [°*C] NA [°*C] NA [°*C] NA [percent] 1981-09-08 1981-09-08 # 6 NA [°*C] NA [°*C] NA [°*C] NA [percent] 1981-09-13 # geometry # 1 POINT (74.31731 3.352569) # 2 POINT (44.97918 -47.18619) # 3 POINT (105.5783 -31.26274) # 4 POINT (135.2807 74.69014) # 5 POINT (304.6849 -15.07531) # 6 POINT (190.5137 -27.96455) z = z[,,,c(1,3,4,5,6,7,9)] # now, time becomes intervals st_extract(z, traj, time_column = "time_pts") # Simple feature collection with 6 features and 6 fields # geometry type: POINT # dimension: XY # bbox: xmin: 44.97918 ymin: -47.18619 xmax: 304.6849 ymax: 74.69014 # CRS: NA # sst anom err ice time time_pts # 1 NA [°*C] NA [°*C] NA [°*C] NA [percent] 1981-09-02 # 2 5.07 [°*C] 0.44 [°*C] 0.38 [°*C] NA [percent] 1981-09-03 1981-09-03 # 3 15.60 [°*C] -0.97 [°*C] 0.18 [°*C] NA [percent] 1981-09-04 1981-09-04 # 4 1.44 [°*C] -0.90 [°*C] 0.31 [°*C] NA [percent] 1981-09-06 1981-09-06 # 5 NA [°*C] NA [°*C] NA [°*C] NA [percent] 1981-09-08 # 6 NA [°*C] NA [°*C] NA [°*C] NA [percent] 1981-09-13 # geometry # 1 POINT (74.31731 3.352569) # 2 POINT (44.97918 -47.18619) # 3 POINT (105.5783 -31.26274) # 4 POINT (135.2807 74.69014) # 5 POINT (304.6849 -15.07531) # 6 POINT (190.5137 -27.96455) # verify that with time$points is TRUE, no matches are made: z <- y %>% adrop() attr(z, "dimensions")$time$point = TRUE traj$time_pts = as.POSIXct(traj$time_pts) + 1 # one second later, point=TRUE: no matches st_extract(z, traj, time_column = "time_pts") # Simple feature collection with 6 features and 6 fields # geometry type: POINT # dimension: XY # bbox: xmin: 44.97918 ymin: -47.18619 xmax: 304.6849 ymax: 74.69014 # CRS: NA # sst anom err ice time time_pts # 1 NA [°*C] NA [°*C] NA [°*C] NA [percent] 1981-09-02 02:00:01 # 2 NA [°*C] NA [°*C] NA [°*C] NA [percent] 1981-09-03 02:00:01 # 3 NA [°*C] NA [°*C] NA [°*C] NA [percent] 1981-09-04 02:00:01 # 4 NA [°*C] NA [°*C] NA [°*C] NA [percent] 1981-09-06 02:00:01 # 5 NA [°*C] NA [°*C] NA [°*C] NA [percent] 1981-09-08 02:00:01 # 6 NA [°*C] NA [°*C] NA [°*C] NA [percent] 1981-09-13 02:00:01 # geometry # 1 POINT (74.31731 3.352569) # 2 POINT (44.97918 -47.18619) # 3 POINT (105.5783 -31.26274) # 4 POINT (135.2807 74.69014) # 5 POINT (304.6849 -15.07531) # 6 POINT (190.5137 -27.96455) # verify that with time$points is TRUE, and time_col DATE, matches are made: traj$date = as.Date(traj$time_pts) st_extract(z, traj, time_column = "date") # Simple feature collection with 6 features and 6 fields # geometry type: POINT # dimension: XY # bbox: xmin: 44.97918 ymin: -47.18619 xmax: 304.6849 ymax: 74.69014 # CRS: NA # sst anom err ice time date # 1 28.45 [°*C] -0.05 [°*C] 0.35 [°*C] NA [percent] 1981-09-02 1981-09-02 # 2 5.07 [°*C] 0.44 [°*C] 0.38 [°*C] NA [percent] 1981-09-03 1981-09-03 # 3 15.60 [°*C] -0.97 [°*C] 0.18 [°*C] NA [percent] 1981-09-04 1981-09-04 # 4 1.44 [°*C] -0.90 [°*C] 0.31 [°*C] NA [percent] 1981-09-06 1981-09-06 # 5 NA [°*C] NA [°*C] NA [°*C] NA [percent] 1981-09-08 1981-09-08 # 6 NA [°*C] NA [°*C] NA [°*C] NA [percent] 1981-09-13 # geometry # 1 POINT (74.31731 3.352569) # 2 POINT (44.97918 -47.18619) # 3 POINT (105.5783 -31.26274) # 4 POINT (135.2807 74.69014) # 5 POINT (304.6849 -15.07531) # 6 POINT (190.5137 -27.96455) # verify that if the cube has Date time stamps, and pts are POSIXct or Date, matches are made: z <- y %>% adrop() z = st_set_dimensions(z, "time", as.Date(attr(z, "dimensions")$time$offset)+0:8) st_dimensions(z) # from to offset delta refsys point values x/y # x 1 1440 0 0.25 NA NA NULL [x] # y 1 720 90 -0.25 NA NA NULL [y] # time 1 9 1981-09-01 1 days Date NA NULL st_extract(z, traj, time_column = "time_pts") # Simple feature collection with 6 features and 6 fields # geometry type: POINT # dimension: XY # bbox: xmin: 44.97918 ymin: -47.18619 xmax: 304.6849 ymax: 74.69014 # CRS: NA # sst anom err ice time # 1 28.45 [°*C] -0.05 [°*C] 0.35 [°*C] NA [percent] 1981-09-02 # 2 5.07 [°*C] 0.44 [°*C] 0.38 [°*C] NA [percent] 1981-09-03 # 3 15.60 [°*C] -0.97 [°*C] 0.18 [°*C] NA [percent] 1981-09-04 # 4 1.44 [°*C] -0.90 [°*C] 0.31 [°*C] NA [percent] 1981-09-06 # 5 NA [°*C] NA [°*C] NA [°*C] NA [percent] 1981-09-08 # 6 NA [°*C] NA [°*C] NA [°*C] NA [percent] # time_pts geometry # 1 1981-09-02 02:00:01 POINT (74.31731 3.352569) # 2 1981-09-03 02:00:01 POINT (44.97918 -47.18619) # 3 1981-09-04 02:00:01 POINT (105.5783 -31.26274) # 4 1981-09-06 02:00:01 POINT (135.2807 74.69014) # 5 1981-09-08 02:00:01 POINT (304.6849 -15.07531) # 6 1981-09-13 02:00:01 POINT (190.5137 -27.96455) st_extract(z, traj, time_column = "date") # Simple feature collection with 6 features and 6 fields # geometry type: POINT # dimension: XY # bbox: xmin: 44.97918 ymin: -47.18619 xmax: 304.6849 ymax: 74.69014 # CRS: NA # sst anom err ice time date # 1 28.45 [°*C] -0.05 [°*C] 0.35 [°*C] NA [percent] 1981-09-02 1981-09-02 # 2 5.07 [°*C] 0.44 [°*C] 0.38 [°*C] NA [percent] 1981-09-03 1981-09-03 # 3 15.60 [°*C] -0.97 [°*C] 0.18 [°*C] NA [percent] 1981-09-04 1981-09-04 # 4 1.44 [°*C] -0.90 [°*C] 0.31 [°*C] NA [percent] 1981-09-06 1981-09-06 # 5 NA [°*C] NA [°*C] NA [°*C] NA [percent] 1981-09-08 1981-09-08 # 6 NA [°*C] NA [°*C] NA [°*C] NA [percent] 1981-09-13 # geometry # 1 POINT (74.31731 3.352569) # 2 POINT (44.97918 -47.18619) # 3 POINT (105.5783 -31.26274) # 4 POINT (135.2807 74.69014) # 5 POINT (304.6849 -15.07531) # 6 POINT (190.5137 -27.96455) # use POSIXct time slices, move the traj time points half a day, and interpolate time: z = st_set_dimensions(z, "time", as.POSIXct(attr(z, "dimensions")$time$offset)+(0:8)*3600*24) traj$time_pts = traj$time_pts + 12*3600 st_extract(z, traj, time_column = "time_pts", interpolate_time = TRUE) # Simple feature collection with 6 features and 6 fields # geometry type: POINT # dimension: XY # bbox: xmin: 44.97918 ymin: -47.18619 xmax: 304.6849 ymax: 74.69014 # CRS: NA # sst anom err ice time time_pts # 1 28.5650020 0.06500266 0.2599979 NA 1981-09-01 02:00:00 1981-09-02 14:00:01 # 2 5.0649998 0.42499964 0.3900002 NA 1981-09-02 02:00:00 1981-09-03 14:00:01 # 3 15.4749968 -1.09000275 0.2800023 NA 1981-09-03 02:00:00 1981-09-04 14:00:01 # 4 0.8399861 -1.50001386 0.3100000 NA 1981-09-05 02:00:00 1981-09-06 14:00:01 # 5 NA NA NA NA 1981-09-07 02:00:00 1981-09-08 14:00:01 # 6 NA NA NA NA 1981-09-13 14:00:01 # geometry # 1 POINT (74.31731 3.352569) # 2 POINT (44.97918 -47.18619) # 3 POINT (105.5783 -31.26274) # 4 POINT (135.2807 74.69014) # 5 POINT (304.6849 -15.07531) # 6 POINT (190.5137 -27.96455) ```

Your wishlist number 1 is there, try (after the example above)

z[,,,c(1,4,5)]

and look what happens with the time dimension.

bniebuhr commented 3 years ago

Wow, that's impressive! Thank you so much for that!

That is what I thought. For time, the field point is NA, indicating that (as default) the time dimension does not relate to time intervals, but to time instances. Trajectory points that do not coincide exactly with these time instances will find no match.

Ok, now I got it. Is this documented somewhere, just for reference?

I followed and reproduced the examples here, but in the end there were differences. Below I start the same way but jump to the ending example.

library(stars)
library(dplyr)
library(abind)
library(lubridate)

# 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"
)

set.seed(131)
# see the second vignette:
# 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:
#   sst [°*C]       anom [°*C]      err [°*C]     ice [percent]  
# Min.   :-1.80   Min.   :-4.69   Min.   :0.110   Min.   :0.010  
# 1st Qu.:-1.19   1st Qu.:-0.06   1st Qu.:0.300   1st Qu.:0.730  
# Median :-1.05   Median : 0.52   Median :0.300   Median :0.830  
# Mean   :-0.32   Mean   : 0.23   Mean   :0.295   Mean   :0.766  
# 3rd Qu.:-0.20   3rd Qu.: 0.71   3rd Qu.:0.300   3rd Qu.:0.870  
# Max.   : 9.36   Max.   : 3.70   Max.   :0.480   Max.   :1.000  
# NA's   :13360   NA's   :13360   NA's   :13360   NA's   :27377  
# dimension(s):
#   from   to         offset  delta  refsys point values x/y
# x       1 1440              0   0.25      NA    NA   NULL [x]
# y       1  720             90  -0.25      NA    NA   NULL [y]
# zlev    1    1          0 [m]     NA      NA    NA   NULL    
# time    1    9 1981-09-01 UTC 1 days POSIXct    NA   NULL    
z <- y %>% adrop()

# create trajectory
# points somewhere in in the World
traj <- dplyr::tibble(x = runif(6, 0, 360),
                      y = runif(6, -90, 90),
                      time_pts = lubridate::ymd("19810901") + 
                        lubridate::days(c(1,2,3,5,7,12))) %>% 
  sf::st_as_sf(coords = c("x", "y"), remove = F)
traj
# Simple feature collection with 6 features and 3 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 44.97918 ymin: -47.18619 xmax: 304.6849 ymax: 74.69014
# CRS:            NA
# # A tibble: 6 x 4
# x      y time_pts               geometry
# <dbl>  <dbl> <date>                  <POINT>
#   1  74.3   3.35 1981-09-02  (74.31731 3.352569)
# 2  45.0 -47.2  1981-09-03 (44.97918 -47.18619)
# 3 106.  -31.3  1981-09-04 (105.5783 -31.26274)
# 4 135.   74.7  1981-09-06  (135.2807 74.69014)
# 5 305.  -15.1  1981-09-08 (304.6849 -15.07531)
# 6 191.  -28.0  1981-09-13 (190.5137 -27.96455)

# for the example to work, I understand I need 
attr(z, "dimensions")$time$point <- FALSE

# use POSIXct time slices, move the traj time points half a day, and interpolate time
z = st_set_dimensions(z, "time", as.POSIXct(attr(z, "dimensions")$time$offset)+(0:8)*3600*24)
traj$time_pts = as.POSIXct(traj$time_pts) + 12*3600

# no interpolation
st_extract(z, traj, time_column = "time_pts")
# Simple feature collection with 6 features and 6 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 44.97918 ymin: -47.18619 xmax: 304.6849 ymax: 74.69014
# CRS:            NA
# sst        anom        err          ice       time            time_pts
# 1 28.45 [°*C] -0.05 [°*C] 0.35 [°*C] NA [percent] 1981-09-02 1981-09-02 14:00:00
# 2  5.07 [°*C]  0.44 [°*C] 0.38 [°*C] NA [percent] 1981-09-03 1981-09-03 14:00:00
# 3 15.60 [°*C] -0.97 [°*C] 0.18 [°*C] NA [percent] 1981-09-04 1981-09-04 14:00:00
# 4  1.44 [°*C] -0.90 [°*C] 0.31 [°*C] NA [percent] 1981-09-06 1981-09-06 14:00:00
# 5    NA [°*C]    NA [°*C]   NA [°*C] NA [percent] 1981-09-08 1981-09-08 14:00:00
# 6    NA [°*C]    NA [°*C]   NA [°*C] NA [percent]       <NA> 1981-09-13 14:00:00
# geometry
# 1  POINT (74.31731 3.352569)
# 2 POINT (44.97918 -47.18619)
# 3 POINT (105.5783 -31.26274)
# 4  POINT (135.2807 74.69014)
# 5 POINT (304.6849 -15.07531)
# 6 POINT (190.5137 -27.96455)

# bilinear interpolation
st_extract(z, traj, time_column = "time_pts", interpolate_time = T, bilinear = T)
# Simple feature collection with 6 features and 3 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 44.97918 ymin: -47.18619 xmax: 304.6849 ymax: 74.69014
# CRS:            NA
# sst       time            time_pts                   geometry
# 1 28.504345 1981-09-02 1981-09-02 14:00:00  POINT (74.31731 3.352569)
# 2  5.052837 1981-09-03 1981-09-03 14:00:00 POINT (44.97918 -47.18619)
# 3 15.878668 1981-09-04 1981-09-04 14:00:00 POINT (105.5783 -31.26274)
# 4  1.422091 1981-09-06 1981-09-06 14:00:00  POINT (135.2807 74.69014)
# 5       NaN 1981-09-08 1981-09-08 14:00:00 POINT (304.6849 -15.07531)
# 6        NA       <NA> 1981-09-13 14:00:00 POINT (190.5137 -27.96455)
# Warning message:
#   In write_stars.stars(x, fname, NA_value = NA_value) :
#   all but first attribute are ignored

As you can see, in this last example I get only values for one of the variables in the z stars object, sst, not all.

Also, I did not get exactly what this line does:

z = st_set_dimensions(z, "time", as.POSIXct(attr(z, "dimensions")$time$offset)+(0:8)*3600*24)

Wasn't this the default in the beginning?

Now I continue to understand what happens if we have gaps in time.

# see what happens if we have gaps in the stars data
(z <- z[,,,c(1,4,8)])
# stars object with 3 dimensions and 4 attributes
# attribute(s), summary of first 1e+05 cells:
#   sst [°*C]       anom [°*C]      err [°*C]     ice [percent]  
# Min.   :-1.80   Min.   :-4.69   Min.   :0.110   Min.   :0.010  
# 1st Qu.:-1.19   1st Qu.:-0.06   1st Qu.:0.300   1st Qu.:0.730  
# Median :-1.05   Median : 0.52   Median :0.300   Median :0.830  
# Mean   :-0.32   Mean   : 0.23   Mean   :0.295   Mean   :0.766  
# 3rd Qu.:-0.20   3rd Qu.: 0.71   3rd Qu.:0.300   3rd Qu.:0.870  
# Max.   : 9.36   Max.   : 3.70   Max.   :0.480   Max.   :1.000  
# NA's   :13360   NA's   :13360   NA's   :13360   NA's   :27377  
# dimension(s):
#   from   to offset delta  refsys point
# x       1 1440      0  0.25      NA    NA
# y       1  720     90 -0.25      NA    NA
# time    1    3     NA    NA POSIXct FALSE
# values x/y
# x                                                                         NULL [x]
# y                                                                         NULL [y]
# time [1981-09-01,1981-09-02), [1981-09-04,1981-09-05), [1981-09-08,1981-09-09) 

# as I understand, the default now is already the nearest neighbor, right?
st_extract(z, traj, time_column = "time_pts", interpolate_time = F, bilinear = F)
# Simple feature collection with 6 features and 6 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 44.97918 ymin: -47.18619 xmax: 304.6849 ymax: 74.69014
# CRS:            NA
# sst        anom        err          ice       time            time_pts                   geometry
# 1 28.68 [°*C]  0.18 [°*C] 0.17 [°*C] NA [percent] 1981-09-01 1981-09-02 14:00:00  POINT (74.31731 3.352569)
# 2  5.09 [°*C]  0.44 [°*C] 0.47 [°*C] NA [percent] 1981-09-01 1981-09-03 14:00:00 POINT (44.97918 -47.18619)
# 3 15.60 [°*C] -0.97 [°*C] 0.18 [°*C] NA [percent] 1981-09-04 1981-09-04 14:00:00 POINT (105.5783 -31.26274)
# 4  0.35 [°*C] -1.98 [°*C] 0.32 [°*C] NA [percent] 1981-09-04 1981-09-06 14:00:00  POINT (135.2807 74.69014)
# 5    NA [°*C]    NA [°*C]   NA [°*C] NA [percent]       <NA> 1981-09-08 14:00:00 POINT (304.6849 -15.07531)
# 6    NA [°*C]    NA [°*C]   NA [°*C] NA [percent]       <NA> 1981-09-13 14:00:00 POINT (190.5137 -27.96455)

# currently, this seem to be the same as
st_extract(z, traj, time_column = "time_pts", interpolate_time = T, bilinear = F)
# Simple feature collection with 6 features and 6 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 44.97918 ymin: -47.18619 xmax: 304.6849 ymax: 74.69014
# CRS:            NA
# sst        anom        err          ice       time            time_pts                   geometry
# 1 28.68 [°*C]  0.18 [°*C] 0.17 [°*C] NA [percent] 1981-09-01 1981-09-02 14:00:00  POINT (74.31731 3.352569)
# 2  5.09 [°*C]  0.44 [°*C] 0.47 [°*C] NA [percent] 1981-09-01 1981-09-03 14:00:00 POINT (44.97918 -47.18619)
# 3 15.60 [°*C] -0.97 [°*C] 0.18 [°*C] NA [percent] 1981-09-04 1981-09-04 14:00:00 POINT (105.5783 -31.26274)
# 4  0.35 [°*C] -1.98 [°*C] 0.32 [°*C] NA [percent] 1981-09-04 1981-09-06 14:00:00  POINT (135.2807 74.69014)
# 5    NA [°*C]    NA [°*C]   NA [°*C] NA [percent]       <NA> 1981-09-08 14:00:00 POINT (304.6849 -15.07531)
# 6    NA [°*C]    NA [°*C]   NA [°*C] NA [percent]       <NA> 1981-09-13 14:00:00 POINT (190.5137 -27.96455)

Here I have several questions.

Why doesn't the day 1981-09-08 match with the pts record 1981-09-08 14:00:00?
In both calls, why don't the first two records (that matched in time with 1981-09-01) have similar values for the variables? The same for the following two rows that match with 1981-09-04.
As I understand, the two last calls respond the same. Maybe the default st_extract(z, traj, time_column = "time_pts", interpolate_time = F, bilinear = F) could return NA for the days when there is not a time match (e.g. days 1981-09-02, 03, and 06 in the pts object above).
And the call st_extract(z, traj, time_column = "time_pts", interpolate_time = T, bilinear = F) would "interpolate" finding the nearest neighbor date and getting the value of that date.
Finally, in the nearest_neighbor case it could be nice if event values in the extreme (the last record on 2019-09-13) could take the closest value, in this case the values for 2019-09-08,

Does this all make sense?

Again, thanks a lot for all that!

edzer commented 3 years ago

Ok, now I got it. Is this documented somewhere, just for reference?

Vignette 4: stars data model does, but is not very clear what a value NA for points implies.

The interval matching was done wrongly; this should now be fixed, and I added some tooling that allows putting time intervals in a data.frame (list-)column, and format them properly. This now gives, after loading the data:

z <- z[,,,c(1,4,8)]
format(st_dimensions(z)$time$values)
# [1] "[1981-09-01 02:00:00,1981-09-02 02:00:00)"
# [2] "[1981-09-04 02:00:00,1981-09-05 02:00:00)"
# [3] "[1981-09-08 02:00:00,1981-09-09 02:00:00)"
traj$time_pts
# [1] "1981-09-02 14:00:01 CEST" "1981-09-03 14:00:01 CEST"
# [3] "1981-09-04 14:00:01 CEST" "1981-09-06 14:00:01 CEST"
# [5] "1981-09-08 14:00:01 CEST" "1981-09-13 14:00:01 CEST"
st_extract(z, traj, time_column = "time_pts")
# Simple feature collection with 6 features and 6 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 44.97918 ymin: -47.18619 xmax: 304.6849 ymax: 74.69014
# CRS:            NA
#          sst        anom        err          ice
# 1   NA [°*C]    NA [°*C]   NA [°*C] NA [percent]
# 2   NA [°*C]    NA [°*C]   NA [°*C] NA [percent]
# 3 15.6 [°*C] -0.97 [°*C] 0.18 [°*C] NA [percent]
# 4   NA [°*C]    NA [°*C]   NA [°*C] NA [percent]
# 5   NA [°*C]    NA [°*C]   NA [°*C] NA [percent]
# 6   NA [°*C]    NA [°*C]   NA [°*C] NA [percent]
#                                        time            time_pts
# 1                                   [NA,NA) 1981-09-02 14:00:01
# 2                                   [NA,NA) 1981-09-03 14:00:01
# 3 [1981-09-04 02:00:00,1981-09-05 02:00:00) 1981-09-04 14:00:01
# 4                                   [NA,NA) 1981-09-06 14:00:01
# 5 [1981-09-08 02:00:00,1981-09-09 02:00:00) 1981-09-08 14:00:01
# 6                                   [NA,NA) 1981-09-13 14:00:01
#                     geometry
# 1  POINT (74.31731 3.352569)
# 2 POINT (44.97918 -47.18619)
# 3 POINT (105.5783 -31.26274)
# 4  POINT (135.2807 74.69014)
# 5 POINT (304.6849 -15.07531)
# 6 POINT (190.5137 -27.96455)

record 5 matches with space/time, but the space POINT falls in an NA area of the data cube, hence NA attributes returned.

Finally, in the nearest_neighbor case it could be nice if event values in the extreme (the last record on 2019-09-13) could take the closest value, in this case the values for 2019-09-08,

I don't think this is a good idea, you'd go beyond the data cube extent. gdalwarp also doesn't do this AFAIR.

As you can see, in this last example I get only values for one of the variables in the z stars object, sst, not all.

this now also works, units however are dropped (because we write to disk and use GeoTIFF as format); Also NaN values are returned which should not be the case (TBD).

Wasn't this the default in the beginning?

Yes.

bniebuhr commented 3 years ago

Hi @edzer Thanks again for that!

I have tested it but I am unsure about some things.

Loading the data in the same way, and before slicing (with the complete data for the stars object z, all days), I have the result below when I use interpolate_time = TRUE. What does this interpolated value mean? In the first row, for instance, given that the timestamp is 1981-09-02 14:00:00, is sst = 28.504345 the average value of z for this pixel between the day 1981-09-02 and 1981-09-03 (since the record is in the middle of the day)? If not, what does this means?

# bilinear interpolation
st_extract(z, traj, time_column = "time_pts", interpolate_time = T, bilinear = T)
# Simple feature collection with 6 features and 6 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 44.97918 ymin: -47.18619 xmax: 304.6849 ymax: 74.69014
# CRS:            NA
# sst         anom       err ice       time            time_pts                   geometry
# 1 28.504345 -0.002270305 0.3559174 NaN 1981-09-02 1981-09-02 14:00:00  POINT (74.31731 3.352569)
# 2  5.052837  0.456668503 0.3785723 NaN 1981-09-03 1981-09-03 14:00:00 POINT (44.97918 -47.18619)
# 3 15.878668 -0.767320543 0.1800000 NaN 1981-09-04 1981-09-04 14:00:00 POINT (105.5783 -31.26274)
# 4  1.422091 -0.912794179 0.3100000 NaN 1981-09-06 1981-09-06 14:00:00  POINT (135.2807 74.69014)
# 5       NaN          NaN       NaN NaN 1981-09-08 1981-09-08 14:00:00 POINT (304.6849 -15.07531)
# 6        NA           NA        NA  NA       <NA> 1981-09-13 14:00:00 POINT (190.5137 -27.96455)

After slicing z, we have three possibilities

The first one is no interpolation at all.

(z <- z[,,,c(1,4,8)])

format(st_dimensions(z)$time$values)
#  "[1981-09-01,1981-09-02)" "[1981-09-04,1981-09-05)" "[1981-09-08,1981-09-09)"

traj$time_pts
# "1981-09-02 14:00:00 CEST" "1981-09-03 14:00:00 CEST" "1981-09-04 14:00:00 CEST" "1981-09-06 14:00:00 CEST"
# "1981-09-08 14:00:00 CEST" "1981-09-13 14:00:00 CEST"

st_extract(z, traj, time_column = "time_pts", interpolate_time = F, bilinear = F)

In this case, since we have only values for the days 1981-09-04 and 1981-09-08, the extraction retrieves only these two values (and the value for the last one is NaN because it is a position in land, and the variables are at the sea). Ok there.

The second possibility is to interpolate using bilinear interpolation.

st_extract(z, traj, time_column = "time_pts", interpolate_time = T, bilinear = T)
# Simple feature collection with 6 features and 6 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 44.97918 ymin: -47.18619 xmax: 304.6849 ymax: 74.69014
# CRS:            NA
# sst       anom  err ice                    time            time_pts                   geometry
# 1       NA         NA   NA  NA                 [NA,NA) 1981-09-02 14:00:00  POINT (74.31731 3.352569)
# 2       NA         NA   NA  NA                 [NA,NA) 1981-09-03 14:00:00 POINT (44.97918 -47.18619)
# 3 15.87867 -0.7673205 0.18 NaN [1981-09-04,1981-09-05) 1981-09-04 14:00:00 POINT (105.5783 -31.26274)
# 4       NA         NA   NA  NA                 [NA,NA) 1981-09-06 14:00:00  POINT (135.2807 74.69014)
# 5      NaN        NaN  NaN NaN [1981-09-08,1981-09-09) 1981-09-08 14:00:00 POINT (304.6849 -15.07531)
# 6       NA         NA   NA  NA                 [NA,NA) 1981-09-13 14:00:00 POINT (190.5137 -27.96455)

As it is, it interpolates the value for the 3rd row (maybe using what I understood and asked in the question 1), but no values are extracted for the rows 1, 2, and 4. And I understand the powerful use of the interpolation tool would be to fill these values based in an interpolation of the values of z, right? I expected the rows 1 and 2 would be filled with an interpolation of z[,,,1] and z[,,,4] (indices from the original unsliced stars object), and the 4th row would be filled with an interpolation between z[,,,4] and z[,,,8]. Does it make sense to work like that?

Finally, the last possibility would be a nearest neighbor approach (should I call this interpolation?). This feature was present in one of the previous versions, when the first two rows matched the record of z for 1981-09-01 and the next two rows matched 1981-09-04. Now, however, rows 1, 2, and 4 return no time match, and the result is the same whether I use interpolate_time = F, bilinear = F or interpolate_time = T, bilinear = F.

st_extract(z, traj, time_column = "time_pts", interpolate_time = T, bilinear = F)
# Simple feature collection with 6 features and 6 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 44.97918 ymin: -47.18619 xmax: 304.6849 ymax: 74.69014
# CRS:            NA
# sst  anom  err ice                    time            time_pts                   geometry
# 1   NA    NA   NA  NA                 [NA,NA) 1981-09-02 14:00:00  POINT (74.31731 3.352569)
# 2   NA    NA   NA  NA                 [NA,NA) 1981-09-03 14:00:00 POINT (44.97918 -47.18619)
# 3 15.6 -0.97 0.18  NA [1981-09-04,1981-09-05) 1981-09-04 14:00:00 POINT (105.5783 -31.26274)
# 4   NA    NA   NA  NA                 [NA,NA) 1981-09-06 14:00:00  POINT (135.2807 74.69014)
# 5   NA    NA   NA  NA [1981-09-08,1981-09-09) 1981-09-08 14:00:00 POINT (304.6849 -15.07531)
# 6   NA    NA   NA  NA                 [NA,NA) 1981-09-13 14:00:00 POINT (190.5137 -27.96455)

I hope I was clear in my points, please tell me if I should clarify them. Thanks, again!

edzer commented 3 years ago

Yes, but it's the time interpolated value between the bilinearly interpolated spatial values, not "the pixel".
The interpretation of time for z is not time intervals, so we get a match to the interval, no time interpolation. (I understand your confusion, and a warning for this case would be in place.)
No, by taking out time slices while creating z, you get time intervals for the remaining pieces, the missing time slices are now considered "missing" (NA).

bniebuhr commented 3 years ago

Hi, just an update on that.

I am annotating a trajectory with 268,609 features/points, using a stars_proxy object with daily data on snow depth, with 1km spatial resolution with 13 years of data. It worked quite well (elapsed time = 2.434597 mins) in a good computer (Windows 10 64bit, i7 3,2 GHz, 128MB RAM, R version 4.0.3). It did not work on a simpler notebook, however (Windows 10 64bit, i7 2,6 GHz, 16GB RAM, R version 4.0.3), because it reached the memory limit.

edzer commented 3 years ago

What was the object.size() of the trajectory object?

bniebuhr commented 3 years ago

> object.size(mov_track_ua_12h)
132182744 bytes

edzer commented 3 years ago

Yeah, that's pretty big, but I have no idea whether that's the cause.

jbrow247 commented 3 years ago

Hello @edzer , Thank you for working on this issue, it is something I would find useful, though I too am struggling to understand what the time interpolation is doing (or if it is actually using bilinear interpretation versus nearest-neighbor.

Following from the above example...

z <- y %>% adrop()
attr(z, "dimensions")$time$point <- FALSE ## points can't be used for interpolation, set to interval
### same location but at different times
traj <- dplyr::tibble(x = 181,
                      y = 2,
                      ts = c(lubridate::ymd_h("19810901 00") + days(c(1,2,3, 7,8)),
                             lubridate::ymd_h("19810901 00") + days(c(1,2,7,8)) + hours(12))) %>% 
  arrange(ts) %>%
  sf::st_as_sf(coords = c("x", "y"), remove = F)

Traj looks like this:

Simple feature collection with 9 features and 3 fields geometry type: POINT dimension: XY bbox: xmin: 181 ymin: 2 xmax: 181 ymax: 2 CRS: NA x y ts geometry dbl dbl dttm POINT 1 181 2 1981-09-02 00:00:00 (181 2) 2 181 2 1981-09-02 12:00:00 (181 2) 3 181 2 1981-09-03 00:00:00 (181 2) 4 181 2 1981-09-03 12:00:00 (181 2) 5 181 2 1981-09-04 00:00:00 (181 2) 6 181 2 1981-09-08 00:00:00 (181 2) 7 181 2 1981-09-08 12:00:00 (181 2) 8 181 2 1981-09-09 00:00:00 (181 2) 9 181 2 1981-09-09 12:00:00 (181 2)

I would expect that if bilinear interpolation in time is functioning, where hour(ts) = 0, this should result in uninterpolated values from z, and where hour(ts) = 12 should be the average of the ts on the day previous and following (or at least an intermediate value). The final row (9), should return NAs, because it is off the grid, where as row 8 should work, because it is on the boarder of the grid.

After applying st_extract with interpolate_time = T: st_extract(z, traj, bilinear = T, time_column = "ts", interpolate_time = T) Output:

Simple feature collection with 9 features and 6 fields geometry type: POINT dimension: XY bbox: xmin: 181 ymin: 2 xmax: 181 ymax: 2 CRS: NA sst anom err ice time ts geometry 1 28.8725 0.1050 0.1200 NaN 1981-09-02 1981-09-02 00:00:00 POINT (181 2) 2 28.8725 0.1050 0.1200 NaN 1981-09-02 1981-09-02 12:00:00 POINT (181 2) 3 29.1825 0.4150 0.1200 NaN 1981-09-03 1981-09-03 00:00:00 POINT (181 2) 4 29.1825 0.4150 0.1200 NaN 1981-09-03 1981-09-03 12:00:00 POINT (181 2) 5 28.9100 0.1375 0.1275 NaN 1981-09-04 1981-09-04 00:00:00 POINT (181 2) 6 29.4000 0.6200 0.3050 NaN 1981-09-08 1981-09-08 00:00:00 POINT (181 2) 7 29.4000 0.6200 0.3050 NaN 1981-09-08 1981-09-08 12:00:00 POINT (181 2) 8 NA NA NA NA 1981-09-09 00:00:00 POINT (181 2) 9 NA NA NA NA 1981-09-09 12:00:00 POINT (181 2)

What I instead get is values obtained from the intermediate time stamps (rows 2, 4, 7), matching proceeding day in z (which matches the returned 'time' column), rather than interpolating between the proceeding and following. Further, row 8, which should have a match (though on the border) in z, is returning NA.

I'm not sure if this is my misunderstanding of how temporal interpolation should work, or an issue with the function.

Thanks for your help! (and please let me know if there is better way to report the data & code above!)

edzer commented 3 years ago

@jbrow247 please tidy your issue so that it is clear where the R code is that I should copy and paste.

jbrow247 commented 3 years ago

@jbrow247 please tidy your issue so that it is clear where the R code is that I should copy and paste.

@edzer Sorry, I've edited the formatting above. Hopefully this is better.

also - my stars package version is 0.5-1 and R version 4.0.2

edzer commented 3 years ago

Thanks! This now gives

library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.0
library(dplyr)

# Attaching package: ‘dplyr’

# The following objects are masked from ‘package:stats’:

#     filter, lag

# The following objects are masked from ‘package:base’:

#     intersect, setdiff, setequal, union

library(abind)
library(lubridate)

# Attaching package: ‘lubridate’

# The following objects are masked from ‘package:base’:

#     date, intersect, setdiff, union

# 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:
# 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:
#    sst [°*C]       anom [°*C]      err [°*C]     ice [percent]  
#  Min.   :-1.80   Min.   :-4.69   Min.   :0.110   Min.   :0.010  
#  1st Qu.:-1.19   1st Qu.:-0.06   1st Qu.:0.300   1st Qu.:0.730  
#  Median :-1.05   Median : 0.52   Median :0.300   Median :0.830  
#  Mean   :-0.32   Mean   : 0.23   Mean   :0.295   Mean   :0.766  
#  3rd Qu.:-0.20   3rd Qu.: 0.71   3rd Qu.:0.300   3rd Qu.:0.870  
#  Max.   : 9.36   Max.   : 3.70   Max.   :0.480   Max.   :1.000  
#  NA's   :13360   NA's   :13360   NA's   :13360   NA's   :27377  
# dimension(s):
#      from   to         offset  delta  refsys point values x/y
# x       1 1440              0   0.25      NA    NA   NULL [x]
# y       1  720             90  -0.25      NA    NA   NULL [y]
# zlev    1    1          0 [m]     NA      NA    NA   NULL    
# time    1    9 1981-09-01 UTC 1 days POSIXct    NA   NULL    

z <- y %>% adrop()
attr(z, "dimensions")$time$point <- FALSE ## points can't be used for interpolation, set to interval
z
# stars object with 3 dimensions and 4 attributes
# attribute(s), summary of first 1e+05 cells:
#    sst [°*C]       anom [°*C]      err [°*C]     ice [percent]  
#  Min.   :-1.80   Min.   :-4.69   Min.   :0.110   Min.   :0.010  
#  1st Qu.:-1.19   1st Qu.:-0.06   1st Qu.:0.300   1st Qu.:0.730  
#  Median :-1.05   Median : 0.52   Median :0.300   Median :0.830  
#  Mean   :-0.32   Mean   : 0.23   Mean   :0.295   Mean   :0.766  
#  3rd Qu.:-0.20   3rd Qu.: 0.71   3rd Qu.:0.300   3rd Qu.:0.870  
#  Max.   : 9.36   Max.   : 3.70   Max.   :0.480   Max.   :1.000  
#  NA's   :13360   NA's   :13360   NA's   :13360   NA's   :27377  
# dimension(s):
#      from   to         offset  delta  refsys point values x/y
# x       1 1440              0   0.25      NA    NA   NULL [x]
# y       1  720             90  -0.25      NA    NA   NULL [y]
# time    1    9 1981-09-01 UTC 1 days POSIXct FALSE   NULL    

### same location but at different times
traj <- dplyr::tibble(x = 181,
                      y = 2,
                      ts = c(lubridate::ymd_h("19810901 00") + days(c(1,2,3, 7,8)),
                             lubridate::ymd_h("19810901 00") + days(c(1,2,7,8)) + hours(12))) %>%
  arrange(ts) %>%
  sf::st_as_sf(coords = c("x", "y"), remove = F)
st_extract(z, traj, bilinear = T, time_column = "ts", interpolate_time = T)
# Simple feature collection with 9 features and 6 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 181 ymin: 2 xmax: 181 ymax: 2
# CRS:            NA
#       sst   anom    err ice       time                  ts      geometry
# 1 28.8725 0.1050 0.1200 NaN 1981-09-02 1981-09-02 00:00:00 POINT (181 2)
# 2 28.8725 0.1050 0.1200 NaN 1981-09-02 1981-09-02 12:00:00 POINT (181 2)
# 3 29.1825 0.4150 0.1200 NaN 1981-09-03 1981-09-03 00:00:00 POINT (181 2)
# 4 29.1825 0.4150 0.1200 NaN 1981-09-03 1981-09-03 12:00:00 POINT (181 2)
# 5 28.9100 0.1375 0.1275 NaN 1981-09-04 1981-09-04 00:00:00 POINT (181 2)
# 6 29.4000 0.6200 0.3050 NaN 1981-09-08 1981-09-08 00:00:00 POINT (181 2)
# 7 29.4000 0.6200 0.3050 NaN 1981-09-08 1981-09-08 12:00:00 POINT (181 2)
# 8 29.4000 0.6200 0.3050 NaN 1981-09-08 1981-09-09 00:00:00 POINT (181 2)
# 9      NA     NA     NA  NA       <NA> 1981-09-09 12:00:00 POINT (181 2)

attr(z, "dimensions")$time$point <- TRUE ## time points: are used for interpolation
st_extract(z, traj, bilinear = T, time_column = "ts", interpolate_time = T)
# Simple feature collection with 9 features and 6 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 181 ymin: 2 xmax: 181 ymax: 2
# CRS:            NA
#        sst    anom     err ice       time                  ts      geometry
# 1 28.87250 0.10500 0.12000 NaN 1981-09-02 1981-09-02 00:00:00 POINT (181 2)
# 2 29.02750 0.26000 0.12000 NaN 1981-09-02 1981-09-02 12:00:00 POINT (181 2)
# 3 29.18250 0.41500 0.12000 NaN 1981-09-03 1981-09-03 00:00:00 POINT (181 2)
# 4 29.04625 0.27625 0.12375 NaN 1981-09-03 1981-09-03 12:00:00 POINT (181 2)
# 5 28.91000 0.13750 0.12750 NaN 1981-09-04 1981-09-04 00:00:00 POINT (181 2)
# 6 29.40000 0.62000 0.30500 NaN 1981-09-08 1981-09-08 00:00:00 POINT (181 2)
# 7 29.38250 0.60375 0.43250 NaN 1981-09-08 1981-09-08 12:00:00 POINT (181 2)
# 8 29.36500 0.58750 0.56000 NaN 1981-09-09 1981-09-09 00:00:00 POINT (181 2)
# 9       NA      NA      NA  NA       <NA> 1981-09-09 12:00:00 POINT (181 2)

jbrow247 commented 3 years ago

@edzer Thanks for the speedy fix! I'm looking forward to trying it on my own data!

edzer commented 3 years ago

Can we close this?

bniebuhr commented 3 years ago

I think so, at least for the time being! Thanks, @edzer

henningte commented 2 years ago

This is now also possible using the sftime package which extents sf objects to handle irregular time information.

Here is the last example, now using an sftime object together with st_extract():

library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1; sf_use_s2() is TRUE
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(abind)
library(lubridate)
#> 
#> Attaching package: 'lubridate'
#> The following objects are masked from 'package:base':
#> 
#>     date, intersect, setdiff, union
library(sftime)

# 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:
# 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 point values x/y
#> x       1 1440              0   0.25      NA    NA   NULL [x]
#> y       1  720             90  -0.25      NA    NA   NULL [y]
#> zlev    1    1          0 [m]     NA      NA    NA   NULL    
#> time    1    9 1981-09-01 UTC 1 days POSIXct    NA   NULL

z <- y %>% adrop()
attr(z, "dimensions")$time$point <- FALSE ## points can't be used for interpolation, set to interval
z
#> stars object with 3 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 point values x/y
#> x       1 1440              0   0.25      NA    NA   NULL [x]
#> y       1  720             90  -0.25      NA    NA   NULL [y]
#> time    1    9 1981-09-01 UTC 1 days POSIXct FALSE   NULL

### same location but at different times (now with sftime object)
traj <- dplyr::tibble(x = 181,
                      y = 2,
                      ts = c(lubridate::ymd_h("19810901 00") + days(c(1,2,3, 7,8)),
                             lubridate::ymd_h("19810901 00") + days(c(1,2,7,8)) + hours(12))) %>%
  arrange(ts) %>%
  sf::st_as_sf(coords = c("x", "y"), remove = F) %>%
  sftime::st_as_sftime(time_column_name = "ts")

st_extract(z, traj, bilinear = T, interpolate_time = T)
#> Simple feature collection with 9 features and 6 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 181 ymin: 2 xmax: 181 ymax: 2
#> CRS:           NA
#>       sst   anom    err ice       time                  ts      geometry
#> 1 28.8725 0.1050 0.1200 NaN 1981-09-02 1981-09-02 00:00:00 POINT (181 2)
#> 2 28.8725 0.1050 0.1200 NaN 1981-09-02 1981-09-02 12:00:00 POINT (181 2)
#> 3 29.1825 0.4150 0.1200 NaN 1981-09-03 1981-09-03 00:00:00 POINT (181 2)
#> 4 29.1825 0.4150 0.1200 NaN 1981-09-03 1981-09-03 12:00:00 POINT (181 2)
#> 5 28.9100 0.1375 0.1275 NaN 1981-09-04 1981-09-04 00:00:00 POINT (181 2)
#> 6 29.4000 0.6200 0.3050 NaN 1981-09-08 1981-09-08 00:00:00 POINT (181 2)
#> 7 29.4000 0.6200 0.3050 NaN 1981-09-08 1981-09-08 12:00:00 POINT (181 2)
#> 8 29.4000 0.6200 0.3050 NaN 1981-09-08 1981-09-09 00:00:00 POINT (181 2)
#> 9      NA     NA     NA  NA       <NA> 1981-09-09 12:00:00 POINT (181 2)

attr(z, "dimensions")$time$point <- TRUE ## time points: are used for interpolation
st_extract(z, traj, bilinear = T, time_column = "ts", interpolate_time = T)
#> Simple feature collection with 9 features and 6 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 181 ymin: 2 xmax: 181 ymax: 2
#> CRS:           NA
#>        sst    anom     err ice       time                  ts      geometry
#> 1 28.87250 0.10500 0.12000 NaN 1981-09-02 1981-09-02 00:00:00 POINT (181 2)
#> 2 29.02750 0.26000 0.12000 NaN 1981-09-02 1981-09-02 12:00:00 POINT (181 2)
#> 3 29.18250 0.41500 0.12000 NaN 1981-09-03 1981-09-03 00:00:00 POINT (181 2)
#> 4 29.04625 0.27625 0.12375 NaN 1981-09-03 1981-09-03 12:00:00 POINT (181 2)
#> 5 28.91000 0.13750 0.12750 NaN 1981-09-04 1981-09-04 00:00:00 POINT (181 2)
#> 6 29.40000 0.62000 0.30500 NaN 1981-09-08 1981-09-08 00:00:00 POINT (181 2)
#> 7 29.38250 0.60375 0.43250 NaN 1981-09-08 1981-09-08 12:00:00 POINT (181 2)
#> 8 29.36500 0.58750 0.56000 NaN 1981-09-09 1981-09-09 00:00:00 POINT (181 2)
#> 9       NA      NA      NA  NA       <NA> 1981-09-09 12:00:00 POINT (181 2)

^{Created on 2022-04-29 by the reprex package (v2.0.1)}

Session info

``` r sessioninfo::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 4.2.0 (2022-04-22 ucrt) #> os Windows 10 x64 (build 19044) #> system x86_64, mingw32 #> ui RTerm #> language (EN) #> collate German_Germany.utf8 #> ctype German_Germany.utf8 #> tz Europe/Berlin #> date 2022-04-29 #> pandoc 2.17.1.1 @ C:/Program Files/RStudio/bin/quarto/bin/ (via rmarkdown) #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date (UTC) lib source #> abind * 1.4-5 2016-07-21 [1] CRAN (R 4.2.0) #> class 7.3-20 2022-01-16 [2] CRAN (R 4.2.0) #> classInt 0.4-3 2020-04-07 [1] CRAN (R 4.2.0) #> cli 3.3.0 2022-04-25 [1] CRAN (R 4.2.0) #> crayon 1.5.1 2022-03-26 [1] CRAN (R 4.2.0) #> DBI 1.1.2 2021-12-20 [1] CRAN (R 4.2.0) #> digest 0.6.29 2021-12-01 [1] CRAN (R 4.2.0) #> dplyr * 1.0.8 2022-02-08 [1] CRAN (R 4.2.0) #> e1071 1.7-9 2021-09-16 [1] CRAN (R 4.2.0) #> ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.2.0) #> evaluate 0.15 2022-02-18 [1] CRAN (R 4.2.0) #> fansi 1.0.3 2022-03-24 [1] CRAN (R 4.2.0) #> fastmap 1.1.0 2021-01-25 [1] CRAN (R 4.2.0) #> fs 1.5.2 2021-12-08 [1] CRAN (R 4.2.0) #> generics 0.1.2 2022-01-31 [1] CRAN (R 4.2.0) #> glue 1.6.2 2022-02-24 [1] CRAN (R 4.2.0) #> highr 0.9 2021-04-16 [1] CRAN (R 4.2.0) #> htmltools 0.5.2 2021-08-25 [1] CRAN (R 4.2.0) #> KernSmooth 2.23-20 2021-05-03 [2] CRAN (R 4.2.0) #> knitr 1.38 2022-03-25 [1] CRAN (R 4.2.0) #> lifecycle 1.0.1 2021-09-24 [1] CRAN (R 4.2.0) #> lubridate * 1.8.0 2021-10-07 [1] CRAN (R 4.2.0) #> lwgeom 0.2-8 2021-10-06 [1] CRAN (R 4.2.0) #> magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.2.0) #> pillar 1.7.0 2022-02-01 [1] CRAN (R 4.2.0) #> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.2.0) #> proxy 0.4-26 2021-06-07 [1] CRAN (R 4.2.0) #> purrr 0.3.4 2020-04-17 [1] CRAN (R 4.2.0) #> R.cache 0.15.0 2021-04-30 [1] CRAN (R 4.2.0) #> R.methodsS3 1.8.1 2020-08-26 [1] CRAN (R 4.2.0) #> R.oo 1.24.0 2020-08-26 [1] CRAN (R 4.2.0) #> R.utils 2.11.0 2021-09-26 [1] CRAN (R 4.2.0) #> R6 2.5.1 2021-08-19 [1] CRAN (R 4.2.0) #> Rcpp 1.0.8.3 2022-03-17 [1] CRAN (R 4.2.0) #> reprex 2.0.1 2021-08-05 [1] CRAN (R 4.2.0) #> rlang 1.0.2 2022-03-04 [1] CRAN (R 4.2.0) #> rmarkdown 2.13 2022-03-10 [1] CRAN (R 4.2.0) #> rstudioapi 0.13 2020-11-12 [1] CRAN (R 4.2.0) #> sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.2.0) #> sf * 1.0-8 2022-04-29 [1] Github (r-spatial/sf@1281482) #> sftime * 0.2.0.9000 2022-04-29 [1] local #> stars * 0.5-5 2021-12-19 [1] CRAN (R 4.2.0) #> stringi 1.7.6 2021-11-29 [1] CRAN (R 4.2.0) #> stringr 1.4.0 2019-02-10 [1] CRAN (R 4.2.0) #> styler 1.7.0 2022-03-13 [1] CRAN (R 4.2.0) #> tibble 3.1.6 2021-11-07 [1] CRAN (R 4.2.0) #> tidyselect 1.1.2 2022-02-21 [1] CRAN (R 4.2.0) #> units 0.8-0 2022-02-05 [1] CRAN (R 4.2.0) #> utf8 1.2.2 2021-07-24 [1] CRAN (R 4.2.0) #> vctrs 0.4.1 2022-04-13 [1] CRAN (R 4.2.0) #> withr 2.5.0 2022-03-03 [1] CRAN (R 4.2.0) #> xfun 0.30 2022-03-02 [1] CRAN (R 4.2.0) #> yaml 2.3.5 2022-02-21 [1] CRAN (R 4.2.0) #> #> [1] C:/Users/henni/AppData/Local/R/win-library/4.2 #> [2] C:/Program Files/R/R-4.2.0/library #> #> ────────────────────────────────────────────────────────────────────────────── ```