将spTimer包的NYdata转换为spacetime包STFDF格式求助

问题背景 #

我有一个来自spTimer包的数据集NYdata，结构如下：

> str(NYdata)
'data.frame':   1736 obs. of  10 variables:
 $ s.index  : int  1 1 1 1 1 1 1 1 1 1 ...
 $ Longitude: num  -73.8 -73.8 -73.8 -73.8 -73.8 ...
 $ Latitude : num  42.7 42.7 42.7 42.7 42.7 ...
 $ Year     : int  2006 2006 2006 2006 2006 2006 2006 2006 2006 2006 ...
 $ Month    : int  7 7 7 7 7 7 7 7 7 7 ...
 $ Day      : int  1 2 3 4 5 6 7 8 9 10 ...
 $ o8hrmax  : num  53.9 57.1 72 36.6 42.6 ...
 $ cMAXTMP  : num  27.9 30.1 30 27.9 25.7 ...
 $ WDSP     : num  5.46 8.21 4.46 3.69 4.37 ...
 $ RH       : num  2.77 3.2 3.23 4.36 3.95 ...

需要将其转换为spacetime包中的STFDF格式，目标格式参考spcopula包的EU_RB数据集：

> str(EU_RB)
Formal class 'STFDF' [package "spacetime"] with 4 slots
  ..@ data   :'data.frame': 11834 obs. of  2 variables:
  .. ..$ PM10: num [1:11834] 14 9.7 7.8 21.9 11.2 9 11 6.1 7.4 7.4 ...
  .. ..$ EMEP: num [1:11834] 11.62 5.02 3.94 3.82 7.01 ...
  ..@ sp     :Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  .. .. ..@ data       :'data.frame':   194 obs. of  1 variable:
  .. .. .. ..$ station_altitude: int [1:194] 525 581 918 560 172 117 665 1137 330 330 ...
  .. .. ..@ coords.nrs : num(0) 
  .. .. ..@ coords     : num [1:194, 1:2] 4592866 4761515 4658756 4690954 4799839 ...
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : NULL
  .. .. .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
  .. .. ..@ bbox       : num [1:2, 1:2] 2749697 1647732 6412269 4604814
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
  .. .. .. .. ..$ : chr [1:2] "min" "max"
  .. .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. .. .. ..@ projargs: chr " +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"
  ..@ time   :An xts object on 2005-06-01 / 2005-07-31 containing: 
  Data:    integer [61, 1]
  Columns: ..1
  Index:   POSIXct,POSIXt [61] (TZ: "GMT")
  ..@ endTime: POSIXct[1:61], format: "2005-06-02 03:00:00" "2005-06-03 03:00:00" "2005-06-04 03:00:00" ...

之前尝试的代码未成功：

library(sf)
library(xts)
library(st)

# Create a sf data frame with spatial coordinates
sf_NYdata <- st_as_sf(NYdata, coords = c("Longitude", "Latitude"), crs = 4326)

# Create a POSIXct vector for time index
time_index <- as.POSIXct(paste0(NYdata$Year, "-", NYdata$Month, "-", NYdata$Day))

# Find the minimum time value and calculate the offset
time_offset <- abs(min(time_index)) + 1  # Add 1 to ensure non-negative values

# Add the offset to the time index
time_index <- time_index + time_offset

# Create an xts object for time series data
xts_NYdata <- xts(NYdata[, -c(1:3, 4:6)], order.by = time_index)

# Combine spatial and time components into an STFDF object
STFDF_NYdata <- STFDF(sf_NYdata, xts_NYdata)

# View the resulting STFDF object
head(STFDF_NYdata)

问题分析 #

原代码的核心问题：

• 使用sf对象而非spacetime原生支持的sp类对象
• 时间处理错误，无需添加偏移，且未提取唯一时间点
• 空间部分未去重，每个时间点重复了相同的空间坐标

正确转换代码 #

# 加载依赖包
library(spacetime)
library(sp)
library(dplyr)
library(xts)

# 1. 处理空间部分：提取唯一空间点（按s.index分组，同一s.index对应同一坐标）
unique_stations <- NYdata %>%
  distinct(s.index, .keep_all = TRUE) %>%
  select(s.index, Longitude, Latitude)

# 转换为SpatialPointsDataFrame
sp_points <- SpatialPointsDataFrame(
  coords = unique_stations[, c("Longitude", "Latitude")],
  data = unique_stations[, "s.index", drop = FALSE],
  proj4string = CRS("+init=epsg:4326")  # 指定WGS84坐标系
)

# 2. 处理时间部分：生成唯一时间序列并转为xts对象
NYdata$date <- as.POSIXct(paste(NYdata$Year, NYdata$Month, NYdata$Day, sep = "-"), tz = "GMT")
unique_dates <- unique(NYdata$date)
time_xts <- xts(seq_along(unique_dates), order.by = unique_dates)

# 3. 整理观测数据：保留空间索引、时间索引和变量列
obs_data <- NYdata %>%
  select(s.index, date, o8hrmax, cMAXTMP, WDSP, RH)

# 4. 创建STFDF对象
ny_stfdf <- STFDF(
  sp = sp_points,
  time = time_xts,
  data = obs_data,
  spatialIndex = obs_data$s.index,  # 匹配空间点的索引
  temporalIndex = match(obs_data$date, unique_dates)  # 匹配时间点的索引
)

# 查看转换结果
str(ny_stfdf)

关键说明 #

• 空间部分：必须是无重复的空间点集合，用s.index去重即可，确保每个空间点只出现一次
• 时间部分：需要是唯一的时间序列，用unique()提取后转为xts对象，spacetime要求时间为xts或POSIXct格式
• 数据匹配：通过spatialIndex和temporalIndex参数，将观测数据与对应的空间点、时间点一一对应，保证数据的时空一致性
• 坐标系：指定正确的投影坐标系（这里用WGS84，EPSG:4326），确保空间数据的有效性

内容的提问来源于stack exchange，提问作者Dr. Statistics