Tạo danh sách gạch raster từ tệp netCDF đa biến

Question

Tôi đã làm việc với dữ liệu không gian RCP (Concentration Conwayration Pathway). Đó là một tập dữ liệu gridded tốt đẹp trong định dạng netCDF. Làm thế nào tôi có thể nhận được một danh sách các viên gạch mà mỗi phần tử đại diện cho một biến từ một tập tin netCDF đa biến (bởi biến tôi không có nghĩa là lat, lon, thời gian, chiều sâu ... vv). Đây là những gì Iv'e đã ​​cố gắng làm. Tôi không thể đăng một ví dụ về dữ liệu, nhưng tôi đã thiết lập kịch bản dưới đây để có thể tái sản xuất nếu bạn muốn xem xét nó. Rõ ràng là các câu hỏi được chào đón ... Tôi có thể không thể hiện ngôn ngữ được liên kết với mã một cách suôn sẻ. Chúc mừng.

A: yêu cầu trọn gói

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

B: Thu thập dữ liệu và nhìn vào cấu trúc tập tin 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: Tạo một viên gạch raster cho một biến. Mức tương ứng với năm

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

D: Names để khuôn mặt

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: Tóm tắt tất cả các ô lưới cho mỗi năm một biến "emis_ene", tôi muốn làm điều này cho mỗi biến của tệp netCDF mà tôi đang làm việc.

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: Một gặp gỡ và chào hỏi. Kiểm tra như thế nào E trông

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: Cố gắng để tạo ra một danh sách các viên gạch. Một danh sách các đối tượng được tạo trong phần 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: Hãy loại bỏ rác bạn có thể đã tải ... Xin lỗi

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

Answer 1

(E) bước của bạn có thể được đơn giản hóa bằng cellStats .

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

sumLayers <- sapply(vars, foo) 

sumLayers là kết quả bạn đang tìm kiếm, nếu tôi hiểu chính xác câu hỏi của bạn. Bạn có thể sử dụng gói zoo vì bạn đang xử lý chuỗi thời gian.

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

xyplot(z)