Построение данных GRIB в виде карты с долготой более 360 градусов с использованием R

Я пытаюсь автоматически получить данные о прогнозе высоты волны и построить графики, подобные этому.

Используя 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 "слишком умный", чтобы допустить это (это не просто компенсирует данные, а вместо этого повторяет их часть).

Любые идеи будут с благодарностью. ура

Источник

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 предпочтительнее.

Источник

Я бы загрузил ваши данные в растровый стек из 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 пакет.

Источник
Другие вопросы по тегам