Plotting GRIB data as a map with more than 360 degrees longitude using R
I am trying to automatically grab the wave height prediction data and build plots similar to this one .
Using R, I can load data and plot it with:
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)
However, if you take a closer look at the link posted above, you will notice a difference in subtlety: the graph from the link spans over 360 degrees in longitude.
This is important for what I am doing, as it allows me to easily check for swelling in all oceans within the same area - it is much more difficult if only 360 degrees are shown at a time, as this causes one of the oceans to be cut.
Despite my best efforts, I cannot find a way to plot more than 360 degrees, since the GRIB format is "too smart" to allow this (it is not just to compensate for the data, but to repeat part of it instead).
Any ideas would be greatly appreciated. Greetings
source to share
A more naive approach would be to create a second dataset with a 360 ° offset:
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)
And you can play with xlim
to center it the way you want:
image(grib, col=(tim.colors(15)), attr=1, xlim=c(-270,90))
image(grib2, add=TRUE, col=(tim.colors(15)), attr=1)
Here it works because the data is on a regular grid, so there is no need for interpolation, otherwise @ Spacedman's solution of course should be preferred.
source to share
I would load your data into a raster stack from a package raster
and then use the merge
and functions crop
. Basically you duplicate the raster, nudge it 360 degrees, then merge it with you, and then crop it to your liking. Here's the function:
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))
}
And here's some usage:
> gstack = stack("foo.grb")
> gstack[gstack>100] = NA
> gstack2 = wwrap(gstack,xmax=460)
> plot(gstack2)
> plot(gstack2[[1]])
> plot(gstack2[[61]])
It is probably more efficient to move and crop the raster first and then merge, but this is the beginning and it only takes a few seconds for the raster to run.
If all you care about is the conspiracy, then it would be easier to write a function that displays it twice, once with a shift. But this would have to be done for each group in your raster ...
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)
}
Then follow these steps:
wraplot(gstack[[1]])
Check out all the features of the package raster
.
source to share