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

Question

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

rsubsetlatitude-longitudenetcdfcdo-climate

15

我有一个netCDF文件，希望使用R语言中的“ncdf”包从其中提取经纬度边界定义的子集(即经纬度定义的框)，该文件摘要如下。它具有两个维度（纬度和经度）和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这些数值，非常欢迎您提供任何建议。非常感谢！

- Emily

3个回答

5

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

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

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

- Adrian Tompkins

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])

- Robert Wilson

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Joe W · Accepted Answer

原则上你已经完成了三分之二的工作。当然，你可以使用以下代码创建起始索引：

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包还有其他优势...）

好的。评论不能像我需要的那么长，所以我更新了答案。没关系。让我们逐步回答这些问题。

The which function way will work. I use it myself.
The data will be in a similar format as in the netcf file, but I am not too sure if there is some problem with the 0 meridian (I guess yes). You might have to swap the two halves by doing something like this (replace the corresponding line in the 2nd example)
```
LonIdx <- c(which( ncFile$dim$lon$vals > 355) , which( ncFile$dim$lon$vals < 10) )
```
This changes the order of the coordinate indices so that the Western part comes first and then the Eastern.
Reformatting everything to a 2x3 data frame is possible. Take the data my 2nd code example returns (will be a matrix, [lon x lat]. Also get the values of the coordinates from
```
lon <- ncFile$dim$lon$val[LonIdx]
```
(or how longitude is called in your example, same for lat). Then assemble the matrix using
```
cbind( rep(lat, each=length(lon)), rep(lon,length(lat)), c(myVariable) )
```
The coordinates will of course be the same as in the netcdf file...

您需要对最后的cbind进行一致性检查，因为我只有大约98%的把握没有搞乱坐标。在我桌面上找到的R脚本中，我使用了循环，这是...邪恶的... 这种方法应该（有点？）更快，也更合理。