如何使用 R 中的纬度/经度边界从 netCDF 文件中获取子集

Question

我有一个 netCDF 文件，我希望使用 R 中的“ncdf”包从定义的纬度/经度边界（即纬度/经度定义的框）中提取一个子集。

我的 netCDF 文件的摘要如下。它有两个维度（纬度和经度）和 1 个变量（10U_GDS4_SFC）。它本质上是一个包含风值的纬度/经度网格：

[1] "file example.nc has 2 dimensions:"
[1] "lat_0   Size: 1280"
[1] "lon_1   Size: 2560"
[1] "------------------------"
[1] "file example.nc has 1 variables:"
[1] "float 10U_GDS4_SFC[lon_1,lat_0]  Longname:10 metre U wind component Missval:1e+30"

纬度变量从+90到-90，经度变量从0到360。

我希望使用以下地理角边界提取整个网格的子集：

左下角：纬度：34.5˚，经度：355˚， 左上角：纬度：44.5˚，经度：355˚， 右上角：纬度：44.5˚，经度：12˚， 右下角：纬度：34.5˚，经度：12˚

我知道可以使用get.var.ncdf() 命令提取变量的一部分（示例如下）：

z1 = get.var.ncdf(example.nc, "10U_GDS4_SFC", start=c(11,26), count=c(5,5))

但是，我无法计算出如何合并纬度/经度，因此我最终得到了一个包含变量值的子集空间网格。我是在 R 中使用 netCDF 值的新手，任何建议都将不胜感激。非常感谢！

Answer 1

如果您使用的是 Linux，这可以使用 nctoolkit (https://nctoolkit.readthedocs.io/en/latest/) 轻松实现：

import nctoolkit as nc
data = nc.open_data("example.nc")    
data.clip(lon = [-12, -5], lat = [35.4, 44.5])

Answer 2

也可以使用CDO先从bash命令行中提取区域，然后读取R中的文件：

cdo sellonlatbox,-5,12,34.5,44.5 in.nc out.nc 

我注意到在上面的讨论中存在关于纬度顺序的问题。您还可以使用 CDO 命令“invertlat”为您解决问题。

Answer 3

原则上你是那里的 2/3。您当然可以使用以下方式创建起始索引：

require(ncdf4)

ncFile <- nc_open( MyNetCDF )
LonStartIdx <- which( ncFile$dim$lon$vals == 355)
LatStartIdx <- which( ncFile$dim$lat$vals == 34.5)

对计数执行相同的操作。然后，读取你想要的变量

MyVariable <- ncvar_get( ncFile, varName, start=c( LonStartIdx, LatStartIdx), count=...)

但是据我所知，在你的情况下你很不走运。读/写 netcdf 例程按顺序进行。您的网格环绕，因为您的坐标从经度 0 到 360，并且您对包含零子午线的框感兴趣。

对您来说（假设您没有太多数据）将完整的网格读入 R 会更有意义，然后使用 subset 或使用 which 创建索引并在R.

ncFile <- nc_open( MyNetCDF )
LonIdx <- which( ncFile$dim$lon$vals > 355 | ncFile$dim$lon$vals < 10)
LatIdx <- which( ncFile$dim$lat$vals > 34.5 & ncFile$dim$lat$vals < 44.5)
MyVariable <- ncvar_get( ncFile, varName)[ LonIdx, LatIdx]
nc_close(ncFile)

备注：我更喜欢ncdf4，我发现语法更容易记住（与我忘记的旧netcdf R-package相比还有另一个优势......）

好的。评论不能像我需要的那样长，所以我更新了答案 不用担心。让我们一步一步地回答问题。

• which 函数方式将起作用。我自己用。

• 数据的格式与 netcf 文件中的格式相似，但我不太确定 0 子午线是否有问题（我猜是）。您可能必须通过执行以下操作来交换两半（替换第二个示例中的相应行）

  LonIdx <- c(which( ncFile$dim$lon$vals > 355) , which( ncFile$dim$lon$vals < 10) )
  

  这会改变坐标索引的顺序，使西部先出现，然后是东部。

• 可以将所有内容重新格式化为 2x3 数据框。获取我的第二个代码示例返回的数据（将是一个矩阵，[lon x lat]。还从

  获取坐标值

  lon <- ncFile$dim$lon$val[LonIdx]
  

  （或在您的示例中如何调用经度，lat 相同）。然后使用

  组装矩阵

  cbind( rep(lat, each=length(lon)), rep(lon,length(lat)), c(myVariable) )
  

• 坐标当然会和netcdf文件中的一样……

您需要仔细检查最后一个 cbind，因为我只有大约 98% 的把握确定我没有弄乱坐标。在我在桌面上找到的 R 脚本中，我使用循环，它们是......邪恶的......这应该（有点？）更快，也更明智。