如何使用Terra包从坐标数据获取地形特征（TRI）

拥有一组坐标数据，想要获取地形特征数据（特别是TRI）。已知terra包的terrain()函数可提供此类信息，但该函数需要包含高程数据的SpatRaster对象作为输入，在创建正确的SpatRaster对象时遇到困难。

尝试先用elevatR包获取坐标对应的高程值，但用rast()函数创建SpatRaster对象时失败：

library(ctmm)
library(dplyr)
library(elevatR)
library(terra)

# 从ctmm包内置追踪数据创建示例数据
data("buffalo")
Cilla <- buffalo$Cilla
Cilla <- data.frame(Cilla)
lonlat <- Cilla %>% dplyr::select(2,3)

# 定义经纬度投影
prj_dd <- "+proj=longlat +datum=NAD83"

# 使用elevatR包获取每个点的高程
coords_E <- as.data.frame(elevatr::get_elev_point(lonlat, prj = prj_dd, src = "aws"))

# 向主数据框添加高程列（单位：米）
lonlat$elev <- coords_E$elevation

# 尝试创建SpatRaster失败
raster <- rast(lonlat, type="xyz")

疑问：哪里出错了？以及如何使用terrain()函数获取坐标数据对应的TRI值？

你的核心问题是：terrain()计算TRI需要连续的高程栅格（包含邻域单元的高程信息），而你用点数据创建的是离散点对应的栅格，每个点是独立单元，没有相邻栅格，无法计算TRI。正确流程是：

• 第一步：获取覆盖所有坐标点的连续高程栅格
• 第二步：基于该栅格计算TRI
• 第三步：将TRI值提取到对应的坐标点上

修改后的代码如下：

library(ctmm)
library(dplyr)
library(elevatr)
library(terra)

# 1. 处理示例数据
data("buffalo")
Cilla <- buffalo$Cilla
Cilla_df <- data.frame(Cilla)
lonlat <- Cilla_df %>% select(2,3)
prj_dd <- "+proj=longlat +datum=NAD83"

# 2. 将点数据转为SpatVector，确定研究区域范围
points_v <- vect(lonlat, geom=c("longitude","latitude"), crs=prj_dd)
# 给研究范围加缓冲，确保覆盖所有点及周边（缓冲距离单位是度，可按需调整）
study_ext <- ext(points_v) %>% expand(0.1, 0.1)

# 3. 获取该区域的连续高程栅格
elev_raster <- get_elev_raster(locations = study_ext, prj = prj_dd, src = "aws", z=10)
# 转为terra的SpatRaster格式
elev_rast <- rast(elev_raster)

# 4. 计算TRI地形粗糙度指数
tri_rast <- terrain(elev_rast, v="tri", unit="m")

# 5. 将TRI值提取到原始坐标点上
lonlat$TRI <- extract(tri_rast, points_v)[,2]

关键说明： #

• 用get_elev_raster()替代get_elev_point()，获取连续的高程栅格而非单个点的高程值，这样terrain()才能利用邻域栅格计算TRI。
• z参数控制栅格分辨率，值越大分辨率越高（数据量也越大），可根据研究需求调整。
• 给研究范围加缓冲expand(0.1,0.1)是为了避免边缘点没有足够邻域计算TRI，缓冲距离可根据实际情况修改。

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