R裁剪光栅数据和设置轴限制 [英] R crop raster data and set axis limits

问题描述

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

现在我只想画一个地图的一个子区域。我已经尝试了crop命令并成功地提取了全局nc文件的子区域。但是当绘图我找不到如何控制轴限制。它绘制了一个比数据区域更大的地图，这样大的空格就出现在两边。

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

$ b $
库（raster）
库（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 =格式（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，snow 2，snow3，seagreen，orange，firebrick），space =rgb）＃colors
plot（tt，col = rgb.palette（200），main = as.expression （粘贴（Temperatura a 2m，DAY，+ 00 UTC，sep =）），axes = T）
dev.off（）
pre>

提供此输出。

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

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

解决方案

下载数据文件后，我可以直接用
光栅读取。我选择乐队221（如果我没有错），这是你
需要根据
此表格：

 光栅）
 t2mc<  -  raster（'gfs.grb'，band = 221）
 
> t2mc 
类别：RasterLayer 
带：221（315个带）
尺寸：361,720,259920（nrow，ncol，ncell）
分辨率：0.5,0.5（x， y）
范围：-0.25,359.75，-90.25,90.25（xmin，xmax，ymin，ymax）
 coord。 REF。 ：+ proj = longlat + a = 6371229 + b = 6371229 + no_defs 
数据源：/home/oscar/gfs.grb 
名称：gfs 
  pre> 
 
 
您不需要整个范围，因此您可以使用 crop 获取
所需的范围：
  e<  -  extent（-40,40,20,90）
 tt<  -  crop t2mc，e）
  
我试图显示 tt  raster with  plot  without 
 success。但是，如果您使用
不同的程度（89.5而不是90），则可以正确使用 spplot ：
  e<  - 范围（-40,40,20,89.5）
 tt<  -  crop（t2mc，e）
 
 spplot tt）
  
现在我们必须添加管理界限：
 库（maps）
库（mapdata）
库（maptools）
 
 ext < -  as.vector e）
 boundary<  -  map（'worldHires'，
 xlim = ext [1：2]，ylim = ext [3：4]，
 plot = FALSE）
边界<  -  map2SpatialLines（边界，
 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'，boundary，lwd = 0.5））
  
  
 
 
 
如果您更喜欢 latticeExtra :: layer 方法，您可以通过此代码实现
a相似的结果：
  library（rasterVis）
 levelplot（tt，col.regions = rgb.palette，
 colorkey = list（height = .3））+ 
 layer（sp.lines（boundary，lwd = 0.5））
  
 
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.

解决方案
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))


                        
这篇关于R裁剪光栅数据和设置轴限制的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！