2013-06-20 14 views
7

Con il vostro aiuto in un altro thread sono riuscito a tracciare alcune mappe globali. Innanzitutto converto i dati meteorologici GRIB2 in Netcdf e quindi tracciamo le mappe globali.R ritaglio dati raster e impostazione limiti asse

Ora voglio tracciare solo una subregione della mappa. Ho provato il comando di ritaglio e ho estratto con successo la sottoregione del file globale nc. Ma quando sto disegnando non riesco a trovare il modo di controllare i limiti degli assi. Rappresenta una mappa più grande della regione dati, in modo che gli spazi bianchi appaiano su entrambi i lati.

Questo è lo script che sto usando per tracciare mappe

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

che danno questa uscita.

cropped plot

deve essere un semplice, ma io sono un utente R semplice. Grazie in anticipo.

MODIFICA: Nuova uscita quando si aggiunge xlim = c (-40,40), ylim = c (20,90) come suggerito. Sembra che non risolva il problema. Ma giocare con le dimensioni x, y del file png di output sembra promettente in quanto posso regolare la dimensione per adattarla alla mappa. Per essere sicuro che deve essere un'altra soluzione, quella giusta che non riesco a trovare.

enter image description here

+0

sufficiente aggiungere una 'xlim' e un' ylim' ai vostri comandi 'plot', per esempio 'trama (...., xlim = c (-10,30), ylim = c (30, 80))' E bella trama, a proposito, +1 –

+0

Forse date un'occhiata al pacchetto googleVis. Non sono sicuro che sarà d'aiuto, ma è abbastanza pulito. Contiene le funzioni IntensityMap, GeoMap e Map che potrebbero aiutare. –

+0

Ciao @ SimonO101 Questo è stato il mio primo tentativo prima di guardare il raccolto, non sono sicuro di aver provato entrambi. Non al lavoro ora, ci proverò. Grazie mille. – pacomet

risposta

9

Dopo aver scaricato il file di dati, posso leggere direttamente con raster. Scelgo di banda 221 che (se non sbaglio) è quello che si bisogno secondo le 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 

Non hanno bisogno tutta l'estensione in modo da utilizzare crop per ottenere la misura desiderata :

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

ho cercato di visualizzare il tt raster con plot senza successo. Tuttavia, funziona correttamente con spplot se si utilizza un diversa misura (89,5 invece di 90):

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

spplot(tt) 

Ora dobbiamo aggiungere i confini amministrativi:

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 cambiare la tavolozza:

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

Se si preferisce l'approccio latticeExtra::layer, è possibile ottenere un risultato simile con questo codice:

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

Gracias @ oscar-perpinan Il tuo suggerimento funziona bene. Cercherò altre opzioni spplot per l'asse, il titolo e così via ... Grazie – pacomet

+0

Cercherò ancora una soluzione con la trama. – pacomet

Problemi correlati