R crop raster data and set axis limits

问题:

With your help in another thread I have managed to plot some global maps. First I convert meteorological GRIB2 data to Netcdf and then plot the global maps.

Now I want to plot just a subregion of the map. I have tried crop command and succesfully extracted the subregion of the global nc file. But when plotting I can't find how to control axis limits. It plots a map bigger than data region so big white spaces appear on both sides.

This is the script I'm using to plot maps

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()

that give this output.

It has to be a simple one but I am a simple R user. Thanks in advance.

EDIT: New output when adding xlim=c(-40,40),ylim=c(20,90) as suggested. It seems it does not fix the problem. But playing with x,y size of the output png file looks promising as I can adjust size to fit the map.For sure it has to be another solution, the right one I can't find.

回答1:

After downloading the data file, I can read directly with raster. I choose band 221 that (if I am not wrong) it is what you need according to this table:

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 

You don't need the whole extent so you use crop to get the desired extent:

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

I have tried to display the tt raster with plot without success. However, it works correctly with spplot if you use a different extent (89.5 instead of 90):

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

spplot(tt)

Now we have to add the administrative boundaries:

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)))

and change the palette:

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))

If you prefer the latticeExtra::layer approach, you can achieve a similar result with this code:

library(rasterVis)
levelplot(tt, col.regions=rgb.palette,
          colorkey=list(height=.3)) +
    layer(sp.lines(boundaries, lwd=0.5))