如何获取指定经纬度范围的GeoTIFF高程数据或裁剪西班牙DTM文件？

我来帮你搞定这个问题！你已经拿到了覆盖全西班牙的20m分辨率DTM，现在只需要提取你在leaflet里框选的那一小块区域就行。这里分两种常用方法，都是用R工具处理，你可以选适合自己的：

方法一：用terra包（推荐，处理大文件更高效） #

terra是raster包的升级版，对大栅格文件的读写和裁剪速度更快，内存占用也更合理。步骤如下：

1. 先确保你装了terra包（没装的话先跑install.packages("terra")）
2. 执行以下代码：

library(terra)

# 读取你的全西班牙DTM文件（注意路径要正确）
dtm_full <- rast("DTM Spain_Mainland (2019) 20m.tif")

# 定义你的经纬度边界框（注意顺序是xmin, xmax, ymin, ymax）
bbox_wgs <- ext(-3.7525599, -3.6525599, 40.4065001, 40.4965001)
# 给这个边界框设置WGS84投影（也就是你leaflet用的经纬度投影）
crs(bbox_wgs) <- "EPSG:4326"

# 把经纬度边界框转换成DTM的UTM投影（Zone30，和你的DTM投影一致）
bbox_utm <- project(bbox_wgs, crs(dtm_full))

# 裁剪DTM到目标范围
dtm_cropped <- crop(dtm_full, bbox_utm)

# 把原数据的NoData值（-32767）设置正确，避免后续渲染出问题
NAvalue(dtm_cropped) <- -32767

# 保存裁剪后的小文件，方便后续用rayshader处理
writeRaster(dtm_cropped, "dtm_my_area.tif", overwrite = TRUE)

方法二：用whitebox包（基于GDAL，适合超大型文件） #

如果你觉得加载2.4GB的文件太占内存，可以用whitebox的GDAL底层工具来裁剪，不需要把整个文件读进内存：

library(whitebox)
library(terra)

# 先获取UTM格式的边界框坐标（用terra来转换）
dtm_full <- rast("DTM Spain_Mainland (2019) 20m.tif")
bbox_wgs <- ext(-3.7525599, -3.6525599, 40.4065001, 40.4965001)
crs(bbox_wgs) <- "EPSG:4326"
bbox_utm <- project(bbox_wgs, crs(dtm_full))

# 提取UTM坐标的四个值
x_min <- bbox_utm[1]
x_max <- bbox_utm[2]
y_min <- bbox_utm[3]
y_max <- bbox_utm[4]

# 用whitebox的crop_raster工具裁剪
crop_raster(
  input = "DTM Spain_Mainland (2019) 20m.tif",
  output = "dtm_my_area_wb.tif",
  x_min = x_min,
  x_max = x_max,
  y_min = y_min,
  y_max = y_max,
  nodata = -32767
)

为什么要转换投影？ #

你的DTM用的是UTM Zone30投影（单位是米），而你定义的边界框是WGS84经纬度（单位是度），直接用经纬度裁剪会因为坐标系统不匹配而出错，所以必须先把边界框转换成和DTM一致的投影。

裁剪后怎么用rayshader？ #

拿到裁剪后的dtm_my_area.tif，就可以按照你参考的博客流程继续了，比如：

library(rayshader)

# 读取裁剪后的DTM
dtm_raster <- raster("dtm_my_area.tif")
# 转换成rayshader需要的矩阵格式
elmat <- raster_to_matrix(dtm_raster)

# 生成3D地形渲染
elmat %>%
  sphere_shade(texture = "desert") %>%
  plot_3d(elmat, zscale = 20, fov = 60, theta = 25, zoom = 0.75, phi = 45)

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