如何通过经纬度索引R语言中的RasterBrick对象

嘿，刚好我之前也处理过类似的RasterBrick经纬度提取需求，给你两个亲测好用的思路，先结合你的数据背景整理下：

你的RasterBrick对象信息（方便参考）：

> x
class      : RasterBrick 
dimensions : 360, 720, 259200, 110  (nrow, ncol, ncell, nlayers)
resolution : 0.5, 0.5  (x, y)
extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref.: +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source: /Users/kirkw/R_code/rgis/test.grd 
names      : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7, layer.8, layer.9, layer.10, layer.11, layer.12, layer.13, layer.14, layer.15, ... 
min values : -57.59665, -56.89876, -57.92206, -55.95920, -57.09456, -55.45393, -55.38048, -57.27017, -56.05618, -56.92603, -57.68953, -55.23584, -56.70047, -56.50427, -56.72116, ... 
max values : 30.52231, 30.49426, 29.96689, 30.05002, 30.13025, 30.08621, 29.73001, 30.27164, 30.00519, 30.15961, 30.15723, 30.53934, 30.00043, 30.47162, 30.45111, ... 

你已经能通过行列号提取单元格的时间序列：

x[1, 1, ]
# 输出示例：
layer.1  layer.2  layer.3  layer.4  layer.5  layer.6  layer.7  layer.8  layer.9 
-14.92191 -16.13638 -14.39139 -15.55865 -14.23444 -14.64407 -13.99429 -13.81390 -14.96927 
layer.10 layer.11 layer.12 layer.13 layer.14 layer.15 layer.16 layer.17 layer.18 
-14.28351 -15.20233 -15.07458 -13.19043 -14.15805 -13.78281 -14.42188 -14.89151 -14.91458 
...

接下来给你两种可靠的经纬度提取方法：

方法1：直接用extract()函数（最直观） #

raster包的extract()函数原生支持传入经纬度坐标，直接提取对应位置的全图层时间序列，步骤很简单：

# 定义你要提取的目标经纬度，比如北京的经纬度：116.4°E，39.9°N
target_coords <- data.frame(lon = 116.4, lat = 39.9)

# 提取时间序列
ts_data <- extract(x, target_coords)

# 查看结果，会得到一个包含所有110个图层值的向量
ts_data

因为你的RasterBrick的坐标参考系（CRS）已经是WGS84（和常用的GPS经纬度一致），所以直接用原始经纬度就行，不用做坐标转换。如果你的数据CRS不一样，记得先把目标经纬度转成对应CRS再传入。

方法2：先转单元格索引再提取（和你现有逻辑衔接） #

如果你习惯用索引提取的方式，可以先通过cellFromXY()把经纬度转换成对应的单元格编号，再用你熟悉的方式提取，逻辑和你用行列号提取一致：

# 定义目标经纬度
lon <- 116.4
lat <- 39.9

# 获取该经纬度对应的单元格索引
cell_idx <- cellFromXY(x, c(lon, lat))

# 提取该单元格的时间序列
ts_data <- x[cell_idx, ]

# 查看结果
ts_data

这个方法的优势是如果需要批量提取多个位置，先批量获取单元格索引会更高效，避免重复的坐标匹配计算。

小提醒 #

• 确保你的经纬度在数据的范围里（你的数据是全球范围：经度-180到180，纬度-90到90），超出范围会返回NA
• 如果经纬度刚好落在单元格边界上，cellFromXY()会自动匹配到对应的单元格，不用你手动计算行列号，非常省心

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