根据 Joe Wheatley ( http://joewheatley.net/ncep-global-forecast-system/ ) 的精彩帖子,我设法制作了一张全球温度图。但是,我尝试使用 maptools 包来绘制国家边界,而不是只绘制海岸线。当只绘制东半球国家边界时,问题就来了。我应该错过一些我无法弄清楚的东西,仍在stackoverflow和谷歌上寻找。希望你能帮忙。
这是我正在使用的代码(大部分来自 Joe 的帖子)
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()
这是产生的图像
这是一张全球地图,我应该在 image.plot 中使用 xlim 和 ylim 来提取区域(即欧洲)
编辑:为 temp.nc 文件添加了 url
http://ubuntuone.com/29DKAeRjUCiCzLblgfSLc9
任何帮助将不胜感激,谢谢