我正在尝试使用R函数从ECMWF ERA-40网站读取带有波高的GRIB文件wavedata.grib . 这是我的源代码，直到现在：

mylat = 43.75
mylong = 331.25

# read the GRIB file
library(rgdal)
library(sp)
gribfile<-"wavedata.grib"
grib <- readGDAL(gribfile)

summary = GDALinfo(gribfile,silent=TRUE)
save(summary, file="summary.txt",ascii = TRUE)
# >names(summary): rows columns bands ll.x ll.y res.x res.y oblique.x oblique.y
rows = summary[["rows"]]
columns = summary[["columns"]]
bands = summary[["bands"]]

# z=geometry(grib)
# Grid topology:
# cellcentre.offset cellsize cells.dim
# x            326.25      2.5        13
# y             28.75      2.5         7
# SpatialPoints:
#           x     y
# [1,] 326.25 43.75
# [2,] 328.75 43.75
# [3,] 331.25 43.75

myframe<-t(data.frame(grib))

# myframe[bands+1,3]=331.25 myframe[bands+2,3]=43.75
# myframe[1,3]=2.162918 myframe[2,3]=2.427078 myframe[3,3]=2.211989
# These values should match the values read by Degrib (see below)

# degrib.exe wavedata.grib -P -pnt 43.75,331.25 -Interp 1 > wavedata.txt
# element, unit, refTime, validTime, (43.750000,331.250000)
# SWH, [m], 195709010000, 195709010000, 2.147
# SWH, [m], 195709020000, 195709020000, 2.159
# SWH, [m], 195709030000, 195709030000, 1.931

lines = rows * columns
mycol = 0
for (i in 1:lines) {
if (mylat==myframe[bands+2,i] & mylong==myframe[bands+1,i]) {mycol = i+1}
}
# notice mycol = i+1 in order to get values in column to the right

myvector <- as.numeric(myframe[,mycol])

sink("output.txt")
cat("lat:",myframe[bands+2,mycol],"long:",myframe[bands+1,mycol],"\n")
for (i in 1:bands) { cat(myvector[i],"\n") }
sink()

wavedata.grib文件在1957-09-01至2002-08-31期间记录了SWH值 . 每个频带指的是一对纬度/经度，并且在每天的00h处具有一系列16346个SWH值（1个频带=在特定纬度/长度处的16346个值） .

myframe的尺寸为16438 x 91.注意91 = 7rows x 13列 . 数字16438几乎等于频段数 . 额外的2行/波段是长和纬度值，所有其他列应该是对应于16436波段的波高 .

问题是我想在lat / long = 43.75,331.25处提取SWH（波形高度），但是它们与我在同一纬度/经度上使用Degrib实用程序读取文件的值不匹配 .

此外，我想要的正确值（2.147,2.159,1.931，...）在第4列而不是myframe的第3列，即使myframe [16438,3] = 43.75（lat）和myframe [16437,3] = 331.25（长） . 为什么是这样？我想知道myframe [i，j]值实际上对应于哪个lat / long，或者进程中是否存在某些数据导入错误 . 我假设Degrib没有错误 .

如果我想在网格点之间提取值，是否有任何R例程可以在矩阵中轻松插值？更一般地说，我需要帮助编写一个有效的R函数来提取波高，如下所示：

SWH <- function (latitude, longitude, date/time)

请帮忙 .

r function grib