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

Posted 2023-03-11

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

【英文标题】：How to take a subset from a netCDF file using latitude/longitude boundaries in R

【发布时间】：2014-02-12 08:15:38

【问题描述】：

我有一个 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 值的新手，任何建议都将不胜感激。非常感谢！

【参考方案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])

【参考方案2】：

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

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

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

【参考方案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 脚本中，我使用循环，它们是......邪恶的......这应该（有点？）更快，也更明智。

【讨论】：

@Joe 感谢您的回复。是的，0-360 经度是个问题。我需要每个变量值都有一个经纬度参考，所以我不认为你的哪个选项会起作用？您是否知道以与原始 netCDF 相同的格式或在具有以下列的 ?x3 数据帧中结束该框的方法：lat、long、variable。我还需要该框在坐标值中具有与原始 netCDF 相同的间距，因为它是一个网格。任何进一步的建议将不胜感激。为我有限的编码知识道歉。

我试图将所有内容都放入评论中，但它们的长度不够长。所以我更新了答案。希望有帮助！

@Joe 感谢您的更新。我可以看到这应该如何运作良好。我目前无法让它与我的文件一起使用，但我认为这可能与文件有关，而不是与代码有关！它给了我一个很好的起点！非常感谢

你卡在哪里了？也许我可以澄清一下答案。祝你好运！

@乔：谢谢！好像是纬度的问题。纬度数据从 +90 到 -90 运行，即使在使用 >34.5 和