10

在另一个线程中的帮助下，我设法绘制了一些全球地图。首先，我将气象 GRIB2 数据转换为 Netcdf，然后绘制全球地图。

现在我只想绘制地图的一个子区域。我尝试了crop命令并成功提取了全局nc文件的子区域。但是在绘图时我找不到如何控制轴限制。它绘制了一个比数据区域更大的地图，因此两边都会出现大的空白。

这是我用来绘制地图的脚本

library("ncdf")
library("raster")
library("maptools")

DIA=format(Sys.time(), "%Y%m%d00") # Data d'avui
url=sprintf("ftp://ftp.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/gfs.%s/gfs.t00z.pgrb2f00", DIA) # Ruta del ftp
loc=file.path(sprintf("%s",url))
download.file(loc,"gfs.grb",mode="wb")

system("/usr/bin/grib2/wgrib2/wgrib2 -s gfs.grb | grep :TMP: | /usr/bin/grib2/wgrib2/wgrib2 -i gfs.grb -netcdf temp.nc",intern=T)

t2m <- raster("temp.nc", varname = "TMP_2maboveground")
rt2m <- rotate(t2m)
t2mc=rt2m-273.15

DAY=format(Sys.time(), "%Y%m%d") # Data d'avui

e=extent(-40,40,20,90)
tt=crop(t2mc,e)

png(filename="gfs.png",width=700,height=600,bg="white")    
    rgb.palette <- colorRampPalette(c("snow1","snow2","snow3","seagreen","orange","firebrick"), space = "rgb")#colors
    plot(tt,col=rgb.palette(200),main=as.expression(paste("Temperatura a 2m ",DAY," + 00 UTC",sep="")),axes=T)
dev.off()

给出这个输出。

裁剪图

它必须是一个简单的，但我是一个简单的 R 用户。提前致谢。

编辑：按照建议添加 xlim=c(-40,40),ylim=c(20,90) 时的新输出。似乎它不能解决问题。但是使用输出 png 文件的 x,y 大小看起来很有希望，因为我可以调整大小以适应地图。肯定它必须是另一种解决方案，我找不到正确的解决方案。

在此处输入图像描述

4

1 回答 1

9

下载数据文件后，我可以直接用光栅读取。我选择了 221 波段（如果我没记错的话），根据 这张表，这就是你所需要的：

library("raster")
t2mc <- raster('gfs.grb', band=221)

> t2mc
class       : RasterLayer 
band        : 221  (of  315  bands)
dimensions  : 361, 720, 259920  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -0.25, 359.75, -90.25, 90.25  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +a=6371229 +b=6371229 +no_defs 
data source : /home/oscar/gfs.grb 
names       : gfs

您不需要整个范围，因此您可以使用它crop来获得所需的范围：

e <- extent(-40,40,20,90)
tt <- crop(t2mc,e)

我试图显示tt光栅plot但没有成功。但是，spplot如果您使用不同的范围（89.5 而不是 90），它可以正常工作：

e <- extent(-40,40,20,89.5)
tt <- crop(t2mc,e)

spplot(tt)

现在我们必须添加行政边界：

library(maps)
library(mapdata)
library(maptools)

ext <- as.vector(e)
boundaries <- map('worldHires',
                  xlim=ext[1:2], ylim=ext[3:4],
                  plot=FALSE)
boundaries <- map2SpatialLines(boundaries,
                               proj4string=CRS(projection(tt)))

并更改调色板：

rgb.palette <- colorRampPalette(c("snow1","snow2","snow3","seagreen","orange","firebrick"),
                                space = "rgb")

spplot(tt, col.regions=rgb.palette,
       colorkey=list(height=0.3),
       sp.layout=list('sp.lines', boundaries, lwd=0.5))

spplot 结果

如果您更喜欢这种latticeExtra::layer方法，您可以使用以下代码获得类似的结果：

library(rasterVis)
levelplot(tt, col.regions=rgb.palette,
          colorkey=list(height=.3)) +
    layer(sp.lines(boundaries, lwd=0.5))
于 2013-06-20T23:45:33.007 回答