如何生成空间矩形坐标组合分块下载EMODnetWCS地形数据？

问题背景 #

因内存限制，计划分块从EMODnetWCS下载地形数据后合并。目标范围为经度-15至13，纬度35至63，尝试用循环将区域划分为1°×1°矩形下载，但仅获取到对角线区域，需实现所有矩形的正确坐标定义以完成批量下载。

错误原因 #

原代码将经度序列x和纬度序列y直接合并为数据框，循环时仅取对应位置的x、y值生成边界框，这相当于只遍历了对角线方向的块，没有覆盖所有x区间与y区间的组合，因此只下载了对角线区域的数据。

修正后的代码 #

library(EMODnetWCS)
library(raster)

# 初始化WCS客户端
wcs <- emdn_init_wcs_client(service = "bathymetry")

# 定义目标范围的经纬度序列
x_seq <- seq(-15, 13, 1)  # 经度区间端点
y_seq <- seq(35, 63, 1)   # 纬度区间端点

my_list <- list()
index <- 1

# 嵌套循环遍历所有1°×1°的矩形块
for (i in 1:(length(x_seq)-1)) {
  xmin <- x_seq[i]
  xmax <- x_seq[i+1]
  
  for (j in 1:(length(y_seq)-1)) {
    ymin <- y_seq[j]
    ymax <- y_seq[j+1]
    
    # 构造当前块的边界框
    bbox <- c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax)
    cat("正在下载块：", bbox, "\n")
    
    # 下载当前块的地形数据
    cov <- emdn_get_coverage(wcs, 
                             coverage_id = "emodnet__mean",
                             bbox = bbox,
                             nil_values_as_na = TRUE)
    
    # 转换为raster对象并存入列表
    rast <- raster(cov)
    my_list[[index]] <- rast
    index <- index + 1
  }
}

# 合并所有分块数据
merged_raster <- do.call(merge, my_list)

代码说明 #

• 使用嵌套循环：外层循环遍历所有经度区间，内层循环遍历所有纬度区间，确保覆盖目标范围内的每一个1°×1°矩形块。
• 新增index变量：用于在列表中依次存储每个下载的栅格对象，避免嵌套循环的索引冲突。
• 最后通过do.call(merge, my_list)将所有分块栅格合并为完整的地形数据。

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