2014-12-01 9 views
6

Pracowałem z danymi przestrzennymi RCP (Representative Concentration Pathway). To ładny gridowany zestaw danych w formacie netCDF. Jak mogę uzyskać listę cegieł, gdzie każdy element reprezentuje jedną zmienną z wielowymiarowego pliku netCDF (przez zmienną nie mam na myśli lat, lon, czas, głębokość ... itd.). Tak właśnie postąpił Iv'e. Nie mogę podać przykładu danych, ale skonfigurowałem poniższy skrypt, aby był możliwy do odtworzenia, jeśli chcesz się do niego zaglądać. Oczywiście pytania mile widziane ... Mogłem nie wyrazić płynnie języka związanego z kodem. Twoje zdrowie.Tworzenie listy klatek rastrowych z wielowymiarowego pliku netCDF

A: Wymagania dotyczące pakietu

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

B: Zbieranie danych i spojrzymy na strukturę plików 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: Tworzenie cegłę rastra dla jednej zmiennej. Poziomy odpowiadają lat

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

D: nazwisk do twarzy

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: Podsumowanie wszystkich komórkach sieci dla każdego roku jedna zmienna "emis_ene", chcę zrób to dla każdej zmiennej z pliku netCDF, nad którym pracuję.

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: Kolejny spotkać i pozdrowić. Sprawdź, jak wygląda 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: próba stworzenia listy cegieł. Lista obiektów utworzonych w części 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: Pozbądź się śmieci może pobrałeś ... Niestety

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

Odpowiedz

4

Państwa (E) krok można uprościć stosując cellStats .

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

sumLayers <- sapply(vars, foo) 

sumLayers jest wynikiem szukasz, jeśli dobrze zrozumiałam Twoje pytanie.

Co więcej, możesz użyć pakietu zoo, ponieważ masz do czynienia z szeregami czasowymi.

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

xyplot(z) 

time series

+0

Hi Oscar, jest to dokładnie to, czego szukałem do zrobienia i dziękuję za dostarczenie obie strony do przodu. (świetna praca na rastrzeVis BTW) ... Najlepsze mile – miles2know

+1

Dobrze. Miło, że mogłem pomóc. Dziękujemy za opinię na temat rasterVis. –