R：从深度级别的多维数组中获取特定深度

Question

我需要根据 HYCOM 数据计算 25 度等温线深度，该数据在选定区域具有 33 个海洋深度级别的温度。

我使用下面的 netcdf 子集工具下载数据 link

http://ncss.hycom.org/thredds/ncss/GLBa0.08/latest/temp?var=temperature&north=25&west=74.1199&east=434.1199&south=-15&horizStride=1&time_start=2016-09-27T00%3A00%3A00Z&time_end=2016-09-27T23%3A59%3A00Z&timeStride=1&vertCoord=&accept=netcdf4

数据集为 netcdf 4 格式并通过 ncdf4 库导入到 R

library(ncdf4)
ncdata <- nc_open(file)
lon <- ncvar_get(ncdata, "Longitude")
lat <- ncvar_get(ncdata, "Latitude")
temp <-ncvar_get(ncdata,"temperature")
str(temp)

num [1:4500, 1:512, 1:33] 24.7 24.6 24.6 24.7 24.7 ...

如何从上面的数组中找到特定温度（25）的深度？然后将其子集化为小区域？

Answer 1

扩展@Richard Telford 的回复，您只需在所有像素的 z 维度上应用一个函数来计算温度与 25 度相交的深度。阈值。

非常简单，然后您可以这样做：

file <- "http://ncss.hycom.org/thredds/ncss/GLBa0.08/latest/temp?var=temperature&north=25&west=74.1199&east=434.1199&south=-15&horizStride=1&time_start=2016-09-27T00%3A00%3A00Z&time_end=2016-09-27T23%3A59%3A00Z&timeStride=1&vertCoord=&accept=netcdf4"
savefile = tempfile(, fileext = "nc4")
download.file(file, savefile)
library(ncdf4)
ncdata <- nc_open(savefile)
lon <- ncvar_get(ncdata, "Longitude")  
lat <- ncvar_get(ncdata, "Latitude") 
temp <-ncvar_get(ncdata,"temperature")

temp <- temp [,,1:10] # subset depths to speed up
depths <- 1:10  # let's define some dummy depths - you want to put actual values, here !

finddepth = function(pixtemp, ...) {
  if (max(pixtemp, na.rm = TRUE) < predtemp$temp) {
    NA    # set to NA if no values >= 25
  } else {
    depth <- tryCatch({
      depth <- approx(pixtemp, depths,predtemp$temp)$y # interpolate using linear (faster)
      # interp  <- loess(depths~pixtemp)  # interpolate using loess  (slower - deals with non-linearity)
      # depth  <- predict(interp, predtemp$temp) # find depth @ temperature
      return(depth)   # send back computed depth
    }, error = function(e) {NA}
    )

  }
}
predtemp  <- data.frame(temp = 25)   # set desired isotherm
iso_depth <- apply(temp, c(1, 2), FUN = finddepth)

哪个（我认为）在 "iso_depth":

中给出了所需的数据

library(lattice)
levelplot(iso_depth, main = "Isotherm depth", col.regions = terrain.colors(250))

（图中白色区域对应点@，因为最高温度<25°，不存在25°等温线）

在这里，我通过 "approx" 使用线性插值来找到深度 @，其中 T > 25 的最后一个深度与 T < 25 的第一个 @ 之间的直线与 T = 25 级别相交.如果线性插值不适合您，您可以取消注释 "depth <- approx(...)" 之后的两行，您将使用 loess 局部二次插值（但是速度较慢并提供更多 NA）。

请注意，为了获得 "meaningful" 值，您必须用正确的深度值替换我设置的虚拟 "depths" 变量。

。另请注意，这非常慢：更复杂的方法可以提供更快的速度。您可以通过实施并行处理来提高速度。

@sudheera 之前的迭代崩溃了，因为有些像素只有 3 个非 NA 值，并且 try catch 构造有错误。

HTH

洛伦佐

R：从深度级别的多维数组中获取特定深度

R: obtain a specific depth from multi dimension array of depth levels

r

multidimensional-array

interpolation

netcdf

r-raster