R 関数を使用して、ECMWF ERA-40 Web サイトから波高を含む 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 値がグリッド化されています。各バンドは緯度/経度のペアを参照し、各日の 00 時に一連の 16346 SWH 値を持ちます (1 バンド = 特定の緯度/経度での 16346 値)。
myframe のサイズは 16438 x 91 です。91 = 7 行 x 13 列であることに注意してください。そして、16438という数字はバンド数とほぼ同じです。追加の 2 つの行/バンドは経度と緯度の値で、他のすべての列は 16436 バンドに対応する波高でなければなりません。
問題は、緯度/経度 = 43.75,331.25 で SWH (波の高さ) を抽出したいのですが、同じ緯度/経度で Degrib ユーティリティを使用してファイルを読み取った値と一致しません。
また、必要な正しい値 (2.147、2.159、1.931、...) は myframe の列 3 ではなく列 4 にあります。 331.25 (長い)。どうしてこれなの？myframe[i,j] の値が実際にどの緯度/経度に対応しているか、またはプロセス中にデータ インポート エラーがあるかどうかを知りたいです。Degrib にはエラーがないと仮定しています。
グリッドポイント間の値を抽出したい場合、マトリックス内の値を簡単に補間する R ルーチンはありますか? より一般的には、次のような波の高さを抽出する効果的な R 関数を作成する際に助けが必要です。
SWH <- function (latitude, longitude, date/time)
助けてください。