3

波高予測に関するデータを自動的に取得し、これと同様のプロットを作成しようとしています。

R を使用して、データをダウンロードし、次を使用してプロットできます。

library(rgdal)
library(fields)

ftp.string <- "ftp://polar.ncep.noaa.gov/pub/waves//20130205.t00z/nww3.HTSGW.grb"
#this link may become broken with time, as folders are removed after some time. just edit the date to reflect the most recent day at the time you run these lines

download.file(ftp.string, "foo.grb", mode="wb")

grib <- readGDAL("foo.grb")
is.na(grib$band1) <- grib$band1 > 100
image(grib, col=(tim.colors(15)), attr=1)

ただし、上に投稿したリンクを詳しく見てみると、微妙な違いに気付くでしょう。リンクからのプロットは、経度で 360 度以上にまたがっています。

これは、同じプロット内のすべての海のうねりを簡単に確認できるため、私が行っていることにとって重要です。一度に 360 度しか表示されない場合、海の 1 つが強制的に表示されるため、これははるかに困難です。切る。

最善を尽くしたにもかかわらず、360 度以上をプロットする方法を見つけることができませんでした。これは、GRIB 形式が「賢すぎて」それを許可していないためです (これは、データをオフセットするだけでなく、その一部を繰り返すことになります)。

どんな洞察も大歓迎です。乾杯

4

2 に答える 2

5

rasterパッケージからデータをラスター スタックにロードし、関数mergeとcrop関数を使用します。基本的には、ラスターを複製し、360 度シフトし、それをそれ自体とマージしてから、好みに合わせてトリミングします。ここに関数があります：

require(raster)
wwrap <- function(g,xmax=720){
  gE = extent(g)

  shiftE = extent(360,720,gE@ymin, gE@ymax)
  g2 = g
  extent(g2)=shiftE

  gMerge = merge(g,g2)

  crop(gMerge,extent(0,xmax,gE@ymin, gE@ymax))
}

そして、ここにいくつかの使用法があります：

> gstack = stack("foo.grb")
> gstack[gstack>100] = NA
> gstack2 = wwrap(gstack,xmax=460)
> plot(gstack2)
> plot(gstack2[[1]])
> plot(gstack2[[61]])

おそらく、最初にラスターをシフトしてトリミングしてからマージする方が効率的ですが、これは最初のステップであり、ラスターで実行するのに数秒しかかかりません。

関心があるのがプロットだけである場合は、範囲をシフトして 1 回、2 回プロットする関数を作成する方が簡単かもしれません。しかし、それはラスターの各バンドに対して行う必要があります....

wraplot <- function(g,xmax=720){
  gE = extent(g)
  ## to setup the plot
  worldWrap = extent(0,xmax,gE@ymin, gE@ymax)
  rWrap = raster(nrows=1,ncols=1,ext=worldWrap)
  rWrap[]=NA
  plot(rWrap)
  ## first part
  plot(g,add=TRUE)
  ## setup and plot it again
  shiftE = extent(360,720,gE@ymin, gE@ymax)
  cropE = extent(360,xmax,gE@ymin, gE@ymax)
  extent(g)=shiftE
  g=crop(g,cropE)
  plot(g,add=TRUE)
}

次に、次のようにします。

wraplot(gstack[[1]])

rasterパッケージのすべての機能を確認してください。

于 2013-02-06T08:43:43.337 に答える
2

より単純なアプローチは、360°のオフセットで2番目のデータセットを作成することです。

grib2 <- grib
grib2@bbox[1, ] <- grib2@bbox[1, ] - 360
image(grib, col=(tim.colors(15)), attr=1, xlim=c(-360,360))
image(grib2, add=TRUE, col=(tim.colors(15)), attr=1)

ここに画像の説明を入力してください

そして、あなたはそれをxlimあなたが望むように中央に置くために遊ぶことができます：

image(grib, col=(tim.colors(15)), attr=1, xlim=c(-270,90))
image(grib2, add=TRUE, col=(tim.colors(15)), attr=1)

ここに画像の説明を入力してください

ここでは、データが通常のグリッド上にあるため、補間の必要がないため、これが機能します。それ以外の場合は、もちろん、@Spacedmanソリューションが推奨されます。

于 2013-02-07T09:19:46.883 に答える