2014-12-01 10 views
6

Ho lavorato con i dati spaziali RCP (Representative Concentration Pathway). È un bel set di dati a griglia in formato netCDF. Come posso ottenere un elenco di mattoni in cui ogni elemento rappresenta una variabile da un file netCDF multivariato (per variabile non intendo lat, lon, time, depth ... ecc.). Questo è quello che Iv'e ha provato a fare. Non riesco a pubblicare un esempio dei dati, ma ho impostato lo script qui sotto per essere riproducibile se si vuole guardarci dentro. Ovviamente le domande sono benvenute ... Potrei non aver espresso la lingua associata al codice senza problemi. Saluti.Creazione di un elenco di mattoni raster da un file netCDF multivariato

A: requisiti della confezione

library(sp) 
    library(maptools) 
    library(raster) 
    library(ncdf) 
    library(rgdal) 
    library(rasterVis) 
    library(latticeExtra) 

B: raccogliere dati e guardare la struttura dei file netCDF

td <- tempdir() 
    tf <- tempfile(pattern = "fileZ") 

    download.file("http://tntcat.iiasa.ac.at:8787/RcpDb/download/R85_NOX.zip", tf , mode = 'wb') 
    nc <- unzip(tf , exdir = td) 
    list.files(td) 

## Take a look at the netCDF file structure, beyond this I don't use the ncdf package directly 

    ncFile <- open.ncdf(nc) 
    print(ncFile) 
    vars <- names(ncFile$var)[1:12] # I'll try to use these variable names later to make a list of bricks 

C: Creare un mattone raster per una variabile. Livelli corrispondono agli anni

r85NOXene <- brick(nc, lvar = 3, varname = "emiss_ene") 
    NAvalue(r85NOXene) <- 0 
    dim(r85NOXene) # [1] 360 720 12 

D: I nomi a facce

data(wrld_simpl) # in maptools 
    worldPolys <- SpatialPolygons([email protected]) 
    cTheme <- rasterTheme(region = rev(heat.colors(20))) 
    levelplot(r85NOXene,layers = 4,zscaleLog = 10,main = "2020 NOx Emissions From Power Plants", 
      margin = FALSE, par.settings = cTheme) + layer(sp.polygons(worldPolys)) 

Global NOx Emissions

E: Riassumere tutte le celle della griglia per ogni anno "emis_ene" una variabile, voglio fai questo per ogni variabile del file netCDF con cui sto lavorando.

gVals <- getValues(r85NOXene) 
    dim(gVals) 

    r85NOXeneA <- sapply(1:12,function(x){ mat <- matrix(gVals[,x],nrow=360) 
        matfun <- sum(mat, na.rm = TRUE) # Other conversions are needed, but not for the question 
       return(matfun) 
})

F: Un altro incontrare e salutare. Scopri come E sembra

library(ggplot2) # loaded here because of masking issues with latticeExtra 
    years <- c(2000,2005,seq(2010,2100,by=10)) 
    usNOxDat <- data.frame(years=years,NOx=r85NOXeneA) 
    ggplot(data=usNOxDat,aes(x=years,y=(NOx))) + geom_line() # names to faces again 
    detach(package:ggplot2, unload=TRUE) 

G: Tentativo di creare un elenco di mattoni. Un elenco di oggetti creati in parte C

brickLst <- lapply(1:12,function(x){ tmpBrk <- brick(nc, lvar = 3, varname = vars[x]) 
             NAvalue(tmpBrk) <- 0 
             return(tmpBrk) 

    # I thought a list of bricks would be a good structure to do (E) for each netCDF variable. 
    # This doesn't break but, returns all variables in each element of the list. 
    # I want one variable in each element of the list. 
    # with brick() you can ask for one variable from a netCDF file as I did in (C) 
    # Why can't I loop through the variable names and return on variable for each list element. 
}) 

H: Sbarazzarsi della spazzatura è stato scaricato ... Siamo spiacenti

file.remove(dir(td, pattern = "^fileZ",full.names = TRUE)) 
    file.remove(dir(td, pattern = "^R85",full.names = TRUE)) 
    close(ncFile) 

risposta

4

tuo (E) passo può essere semplificata utilizzando cellStats .

foo <- function(x){ 
    b <- brick(nc, lvar = 3, varname = x) 
    NAvalue(b) <- 0 
    cellStats(b, 'sum') 
} 

sumLayers <- sapply(vars, foo) 

sumLayers è il risultato che stai cercando, se ho capito bene la tua domanda.

Inoltre, è possibile utilizzare il pacchetto zoo perché si tratta di serie temporali.

library(zoo) 
tt <- getZ(r85NOXene) 
z <- zoo(sumLayers, tt) 

xyplot(z) 

time series

+0

Ciao Oscar, questo è esattamente quello che stavo cercando di fare e vi ringrazio per la fornitura di entrambi i modi in avanti. (ottimo lavoro su rasterVis BTW) ... Le migliori miglia – miles2know

+1

Buono. Felice di aiutare. Grazie per il tuo feedback su rasterVis. –

Problemi correlati