通过时间和空间匹配从栅格中提取点

Question

我有一个 strast，其中包含我根据 u 和 v 分量计算出的当前速度信息。

如何根据时间和空间匹配从该 SPATRAST 中提取点到数据帧？

用于重新创建示例数据的代码：

library(stars)
library(lubridate)
library(terra)

# create a datetime vector holding 364 entries at one hour intervals####
datetime <- seq(as.POSIXct("2020-01-01 00:00:00"), as.POSIXct("2020-01-17 00:00:00"), by="hour")
datetime <- datetime[1:364]

# Create a sample rast ####
v <- rast(ncol=37, nrow=74, nlyr=364, extent=ext(-6.3875, -5.4875, 52.2375, 54.0625), 
          crs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0", time=datetime)
v <- init(v, runif)
 
# Create a sample df 
time <- as.POSIXct(c("2020-01-01 15:21:00", "2020-01-14 08:01:00", "2020-01-16 21:19:00"))
lat <- c("53.92265", "52.86578", "53.38290")
long <-c("-6.026689", "-5.819075", "-6.028343")
sighting <- c("1",  "0", "1")
df <- data.frame(sighting, lat, long, time)
DF<- st_as_sf(df, coords = c("long", "lat")) 
DF = st_set_crs(DF, "EPSG:4326") 
#DF = st_transform(DF, 2157) # CRS ITM

我厌倦了使用来自星星的 st_extract:

# save spatrast as a stars object
v.star <-st_as_stars(v)
#v.star = st_transform(v.star, 2157) # CRS ITM

# extract current velocities to DF
DF$cur_vel <- st_extract(v.star, DF, time_column = "time")

但是，这会从重复的 DF 中返回带有时间和几何形状的 DF，以及带有 NA 的 cur_vel 列：

函数好像需要完全匹配？

为了避免这个问题，我使用 data.table 将 DF 的时间与栅格中最近的时间进行匹配

# Drop geom of DF and make a data.table
DF <- st_transform(DF, 4326)
DF <- cbind(as.data.table(st_drop_geometry(DF)),
                      st_coordinates(DF))

# save all the times from the spatrast in a data.table
time <- time(v)
DT = data.table(
  ID = 1:364,
  time = time
)

# Extract the nearest rast times (DT/v) to the times in the DF
DF.rast.time <- DF[, c("RASTtime") :=
                     DT[DF, on = c("time"), roll = 'nearest', .(x.time)]][]

这会产生以下数据框

然后我从栅格中提取当前速度值，但使用 DF 中的新 RASTtimes：

# save DF.rast.time as an sf object again
DF.rast.time<- st_as_sf(DF.rast.time, coords = c("X", "Y")) 
DF.rast.time = st_set_crs(DF.rast.time, "EPSG:4326") 

# extract current velocities to DF
DF.rast.time$cur_vel <- st_extract(v.star, DF.rast.time, time_column = "RASTtime") %>%
  st_as_sf()

这产生了我正在寻找的结果：

但是，是否可以阻止 RASTtime 和几何的复制列出现在这个新数据框中？

此外，是否有更直接的方法来查找数据帧和栅格之间时间上最接近的匹配？

Answer 1

示例数据

library(terra)
datetime <- as.POSIXct("2020-01-01 00:00:00") + (0:363) * 3600
r <- rast(ncol=37, nrow=74, nlyr=364, extent=ext(-6.3875, -5.4875, 52.2375, 54.0625), time=datetime)
set.seed(0)
r <- init(r, runif)
 
time <- as.POSIXct(c("2020-01-01 15:21:00", "2020-01-14 08:01:00", "2020-01-15 21:19:00"))
lat <- c(53.92265, 52.86578, 53.38290)
long <-c(-6.026689, -5.819075, -6.028343)
sighting <- c("1",  "0", "1")
df <- data.frame(sighting, lat, long, time)

SpatRaster 的时间以小时为单位。因此，我将 data.frame 中的时间四舍五入到最接近的小时，然后使用

match

找到感兴趣的图层，并在

extract

中使用它。

m <-  match(round(df$time, units="hours"), time(r))
e <- extract(r, df[, c("long", "lat")], layer=m)
e
#  ID   layer     value
#1  1  lyr.16 0.4852105
#2  2 lyr.321 0.5526752
#3  3 lyr.358 0.8640279

或者像这样，分两步：

ee <- extract(r, df[, c("long", "lat")], ID=FALSE)
ee[cbind(1:length(m), m)]
# [1] 0.4852105 0.5526752 0.8640279

我认为

match

忽略了时区的差异。如果您这样做，则会考虑时区

m <- sapply(df$time, \(i) which.min(abs(time(r)-i)))

这很好，因为它避免了舍入并且更通用。但是您必须确保 SpatRaster 和 data.frame 的设置正确。

例如

datetime <- as.POSIXct("2020-01-01 00:00:00", tz="America/Los_Angeles") + (0:363) * 3600
datetime[1]
# [1] "2020-01-01 PST"

请参阅

Sys.timezone()

了解您的时区

最后，如果您需要将经/纬度数据转换为 SpatRaster 的 crs，您可以这样做

v <- vect(df, c("long", "lat"), crs="+proj=longlat")
pv <- project(v, r)
e <- extract(r, pv, layer=m)

Answer 2

现在使用

stars

包

st_extract

功能变得更加容易，之前答案的作者已更新为以这种方式工作。借用之前答案中的示例数据...

library(terra)
library(stars)

# example data
datetime <- as.POSIXct("2020-01-01 00:00:00") + (0:363) * 3600
r <- rast(ncol = 37, nrow = 74, nlyr = 364, 
          extent = ext(-6.3875, -5.4875, 52.2375, 54.0625), time = datetime, 
          crs = "epsg:4326")
set.seed(0)
r <- init(r, runif)

time <- as.POSIXct(c("2020-01-01 15:21:00", "2020-01-14 08:01:00", "2020-01-15 21:19:00"))
lat <- c(53.92265, 52.86578, 53.38290)
long <-c(-6.026689, -5.819075, -6.028343)
sighting <- c("1",  "0", "1")
df <- data.frame(sighting, lat, long, time)

# cast raster to stars object
str <- st_as_stars(r)
str
#stars object with 3 dimensions and 1 attribute
#3attribute(s):
#              Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
#lyr.1  9.95351e-07 0.2498533 0.5004527 0.4999927 0.7496667 0.9999999
#dimension(s):
#     from  to                  offset    delta  refsys x/y
#x       1  37                  -6.388  0.02432  WGS 84 [x]
#y       1  74                   54.06 -0.02466  WGS 84 [y]
#time    1 364 2020-01-01 05:00:00 UTC  1 hours POSIXct 

# create sf object from lat, long, time
st_sf <- st_as_sf(df, coords = c("long", "lat"), crs = 4326)
st_sf
#Simple feature collection with 3 features and 2 fields
#Geometry type: POINT
#Dimension:     XY
#Bounding box:  xmin: -6.028343 ymin: 52.86578 xmax: -5.819075 ymax: 53.92265
#Geodetic CRS:  WGS 84
#  sighting                time                   geometry
#1        1 2020-01-01 15:21:00 POINT (-6.026689 53.92265)
#2        0 2020-01-14 08:01:00 POINT (-5.819075 52.86578)
#3        1 2020-01-15 21:19:00  POINT (-6.028343 53.3829)

# extract with tim
st_ext <- st_extract(str, st_sf, time_column = "time", interpolate_time = TRUE)
st_ext
#Simple feature collection with 3 features and 2 fields
#Geometry type: POINT
#Dimension:     XY
#Bounding box:  xmin: -6.028343 ymin: 52.86578 xmax: -5.819075 ymax: 53.92265
#Geodetic CRS:  WGS 84
#      lyr.1                time                   geometry
#1 0.4726441 2020-01-01 15:21:00 POINT (-6.026689 53.92265)
#2 0.5549729 2020-01-14 08:01:00 POINT (-5.819075 52.86578)
#3 0.7544850 2020-01-15 21:19:00  POINT (-6.028343 53.3829)