使用R将GRIB数据绘制为具有超过360度经度的地图

Question

我正在尝试自动获取波浪高度预测的数据并构建类似于此的图.

使用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度,这就更难了,因为这会强制导致其中一个海洋切.

尽管我付出了最大的努力,但我找不到绘制超过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包装的所有功能.