Создание списка растровых кирпичей из многомерного файла netCDF

Question

Я работаю с пространственными данными RCP (представительный концентрационный путь). Это чистый сетчатый набор данных в формате netCDF. Как я могу получить список кирпичей, где каждый элемент представляет одну переменную из многомерного файла netCDF (по переменной я не имею в виду lat, lon, time, depth ... и т. Д.). Это то, что Ивэ пытался сделать. Я не могу опубликовать пример данных, но я настроил сценарий ниже, чтобы быть воспроизводимым, если вы хотите посмотреть на него. Очевидно, что вопросы приветствуются ... Возможно, я не мог четко выразить язык, связанный с кодом. Приветствия.

A: Требования к пакетам

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

B: Сбор данных и посмотреть на структуру файла 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: Создание растрового кирпича для одной переменной. Уровни соответствуют годам

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

D: Названия для лица

data(wrld_simpl) # in maptools 
    worldPolys <- SpatialPolygons(wrld_simpl@polygons) 
    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)) 

E: Обобщить все ячейки сетки для каждого года одна переменная "emis_ene", я хочу сделайте это для каждой переменной файла netCDF, с которой я работаю.

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: Еще одна встреча и приветствие. Проверьте, как 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: Попытка создать список кирпичей. Перечень объектов, созданных в части С

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: Избавьтесь от нежелательной вы, возможно, загрузили ... К сожалению

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

Answer 1

Ваш (E) шаг можно упростить, используя cellStats ,

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

sumLayers <- sapply(vars, foo) 

sumLayers результат вы ищете, если я правильно понял ваш вопрос.

Кроме того, вы можете использовать пакет zoo, потому что имеете дело с временными рядами.

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

xyplot(z)