背景

我有一个通过使用 stars 包加载的NetCDS对象所创建的三维数据立方体。这个数据立方体表示一个地理上的二维网格，具有 (x ; y) 索引，第三个维度代表时间。

我的目标是对存储在二维矩阵中的一组定义好的 (x ; y) 坐标，提取时间维度的“核心样本”。函数 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)
# Two-column matrix of x ; y coordinates
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

提问

如何使用坐标矩阵提取时间轴的核心样本？

解决方案

自 stars_0.7-3 起，这个错误不再发生。相反，将返回一个警告：

Warning message:
In st_extract.stars(bcsd_obs, at = pnt_mat) :
incompatible units merged in result matrix: dropping all units

更多细节请参见此 stars commit 。为避免由坐标与NA单元格重合而触发的（无害的）“InterpolateAtPoint()”中n 个错误，请参阅下面的原始答案。

在更新尚未提交到CRAN之前，您需要通过GitHub安装：

remotes::install_github("r-spatial/stars")

packageVersion("stars")
# [1] ‘0.7.3’

library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE

set.seed(42)
pnt_mat <- st_coordinates(st_sample(st_as_sfc(st_bbox(bcsd_obs)), 10))

st_extract(bcsd_obs, at = pnt_mat)
# 36 error(s) in InterpolateAtPoint()
# 36 error(s) in InterpolateAtPoint()
#         foo1  foo2   foo3   foo4  foo5   foo6   foo7   foo8   foo9  foo10 foo11 foo12      bar1     bar2      bar3     bar4     bar5     bar6     bar7     bar8     bar9    bar10     bar11    bar12 
#  [1,]   0.00  0.00   0.00   0.00  0.00   0.00   0.00   0.00   0.00   0.00  0.00  0.00  0.000000 0.000000  0.000000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.000000 0.000000
#  [2,]   0.00  0.00   0.00   0.00  0.00   0.00   0.00   0.00   0.00   0.00  0.00  0.00  0.000000 0.000000  0.000000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.000000 0.000000
#  [3,] 137.46 77.88  73.14 104.22 75.40 111.56 116.36  77.40  42.71  39.91 87.68 41.21  2.075806 2.844286  3.915323 12.69133 15.73952 20.29300 22.74774 21.28564 16.91600 11.64452  8.644333 2.129194
#  [4,]   0.00  0.00   0.00   0.00  0.00   0.00   0.00   0.00   0.00   0.00  0.00  0.00  0.000000 0.000000  0.000000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.000000 0.000000
#  [5,] 193.52 52.43  65.58 106.83 67.35 104.01 114.63 110.89 459.17 168.18 35.45 40.20  8.712420 8.152500  9.455483 17.23667 19.70016 23.68850 27.11065 27.19613 21.59983 16.12742 13.918166 7.284516
#  [6,] 109.50 59.48  66.46  75.85 49.47  34.02 111.89  89.56 247.88  41.20 50.20 73.78  4.259839 4.679286  6.018387 14.15200 17.27065 21.61717 25.55887 24.06065 18.78567 12.23694 10.831834 4.821613
#  [7,] 132.25 37.25 104.77 109.69 48.23  64.80 140.41 120.02 435.49  80.48 41.18 66.32  6.521451 6.417857  7.613387 14.63000 18.95694 22.39833 26.39371 25.57113 20.08700 13.71323 12.386833 5.759677
#  [8,] 151.15 80.34  83.69  44.25 42.11 175.94  87.58  64.54  52.96  55.55 83.27 67.00  8.957742 9.158393  9.275000 17.73917 19.22306 23.57900 26.13306 26.92597 21.46650 16.36113 12.618999 6.779355
#  [9,] 201.49 44.98  62.21 107.74 67.97 111.07 126.18 125.60 559.77 171.96 60.87 30.00  8.719032 8.354285  9.491936 17.15333 19.84419 23.75117 26.98613 27.21161 21.84600 16.34677 14.269667 7.459194
# [10,] 181.53 53.87  50.79  77.73 74.72 115.68 142.59  95.06 682.20 172.46 42.90 24.78 10.136129 9.736964 10.580484 17.32333 20.06032 24.10600 27.43500 26.79452 22.09217 17.08403 14.184667 8.277097

# Warning message:
# In st_extract.stars(bcsd_obs, at = pnt_mat) :
#   incompatible units merged in result matrix: dropping all units

原始答案（为确保准确性更新）

当将矩阵提供给 at 参数于 st_extract() 时，会返回一个矩阵。矩阵中的数值必须为单一类型，而bcsd_obs同时包含Celsius（C）和mm/m（毫米/米）。然而，"n error(s) in InterpolateAtPoint()" 是GDAL的一个错误(?)，以及它处理NA的方式。请注意，使用矩阵返回的错误数量（在下方成功结果中，带有 set.seed(42) 的情况为36）与下方的NA数量相对应。

这并不仅限于在 st_extract() 使用矩阵对象，或在R 中使用。请参阅这个 GDAL GitHub问题，以及这个 stars GitHub issue。

可以通过使用sf对象而不是矩阵来避免此错误。尽管会返回结果，如同这个 stars dev comment 所示，您仍会看到两个（无害的）错误信息。在您的情况下，这是因为 class(bcsd_obs) 是一个 "stars_proxy"/"stars" 对象。

此外，尽管承认bcsd_obs被用于一个最小可复现实例（MRE），为了完整性我还要指出的是bcsd_obs没有CRS。在处理空间数据时，这通常并不理想。

为解决这些问题，您可以可选地将bcsd_obs转换为仅包含 "stars" 的对象，并在运行 st_extract() 之前给bcsd_obs指定一个CRS：

library(stars) # Version 0.7.2
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE

set.seed(42) # For reproducibility

# Assign CRS (assumed to be WGS84/EPSG:4326, you'll need to double-check)
st_crs(bcsd_obs) <- 4326
# Convert to stars
bcsd_obs_loc <- st_as_stars(bcsd_obs)

pnt_mat <- st_sample(st_as_sfc(st_bbox(bcsd_obs_loc)), 10)
st_extract(bcsd_obs_loc, at = pnt_mat)
# stars object with 2 dimensions and 2 attributes
# attribute(s):
#                 Min.   1st Qu.   Median      Mean   3rd Qu.    Max. NA's
# foo [mm/m] 24.779999 53.642499 78.97500 110.75262 117.27500 682.200   36
# bar [C]     2.075806  8.898064 15.93347  15.33476  21.49983  27.435   36
# dimension(s):
#          from to  refsys point
# geometry    1 10  WGS 84  TRUE
# time        1 12 POSIXct    NA
#                                                             values
# geometry POINT (-75.73759 34.86242),...,POINT (-77.86122 35.28557)
# time                                     1999-01-31,...,1999-12-31