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 度，这将更加困难，因为这会强制导致其中一个海洋被切。

尽管我尽了最大的努力，我还是找不到绘制超过 360 度的方法，因为 GRIB 格式“太聪明”而不允许这样做（这不仅仅是抵消数据，而是重复其中的一部分）。

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

4

2 回答 2

5

我会将您的数据从raster包中加载到栅格堆栈中，然后使用mergeandcrop函数。基本上，您复制光栅，将其移动 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软件包的所有功能。

于 2013-02-06T08:43:43.337 回答
2

更天真的方法是创建第二个数据集，偏移为 360°：

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 回答