2013-06-20 8 views
7

別のスレッドであなたの助けを借りて、私はいくつかのグローバルマップをプロットすることができました。最初に、気象GRIB2のデータをNetcdfに変換し、グローバルマップをプロットします。R作画ラスタデータと軸の制限を設定する

今私はマップの部分領域をプロットしたいと思います。作物コマンドを試して、グローバルな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() 

をプロットするために使用しているスクリプトです。

cropped plot

それは簡単なものである必要がありますが、私は、単純なRのユーザーです。前もって感謝します。

EDIT:xlim = c(-40,40)、ylim = c(20,90)を追加すると新しい出力が表示されます。それは問題を解決しないようだ。しかし、出力pngファイルのx、yサイズで遊ぶことは、マップに合うようにサイズを調整できるので有望です。確かに、それは別の解決策でなければなりません。

enter image description here

+0

を'plot(....、xlim = c(-10,30)、ylim = c(30、80))'そしてちょっと素敵なプロット+1 + –

+0

googleVisパッケージを見てみましょう。それが助けになるかどうかわからないが、それはかなりきちんとしている。 IntensityMap、GeoMap、Map関数が含まれています。 –

+0

こんにちは@ SimonO101それは試してみる前に私の最初の試行だった、両方を試してもわからない。今は仕事ではなく、試してみましょう。どうもありがとうございました。 – pacomet

答えて

9

データファイルをダウンロードした後、私は ラスタで直接読み取ることができます。私は this tableによると、それは何である(私は間違っていないよ場合)というバンド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) 

私は 成功せずplotttラスタを表示しようとしています。

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

とパレットを変更します:

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

spplot(tt) 

今、私たちは行政境界線を追加する必要があります。しかし、それはあなたが 異なる程度(89.5代わりの90)を使用する場合spplotで正しく動作します

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 result

あなたがlatticeExtra::layerアプローチを好む場合は、このコードで 同様の結果を達成することができます:ちょうど例えば、あなたの `plot`のコマンドに` xlim`と `ylim`を追加

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

Gracias @ oscar-perpinanあなたの提案はうまくいきます。私は、軸、タイトルなどのための追加のspplotオプションを探すつもりです...ありがとう – pacomet

+0

私はまだプロットとソリューションを探します。 – pacomet

関連する問題