2014-12-01 11 views
6

Ich habe mit den räumlichen Daten des RCP (Representative Concentration Pathway) gearbeitet. Es ist eine nette Griddatei im netCDF-Format. Wie bekomme ich eine Liste von Bausteinen, in denen jedes Element eine Variable aus einer multivariaten netCDF-Datei repräsentiert (mit der Variable I meine ich nicht lat, lon, time, depth ... etc). Das ist, was Iv'e versucht hat. Ich kann kein Beispiel für die Daten veröffentlichen, aber ich habe das folgende Skript so eingerichtet, dass es reproduzierbar ist, wenn Sie es anzeigen möchten. Offensichtlich Fragen willkommen ... Ich habe die Sprache, die mit dem Code verbunden ist, vielleicht nicht glatt ausgedrückt. Prost.Erstellen einer Liste von Rasterbausteinen aus einer multivariaten netCDF-Datei

A: Paket Anforderungen

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

B: Daten sammeln und an der Struktur netCDF Datei

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 aussehen: für eine Variable erstellen ein Raster aus Ziegeln gebaut. Levels entsprechen Jahre

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

D: Namen zu Gesichtern

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: Fassen Sie alle Rasterzellen für jedes Jahr eine Variable "emis_ene", ich will tue dies für jede Variable der netCDF-Datei, mit der ich arbeite.

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: Noch ein Meet & Greet. Prüfen Sie, wie E

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 sieht: Der Versuch, eine Liste der Ziegel zu erstellen. Eine Liste von Objekten in Teil C erstellt

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: Lassen Sie sich von dem Junk befreien Sie heruntergeladen haben könnten ... Sorry

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

Antwort

4

Ihr (E) Schritt kann cellStats vereinfacht wird unter Verwendung von .

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

sumLayers <- sapply(vars, foo) 

sumLayers ist das Ergebnis das Sie suchen, wenn ich deine Frage richtig verstanden.

Darüber hinaus können Sie das Paket zoo verwenden, weil Sie mit Zeitreihen zu tun haben.

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

xyplot(z) 

time series

+0

Hallo Oscar, Das ist genau das, was ich war auf der Suche zu tun, und wir danken Ihnen für beide Wege die Bereitstellung vorwärts. (gute Arbeit an rasterVis BTW) ... Beste Meilen – miles2know

+1

Gut. Froh, dass ich Helfen kann. Vielen Dank für Ihr Feedback zu rasterVis. –