2014-12-01 19 views
6

RCP (Temsilci Konsantrasyon Yolu) uzamsal verileriyle çalışıyorum. NetCDF formatında güzel bir gridded veri kümesi. Her elementin çok değişkenli netCDF dosyasından bir değişkeni temsil ettiği tuğlaların listesini nasıl alabilirim (değişkene göre lat, lon, zaman, derinlik ... vs). Iv'e bunu yapmaya çalıştı. Verilerin bir örneğini yayınlayamıyorum, ancak içeriye bakmak istiyorsanız aşağıdaki betiği tekrar üretilebilir hale getirdim. Açıkça sorulan sorular hoş karşılanır ... Kodla ilgili dili düzgün bir şekilde ifade etmemiş olabilirim. Şerefe.Çok değişkenli netCDF dosyasından raster tuğlaların bir listesini oluşturma

A: Paket gereksinmeleri

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

B: veri toplayın ve netCDF dosya yapısı

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 bakmak: Bir değişken için bir raster tuğla oluşturun. Seviyeleri yıllara

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

D gelmektedir: yüzlerin

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)) 
Adların

Global NOx Emissions

E: Her yıl bir değişken "emis_ene" için tüm ızgara hücreleri özetler, ben istiyorum Ben çalışıyorum netCDF dosyasının her değişken için bunu yapın.

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: Başka karşılamak ve selamlıyorum. tuğla bir liste oluşturmak için girişimi: 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 nasıl göründüğünü kontrol edin. C kısmında oluşturulan nesnelerin bir listesi

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: Maalesef ... indirdiğiniz olabilir önemsiz kurtulun

file.remove(dir(td, pattern = "^fileZ",full.names = TRUE)) 
    file.remove(dir(td, pattern = "^R85",full.names = TRUE)) 
    close(ncFile) 
Kişisel (E) aşama cellStats kullanılarak basitleştirilmiş olabilir

cevap

4

.

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

sumLayers <- sapply(vars, foo) 

sumLayers

Sorunuzu doğru anlamış eğer, aradığınız sonucudur.

Ayrıca, zaman serileriyle uğraştığınız için zoo paketini kullanabilirsiniz.

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

xyplot(z) 

time series

+0

Merhaba Oscar Bu yapmak isteyen ve ileri iki yönde sağlamak için teşekkür tam olarak budur. (rasterVis BTW'de harika bir çalışma) ... En iyi mil – miles2know

+1

İyi. Yardım etmekten memnun oldum. RasterVis hakkındaki görüşleriniz için teşekkür ederiz. –

İlgili konular