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

Question

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

下面是我的netCDF文件的摘要。它有两个维度（纬度和经度）和一个变量（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˚，长：352˚，右上角：纬：44.5˚

我知道可以使用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.subset(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 

我在上面的讨论中注意到有一个关于纬度顺序的问题。您也可以使用数位长命令“invertlat”来为您解决这个问题。

Answer 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或使用哪个创建索引并在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)

备注:我更喜欢< code>ncdf4，我发现它的语法更容易记忆(还有一个比旧的netcdf R-package更好的优点，我已经忘记了...)

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

>

哪种功能方式将起作用。我自己使用它。

数据将采用与 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]

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

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

坐标当然与netcdf文件中的坐标相同。。。

您需要samity检查最后一个cbid，因为我只有大约98%的信心我没有弄乱坐标。在我在桌面上找到的R脚本中，我使用了循环，这是…邪恶的…这应该是（有点？）更快，也更明智。