13

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

4

3 回答 3

8

原则上你是那里的 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 脚本中，我使用循环，它们是......邪恶的......这应该（有点？）更快，也更明智。

于 2014-01-22T11:02:05.130 回答
4

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

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

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

于 2017-11-28T03:32:40.897 回答
0

如果您使用的是 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])
于 2020-08-09T18:27:47.360 回答