使用stars::st_extract()遇单位转换错误，如何提取时间轴岩心样本？

背景 #

我使用stars包加载NetCDS对象创建了一个3D数据立方体，该立方体包含(x; y)索引的地理2D网格，以及代表时间的第三维度。我的目标是从存储在二维矩阵中的一组(x; y)坐标提取时间维度的岩心样本（core samples）。stars::st_extract()函数及其at参数似乎是合适的工具，官方文档对该参数的说明如下：

at
为sf或sfc类几何对象、每行存储坐标点的两列矩阵（用于指定提取x值的位置），或是包含几何和时间维度的stars对象（矢量数据立方体）

问题 #

当我在带时间维度的对象上使用at参数传入二维坐标矩阵时，出现如下错误：

24 error(s) in InterpolateAtPoint()
24 error(s) in InterpolateAtPoint()
Error:
! cannot convert C into mm/m

使用stars包自带的bcsd_obs数据集可复现该问题：

library(stars)
# 存储x ; y坐标的两列矩阵
pnt_mat = st_coordinates(st_sample(st_as_sfc(st_bbox(bcsd_obs)), 10))

st_extract(bcsd_obs, at = pnt_mat)
#> 24 error(s) in InterpolateAtPoint()
#> 24 error(s) in InterpolateAtPoint()
#> Error:
#> ! cannot convert C into mm/m

解决方案 #

这个错误源于st_extract处理带单位的多维数据时，单纯的坐标矩阵缺少坐标系元数据，引发单位转换冲突。只需将坐标矩阵转换为sf/sfc几何对象，即可正确提取每个坐标点对应的时间序列岩心样本：

方法1：直接生成sfc对象 #

library(stars)
library(sf)

# 生成坐标点并转为sfc格式
pnt_sfc = st_as_sfc(st_sample(st_as_sfc(st_bbox(bcsd_obs)), 10))

# 提取时间维度岩心样本
result = st_extract(bcsd_obs, at = pnt_sfc)

# 查看结果：每个点对应一条完整的时间序列
print(result)

方法2：将现有矩阵转为sf对象 #

如果已经有现成的坐标矩阵，可将其转换为sf点对象后传入：

library(stars)
library(sf)

# 存储x ; y坐标的两列矩阵
pnt_mat = st_coordinates(st_sample(st_as_sfc(st_bbox(bcsd_obs)), 10))

# 转为sf对象，指定与原数据集一致的坐标系
pnt_sf = st_as_sf(data.frame(pnt_mat), coords = c("X", "Y"), crs = st_crs(bcsd_obs))

# 提取样本
result = st_extract(bcsd_obs, at = pnt_sf)

原理说明 #

sf/sfc对象会携带完整的坐标系信息，能与原stars数据集的坐标系精准匹配，避免单位转换时的冲突；而单纯的坐标矩阵仅包含数值，缺少元数据支撑，导致函数无法正确解析维度单位，进而报错。

内容的提问来源于stack exchange，提问作者Paul