如何使用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运行.

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

左下角:Lat:34.5˚,长:355˚,左上角:Lat:44.5˚,长:355˚,右上角:Lat:44.5˚,长:12˚,右下角:Lat: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

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

Answer 2

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

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

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