如何使用R中的纬度/经度边界从netCDF文件中获取子集 [英] How to take a subset from a netCDF file using latitude/longitude boundaries in R

问题描述

我有一个netCDF文件，我希望使用R中的"ncdf"包从纬度/经度边界定义的子集(即经纬度定义的经纬度框)中提取子集.

I have a netCDF file that I wish to extract a subset from defined by latitude/longitude boundaries (i.e. a lat/long defined box), using the ‘ncdf’ package in R.

以下是我的netCDF文件的摘要.它具有二维(经度和纬度)和1个变量(10U_GDS4_SFC).它实质上是一个包含风值的经/纬网格:

A summary of my netCDF file is below. It has two dimensions (latitude and longitude) and 1 variable (10U_GDS4_SFC). It is essentially a lat/long grid containing wind values:

[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.

The latitude variable runs from +90 to -90 and the longitude variable runs form 0 to 360.

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

I wish to extract a subset of the overall grid using the following geographical corner boundaries:

左下角:纬度:34.5˚，长线:355˚， 左上角:纬度:44.5˚，长线:355˚， 右上角:纬度:44.5˚，长线:12˚， 右下角:纬度:34.5˚，长度:12˚

bottom left corner: Lat: 34.5˚, Long: 355˚, top left corner: Lat: 44.5˚, Long: 355˚, top right corner: Lat: 44.5˚, Long: 12˚, bottom right corner: Lat: 34.5˚ , Long: 12˚

我知道可以使用get.var.ncdf()命令(下面的示例)提取变量的一部分:

I am aware that parts of a variable can be extracted using the get.var.ncdf() command (example below):

z1 = get.var.ncdf(example.nc, "10U_GDS4_SFC", start=c(11,26), count=c(5,5))

但是，我不知道如何合并经/纬度，因此最终得到包含变量值的子集空间网格.我是R中使用netCDF值的新手，将不胜感激任何建议.非常感谢！

However, I can’t work out how lat/long can be incorporated so that I end up with a subsetted spatial grid containing variable values. I am new to working with netCDF values in R, and any advice would be greatly appreciated. Many thanks!

推荐答案

原则上，您的位置是那里的2/3.当然，您可以使用以下方式创建起始索引:

In principle you are 2/3 of the way there. You can of course create the starting indices using something like this:

require(ncdf4)

ncFile <- nc_open( MyNetCDF )
LonStartIdx <- which( ncFile$dim$lon$vals == 355)
LatStartIdx <- which( ncFile$dim$lat$vals == 34.5)

对计数进行相同的操作.然后，读取所需的变量

Do the same for the counts. Then, read the variable you want

MyVariable <- ncvar_get( ncFile, varName, start=c( LonStartIdx, LatStartIdx), count=...)

但是，就您而言，据我所知，您很不走运.读/写netcdf例程按顺序执行它们的工作.网格会回绕，因为坐标的经度范围是0-360，并且您对包含零子午线的框感兴趣.

However in your case you are out of luck as far as I know. The reading / writing netcdf routines do their stuff sequentially. Your grid wraps around since you have coordinates that go from 0 - 360 in longitude and you are interested in a box that contains the zero meridian.

对于您(假设您没有太多数据)，将整个网格读入R，然后使用subset或使用which创建索引并在其中切出框"会更有意义. R.

For you (assuming you have not too much data) it would make more sense to read in the full grid into R, and then use either subset or create indices using which and cut out your "box" in 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软件包相比，还有另一个优势……)

Remark: I prefer ncdf4, I find the syntax a bit easier to remember (and there was another advantage over the older netcdf R-package that I have forgotten...)

好的.注释的长度不能满足我的需要，因此我更新了答案 不用担心.让我们逐步解决问题.

Ok. Comments cannot be as long as I would need them, so I updated the answer No worries. Let's go through the questions step by step.

• which函数方法将起作用.我自己用.

• 数据将采用与netcf文件类似的格式，但是我不太确定0子午线是否存在问题(我想是的).您可能需要通过执行以下操作来交换两半(替换第二个示例中的相应行)

• 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.

可以将所有内容重新格式化为2x3数据帧.拿第二个代码示例返回的数据(将是一个矩阵[lon x lat].还要从中获取坐标的值

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]

(或示例中的经度调用方式，与lat相同).然后使用

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

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

The coordinates will of course be the same as in the netcdf file...

您需要检查最后一个绑定，因为我只有98％的信心对我没有弄乱坐标.在我在桌面上找到的R脚本中，我使用循环，这是...邪恶的...这应该更快一点，并且也更明智.

You need to samity check the last cbind, as I am only about 98% confident that I have not messed up the coordinates. In the R scripts I found on my desktop I use loops, which are... evil... This should be (a bit?) faster and is also more sensible.

这篇关于如何使用R中的纬度/经度边界从netCDF文件中获取子集的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！