Я пытаюсь прочитать файл GRIB wavedata.grib с высотой волн с веб-сайта ECMWF ERA-40, используя функцию R. Вот мой исходный код до сих пор:
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 содержит сеточные значения SWH в период с 1957-09-01 по 2002-08-31. Каждая полоса относится к паре широта / лонг и имеет серию 16346 значений SWH в 00 часов каждого дня (1 полоса = 16346 значений при определенном широте / лонге).
myframe имеет размеры 16438 x 91. Примечание 91 = 7 рядов x 13 столбцов. А число 16438 практически равно количеству полос. Дополнительные 2 строки / полосы имеют значения long и lat, все остальные столбцы должны иметь высоту волны, соответствующую 16436 полосам.
Проблема в том, что я хочу извлечь SWH (высоту волны) по широте / долготе = 43,75,331.25, но они не совпадают со значениями, которые я получаю, читая файл с помощью утилиты Degrib на том же широте / длине.
Кроме того, нужные значения (2.147, 2.159, 1.931, ...) находятся в столбце 4, а не в столбце 3 myframe, хотя myframe [16438,3] = 43,75 (lat) и myframe [16437, 3] = 331,25 (длинный). Почему это? Я хотел бы знать, какому значению широты / долготы соответствуют значения myframe [i, j] или есть какая-то ошибка импорта данных в процессе. Я предполагаю, что Дегриб не имеет ошибок.
Есть ли подпрограмма R для простой интерполяции значений в матрице, если я хочу извлечь значения между точками сетки? В целом, мне нужна помощь в написании эффективной функции R для извлечения высоты волны, например:
SWH <- function (latitude, longitude, date/time)
Пожалуйста, помогите.