基于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 layers/depth) 砖和一个测深光栅。您想获得海底的 OmegaA 值。这是我的做法：

将 OmegaA 栅格设为与测深仪相同的分辨率和范围
将测深栅格转换为33个0-1的三层栅格砖。例如如果某个特定像素的海底深度为 200 米，则除 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)

然后你将两个光栅砖组合起来。它们具有相同的维度，这应该很容易。输出应该是一个 33 层的栅格砖块，在不是海底的所有地方都为 0，在其他任何地方都为 OmegaA 的值。
将前面得到的砖块的所有图层合并成一个简单的栅格，求和。

这应该有效。如果您在处理栅格砖时遇到问题，您可以将数据制作成基本的 R 数组，它可能更简单。

祝你好运。

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

Query raster brick layer based on another raster in R

r

raster

netcdf4