4

我正在尝试使用 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 年 9 月 1 日至2002年 8 月 31日期间具有网格 SWH 值。每个波段指的是一对纬度/经度，并且在每天的 00h 具有一系列 16346 SWH 值（1 波段 = 特定纬度/经度的 16346 个值）。

myframe 的尺寸为 16438 x 91。注意 91 = 7rows x 13columns。而16438这个数字几乎等于乐队的数量。额外的 2 行/波段是 long 和 lat 值，所有其他列应该是对应于 16436 波段的波高。

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

此外，即使 myframe[16438,3]=43.75 (lat) 和 myframe[16437,3]= 331.25（长）。为什么是这样？我想知道 myframe[i,j] 值实际对应于哪个纬度/经度，或者在此过程中是否存在一些数据导入错误。我假设 Degrib 没有错误。

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

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

请帮忙。

4

0 回答 0