基于R中的另一个栅格查询栅格砖层

Question

我有一个全球海洋学 (OmegaA) 数据的 NetCDF 文件，空间分辨率相对较粗，有 33 个深度级别。我还有一个分辨率更高的全球测深栅格。我的目标是使用从 NetCDF 文件中获取海底 OmegaA 数据，使用测深数据来确定所需的深度。到目前为止我的代码；

    library(raster)
    library(rgdal)
    library(ncdf4)

# Aragonite data. Defaults to CRS WGS84
ncin <- nc_open("C:/..../GLODAPv2.2016b.OmegaA.nc")
ncin.depth <- ncvar_get(ncin, "Depth")# 33 depth levels  

omegaA.brk  <- brick("C:/.../GLODAPv2.2016b.OmegaA.nc")
omegaA.brk  <-rotate(omegaA.bkr)# because netCDF is in Lon 0-360. 

# depth raster.  CRS WGS84
r<-raster("C:/....GEBCO.tif") 

# resample the raster brick to the resolution that matches the bathymetry raster
omegaA.brk  <-resample(omegaA.brk, r, method="bilinear")

# create blank final raster 
omegaA.rast         <- raster(ncol = r@ncols, nrow = r@nrows)
extent(omegaA.rast) <- extent(r)
omegaA.rast[]       <- NA_real_

#  create vector of indices of desired depth values
depth.values<-getValues(r)
depth.values.index<-which(!is.na(depth.values))


# loop to find appropriate raster brick layer, and extract the value at the desired index, and insert into blank raster

for (p in depth.values.index) { 
  dep.index         <-which(abs(ncin.depth+depth.values[p]) == min(abs(ncin.depth+depth.values[p]))) ## this sometimes results in multiple levels being selected

  brk.level         <-omegaA.brk[[dep.index]]  # can be more than on level if multiple layers selected above. 
  omegaA.rast[p]    <-omegaA.brk[[1]][p] ## here I choose the first level if multiple levels have been selected above

  print(paste(p, "of", length(depth.values.index))) # counter to look at progress. 
}

问题：结果是一个栅格，其中应该有数据的地方存在巨大的差距 (NA)。间隙通常具有独特的形状——例如，沿着轮廓，或沿着长直线。我粘贴了一个裁剪的示例。

enter image description here

我认为这可能是因为 1) 出于某种原因，循环中的“which”语句未找到匹配项或 2) 使用“Rotate”时可能会发生我读过的投影错位.

我已尝试确保所有范围、分辨率、单元格数量和 CRS 都相同，它们看起来确实如此。

为了加快这个过程，我将全局砖块和深海栅格裁剪到我感兴趣的区域，再次检查所有空间分辨率等是否匹配 - 为了简单起见，我没有在此处包含这些步骤。

不知所措。欢迎任何帮助！

Answer 1

没有可重现的例子，这类问题很难解决。我不知道你的问题出在哪里，但我会向你介绍我会尝试的方法。也许是好的，也许是坏的，我不知道，但它可能会激励你找到解决问题的方法。

据我了解，您有一块 OmegaA 砖（33 层/深度）和一个测深栅格。你想得到海底的 OmegaA 值。以下是我的做法：

1. 使 OmegaA 栅格具有与水深测量相同的分辨率和范围
2. 将测深栅格转化为33个0-1三层的栅格砖。例如如果某个特定像素的海底位于 200m 处，则除 200 之外的所有深度层上的该像素为 0 和 200 的 1。要对此进行编程，我会走很长的路，例如

：

r_1 <- r
values(r_1) <- values(r)==10  # where 10 is the depth (it could be a range with < or >)
r_2 <- r
values(r_2) <- values(r)==20
...
r_33 <- r
values(r_33) <- values(r)==250
r_brick <- brick(r_1, r_2, ..., r_33)

1. 然后你将你的两个光栅砖都复数。它们具有相同的维度，应该很容易。输出应该是一个 33 层的栅格砖，在不是海底的地方都是 0，在其他地方都是 OmegaA 的值。
2. 将之前获得的砖块的所有层组合成一个简单的栅格和一个总和。

这应该可行。如果您在处理光栅砖时遇到问题，您可以将数据放入基本 R 数组中，这样会更简单。

祝你好运。