2

Following the excelent post from Joe Wheatley (http://joewheatley.net/ncep-global-forecast-system/) I managed to produce a temperature global map. But, instead of only plotting coastline I've tried to use maptools package to plot country borders. The problem comes when only eastern hemisphere country borders are plotted. I should be missing something I can't figure out, still looking for on stackoverflow and google. Hope you can help.

Here is the code I'm using (most coming from Joe's post)

loc=file.path("ftp://ftp.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/gfs.2013052100/gfs.t00z.sfluxgrbf03.grib2")
download.file(loc,"temp.grb",mode="wb")

system("wgrib2 -s temp.grb | grep :TMP: | wgrib2 -i temp.grb -netcdf TMP.nc",intern=T)
system("wgrib2 -s temp.grb | grep :LAND: | wgrib2 -i temp.grb -netcdf temp.nc",intern=T)

library(ncdf)
landFrac <-open.ncdf("LAND.nc")
lon <- get.var.ncdf(landFrac,"longitude")
lat <- get.var.ncdf(landFrac,"latitude")

temp=open.ncdf("TMP.nc")
t2m.mean <- get.var.ncdf(temp,"TMP_2maboveground")


library("fields")
library("sp", lib.loc="/usr/lib/R/site-library")
library("maptools", lib.loc="/usr/lib/R/site-library")

day="DIA"

png(filename="gfs.png",width=1215,height=607,bg="white")

rgb.palette <- colorRampPalette(c("snow1","snow2","snow3","seagreen","orange","firebrick"), space = "rgb")#colors
image.plot(lon,lat,t2m.mean,col=rgb.palette(200),main=as.expression(paste("GFS 24hr Average 2M Temperature",day,"00 UTC",sep="")),axes=T,legend.lab="o C")
data(wrld_simpl)
plot(wrld_simpl, add = TRUE)

dev.off()

and this is the image produced

enter image description here

This is a global map, should I use xlim and ylim in image.plot to extract a region (i.e. Europe)

EDIT: Added url for temp.nc file

http://ubuntuone.com/29DKAeRjUCiCzLblgfSLc9

Any help would be appreciated, thanks

pacomet
  • 5,011
  • 12
  • 59
  • 111

2 Answers2

9

The data set wrld_simpl is aligned on longitudes [-180, 180] but your grid data is clearly [0,360]. Since you have maptools loaded try this to add a modified copy to your plot: 

plot(elide(wrld_simpl, shift = c(360, 0)), add = TRUE)

There are other tools to crop and move/combine and recentre data like this, including ?rotate in the package raster. It really depends on what you need.

Another graphics-alternative is to just use "world2" in the maps package

library(maps)
map("world2", add = TRUE)

Either of those will work to finish your plot started above.

Here's another discussion related to this: Fixing maps library data for Pacific centred (0°-360° longitude) display

Community
  • 1
  • 1
mdsumner
  • 29,099
  • 6
  • 83
  • 91
  • Excellent answer, as usual. – Josh O'Brien May 21 '13 at 12:43
  • Hi, does `plot(elide(wrld_simpl, shift = c(360, 0)), add = TRUE)` substitute `plot(wrld_simpl, add = TRUE)`? Thanks. – pacomet May 21 '13 at 13:19
  • Almost perfect solution. But now, borders are plottted outside the map, just on the right side over the legend. Europe borders appear not just on the left side, inside the map, but to the right. – pacomet May 21 '13 at 13:22
  • It's far from perfect, it's a couple of workarounds. I would get a decent map in the range [0,360] - there are tools in rgdal or maptools to read shapefiles, for example, or work on reconfiguring the raster to [-180, 180]. Using the raster package instead would be a good start, it simplifies a lot. If you can post up your NetCDF file I would elaborate, but I don't have access to wgrib2. – mdsumner May 21 '13 at 20:51
  • You can find the data file, temp.nc, at http://ubuntuone.com/29DKAeRjUCiCzLblgfSLc9 also added the url in my question. I will take a look to raster package, thanks again. – pacomet May 22 '13 at 07:23
3

You can do this:

library(raster)
r <- raster("temp.nc")
r <- rotate(r)
plot(r)
mdsumner
  • 29,099
  • 6
  • 83
  • 91
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63