r - 使用 R 将 GRIB 数据绘制为经度超过 360 度的 map

Question

我正在尝试自动获取有关波高预报的数据并构建类似 this one 的图表.

使用 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 度的经度。

这对我正在做的事情很重要，因为它使我可以轻松地检查同一地 block 内所有海洋上的涌浪 - 如果一次只显示 360 度，这会更加困难，因为这会强制导致以下情况之一被切割的海洋。

尽管我尽了最大的努力，但我找不到一种方法来绘制超过 360 度的图，因为 GRIB 格式“太聪明”而不允许这样做(这不仅仅是偏移数据，而是重复了一部分它)。

任何见解将不胜感激。干杯

Answer 1

我会将您的数据从 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]])

首先移动和裁剪栅格，然后合并可能更有效，但这只是一个开始，只需要几秒钟就可以在栅格上运行。

如果您只关心绘图，那么编写一个函数来绘制它两次可能会更容易，一次是在范围内移动。但这必须对光栅中的每个波段进行....

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 包的所有功能。