Criar arquivo de saida a partir de dados NetCDF

Olá, eu gerei alguns resultados de campo médio de Ozonio, a partir de um arquivo .nc (netCDF) a partir deste meu resultado eu gostaria de gerar um novo arquivo no formato ascii ou .nc porém não consegui gerar esse arquivo de saída. abaixo a estrutura do Script # Carregando biblioteca para manipular arquivos netCDF library(maps) library(RNetCDF) library(fields) #========================================================================== # carregando arquivo e lendo dados - para essa versão do R usa-se o comando # var.get.nc, quando o comando original é get.var.nc (não sei o porque) dados <- open.nc('novembros20052014.nc') # lendo coordenadas espaço-temporal lat <- var.get.nc( dados, 'lat' ) lon <- var.get.nc( dados, 'lon' ) time <- var.get.nc( dados, 'time' ) #========================================================================= # lendo dados coluna total de Ozônio ColumnAmountO3 <- var.get.nc( dados, 'ColumnAmountO3' ) # dimensoes da variavel ColumnAmountO3 dims_ColumnAmountO3 <- dim(ColumnAmountO3) # tornando o arranjo 3D (ColumnAmountO3) em um 2D, organizado em ptos de grade X tempo dim(ColumnAmountO3) <- c( dims_ColumnAmountO3[1]*dims_ColumnAmountO3[2], dims_ColumnAmountO3[3] ) # calculando a média e retornado-a em 2D media_ColumnAmountO3 <- rowMeans( ColumnAmountO3) dim(media_ColumnAmountO3) <- c( dims_ColumnAmountO3[1], dims_ColumnAmountO3[2] ) #========================================================================================================== # longitude varia de 0 a 360, convertendo para -180 a 180, essa conversão é feita para plotagem sobre o mapa for (i in 1:dim(lon)) { if (lon[i]>180) { lon[i] <- lon[i]-360 } } # criando arquivo PNG que receberá o campo com o mapa #png( filename="campo_medio_O3_jan2005.png",width=600,height=800 ) # plotando mapa da America do Sul map( xlim=c(-110,-10), ylim=c(-60,10) ) map.axes() # plotando eixos title( main="Campo médio de Ozônio Novembros" ) # título do gráfico # definindo intervalo de 5 Dobson Units (DU) intervalos = seq( trunc(min(ColumnAmountO3)), trunc(max(ColumnAmountO3)), 5 ) # adicionando campo de pressao ao nivel medio do mar ao mapa # ler mais a respeito da função contour() com help(contour) contour( sort(lon), lat, media_ColumnAmountO3[ order(lon), ], add=T, levels=intervalos, lwd=2, labcex=1.2, col="black" ) # fechando arquivo PNG #dev.off() Desde já agradeço ____________________________________________________________________________ MATEUS DIAS NUNES MESTRANDO DO PROGRAMA DE PÓS-GRADUAÇÃO EM METEOROLOGIA - PPGMET UNIVERSIDADE FEDERAL DE PELOTAS - UFPEL TELEFONE: +55 (53) 81125154 ____________________________________________________________________________

Em geral eu evito usar esses pacotes de baixo nível para abrir netcdf (ncdf4, RNetCDF) porque eles normalmente exigem muito código para acessar os dados que realmente interessam. Se você usar o pacote raster, abrir e escrever um arquivo netcdf pode ser tão fácil quanto: ---------------- library(raster) b <- brick("caminho-do-arquivo.nc") # isso carrega todas as fatias de tempo do arquivo netcdf b1 <- b*10 + 2 # só um exemplo básico: multiplica os dados por 10 e depois soma dois writeRaster(b1, "caminho-do-novo-arquivo.nc", type="CDF", overwrite=T) # salva o novo arquivo ---------------- Esses seriam os passos básicos. Não tenho certeza do que o seu script faz porque eu não tenho acesso ao arquivo que você está usando Mas aparentemente você tem um arquivo com todos os setembros de 2005 a 2014 já separados e precisa fazer a média temporal de todos eles. Nesse caso seria simples: ----------------------- b <- brick("seu arquivo") #abre o seu arquivo b.mean <- calc(b, fun=mean) # faz uma média de todas as camadas (no caso, todos os novembros) writeRaster(b.mean, "novo.arquivo.nc") # escreve o novo arquivo # Caso queria plotar, o leveplot é uma boa opção require(rasterVis) levelplot(b.mean, margin=F) # Você pode também baixar o shapefile da américa do sul (fácil de achar na internet) e plotar o contorno no mapa: south.america <- shapefile(arquivo.shp-que-vc-baixou) levelplot(b.mean, margin=F) + layer(sp.lines(south.america, lwd=0.8, col="black")) Veja um exemplo do mapa: http://s11.postimg.org/clgtdgs9v/map_irrig_present_Page_1.jpg ----------------------- Greetings, -- Thiago V. dos Santos PhD student Land and Atmospheric Science University of Minnesota On Wednesday, February 17, 2016 3:12 PM, Mateus Dias Nunes <nunes.mateusdias@gmail.com> wrote: Olá, eu gerei alguns resultados de campo médio de Ozonio, a partir de um arquivo .nc (netCDF) a partir deste meu resultado eu gostaria de gerar um novo arquivo no formato ascii ou .nc porém não consegui gerar esse arquivo de saída. abaixo a estrutura do Script # Carregando biblioteca para manipular arquivos netCDF library(maps) library(RNetCDF) library(fields) #========================================================================== # carregando arquivo e lendo dados - para essa versão do R usa-se o comando # var.get.nc, quando o comando original é get.var.nc (não sei o porque) dados <- open.nc('novembros20052014.nc') # lendo coordenadas espaço-temporal lat <- var.get.nc( dados, 'lat' ) lon <- var.get.nc( dados, 'lon' ) time <- var.get.nc( dados, 'time' ) #========================================================================= # lendo dados coluna total de Ozônio ColumnAmountO3 <- var.get.nc( dados, 'ColumnAmountO3' ) # dimensoes da variavel ColumnAmountO3 dims_ColumnAmountO3 <- dim(ColumnAmountO3) # tornando o arranjo 3D (ColumnAmountO3) em um 2D, organizado em ptos de grade X tempo dim(ColumnAmountO3) <- c( dims_ColumnAmountO3[1]*dims_ColumnAmountO3[2], dims_ColumnAmountO3[3] ) # calculando a média e retornado-a em 2D media_ColumnAmountO3 <- rowMeans( ColumnAmountO3) dim(media_ColumnAmountO3) <- c( dims_ColumnAmountO3[1], dims_ColumnAmountO3[2] ) #========================================================================================================== # longitude varia de 0 a 360, convertendo para -180 a 180, essa conversão é feita para plotagem sobre o mapa for (i in 1:dim(lon)) { if (lon[i]>180) { lon[i] <- lon[i]-360 } } # criando arquivo PNG que receberá o campo com o mapa #png( filename="campo_medio_O3_jan2005.png",width=600,height=800 ) # plotando mapa da America do Sul map( xlim=c(-110,-10), ylim=c(-60,10) ) map.axes() # plotando eixos title( main="Campo médio de Ozônio Novembros" ) # título do gráfico # definindo intervalo de 5 Dobson Units (DU) intervalos = seq( trunc(min(ColumnAmountO3)), trunc(max(ColumnAmountO3)), 5 ) # adicionando campo de pressao ao nivel medio do mar ao mapa # ler mais a respeito da função contour() com help(contour) contour( sort(lon), lat, media_ColumnAmountO3[ order(lon), ], add=T, levels=intervalos, lwd=2, labcex=1.2, col="black" ) # fechando arquivo PNG #dev.off() Desde já agradeço ____________________________________________________________________________MATEUS DIAS NUNES MESTRANDO DO PROGRAMA DE PÓS-GRADUAÇÃO EM METEOROLOGIA - PPGMET UNIVERSIDADE FEDERAL DE PELOTAS - UFPEL TELEFONE: +55 (53) 81125154 ____________________________________________________________________________ _______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forne� c�igo m�imo reproduz�el.

Senhores, bom dia! Mateus, seria bom você disponibilizar o arquivo de dados (' novembros20052014.nc') pra podermos testar. Mas de antemão posso te dizer que uma opção é transformar a matriz media_ColumnAmountO3 em um objeto {raster} e salvar num dos formatos possíveis (ascii, geotiff, ...). Sugiro adicionar o código abaixo no teu script e testar. Não garanto que vá funcionar, pois sem os dados é um tiro no escuro. Retorne o resultado na lista... ### <code r> ... require(raster) tmp <- raster(media_ColumnAmountO3) r <- t(flip(tmp, 1)) ### para colocar na posição correta! plot(r) extent(r) <- c(range(lon), range(lat)) plot(r) writeFormats() # formatos possíveis writeRaster(r, "nomearq.tif", "GTiff") test <- raster("nomearq.tif") plot(test) ### </code> ================================================ Éder Comunello Agronomist (UEM), MSc in Environ. Sciences (UEM) DSc in Agricultural Systems Engineering (USP/Esalq) Brazilian Agricultural Research Corporation (Embrapa) Dourados, MS, Brazil |<O>| ================================================ GEO, -22.2752, -54.8182, 408m UTC-04:00 / DST: UTC-03:00

Senhores, bom dia! Fiz mais uns testes, com 3 alternativas {raster, ncdf4, RNetCDF} para escrever o arquivo NetCDF a partir de uma matriz. Acredito que seja possível melhorar o script do RNetCDF. Se alguém puder colaborar, será bem-vindo. ### <code r> ### Dados require(raster) data(volcano); volcano <- 10*volcano image(volcano) x <- 100*(1:nrow(volcano)) y <- 100*(1:ncol(volcano)) ########## {ncdf4} ########################################################### # browseURL("https://www.getdatajoy.com/learn/Read_and_Write_NetCDF_from_R") library(ncdf4) nc1_dim1 <- ncdim_def('EW', 'm', as.double(x)) nc1_dim2 <- ncdim_def('SN', 'm', as.double(y)) nc1_varz <- ncvar_def('Elevation', 'm', list(nc1_dim1, nc1_dim2), -1, longname = 'Data about a Volcano') nc1_out <- nc_create('volcano.nc', nc1_varz, force_v4 = TRUE) ncvar_put(nc1_out, nc1_varz, volcano) nc_close(nc1_out) print(nc_open("volcano.nc")) plot(raster("volcano.nc")) ########## {RNetCDF} ######################################################### # browseURL("https://journal.r-project.org/archive/2013-2/michna-woods.pdf") require(RNetCDF) nc2_out <- create.nc("volcano2.nc") dim.def.nc(nc2_out, "EW", length(x)) dim.def.nc(nc2_out, "SN", length(y)) var.def.nc(nc2_out, "EW", "NC_DOUBLE", "EW") var.def.nc(nc2_out, "SN", "NC_DOUBLE", "SN") var.def.nc(nc2_out, "Elevation", "NC_DOUBLE", c("EW", "SN")) att.put.nc(nc2_out, "Elevation", "_FillValue", "NC_DOUBLE", -999.9) var.put.nc(nc2_out, "EW", as.double(x)) var.put.nc(nc2_out, "SN", as.double(y)) var.put.nc(nc2_out, "Elevation", volcano) close.nc(nc2_out) print(nc_open("volcano2.nc")) ### {ncdf4} print.nc(open.nc("volcano2.nc")) ### {RNetCDF} var.get.nc(nc2, "EW") var.get.nc(nc2, "Elevation") plot(raster("volcano2.nc")) ########## {raster} ########################################################## require(raster) tmp <- raster(volcano) r <- t(flip(tmp, 1)) plot(r) extent(r) <- c(range(x), range(y)) plot(r) writeFormats() writeRaster(r, "volcano3.nc", "CDF") print(nc_open("volcano3.nc")) ### {ncdf4} plot(raster("volcano3.nc")) ### </code> ================================================ Éder Comunello Agronomist (UEM), MSc in Environ. Sciences (UEM) DSc in Agricultural Systems Engineering (USP/Esalq) Brazilian Agricultural Research Corporation (Embrapa) Dourados, MS, Brazil |<O>| ================================================ GEO, -22.2752, -54.8182, 408m UTC-04:00 / DST: UTC-03:00

Olá, novamente! Encontrei um erro no código anterior e que está presente também na primeira solução que sugeri. Ao usar {raster} a definição de extent() estava incorreta. Os valores x e y correspondem ao centro do pixel e para o valor de extent() ficar correto deve-se considerar a resolução da imagem. É pouco perceptível na pequena escala, mas pode causar problemas ao aumentar a escala. Segue o código corrigido... ### <code r> ### Dados require(raster) data(volcano); volcano <- 10*volcano image(volcano) x <- 100*(1:nrow(volcano)) y <- 100*(1:ncol(volcano)) ########## {ncdf4} ########################################################### # browseURL("https://www.getdatajoy.com/learn/Read_and_Write_NetCDF_from_R") library(ncdf4) nc1_dim1 <- ncdim_def('EW', 'm', as.double(x)) nc1_dim2 <- ncdim_def('SN', 'm', as.double(y)) nc1_varz <- ncvar_def('Elevation', 'm', list(nc1_dim1, nc1_dim2), -1, longname = 'Data about a Volcano') nc1_out <- nc_create('volcano.nc', nc1_varz, force_v4 = TRUE) ncvar_put(nc1_out, nc1_varz, volcano) nc_close(nc1_out) print(nc_open("volcano.nc")) plot(raster("volcano.nc")) ########## {RNetCDF} ######################################################### # browseURL("https://journal.r-project.org/archive/2013-2/michna-woods.pdf") require(RNetCDF) nc2_out <- create.nc("volcano2.nc") dim.def.nc(nc2_out, "EW", length(x)) dim.def.nc(nc2_out, "SN", length(y)) var.def.nc(nc2_out, "EW", "NC_DOUBLE", "EW") var.def.nc(nc2_out, "SN", "NC_DOUBLE", "SN") var.def.nc(nc2_out, "Elevation", "NC_DOUBLE", c("EW", "SN")) att.put.nc(nc2_out, "Elevation", "_FillValue", "NC_DOUBLE", -999.9) var.put.nc(nc2_out, "EW", as.double(x)) var.put.nc(nc2_out, "SN", as.double(y)) var.put.nc(nc2_out, "Elevation", volcano) close.nc(nc2_out) print(nc_open("volcano2.nc")) ### {ncdf4} print.nc(open.nc("volcano2.nc")) ### {RNetCDF} var.get.nc(nc2, "EW") var.get.nc(nc2, "Elevation") plot(raster("volcano2.nc")) ########## {raster} ########################################################## require(raster) tmp <- raster(volcano) r <- t(flip(tmp, 1)) plot(r) resx <- round(mean(x[-1]-x[-length(x)]), 6) # 100 resy <- round(mean(y[-1]-y[-length(y)]), 6) # 100 extent(r) <- c(range(x), range(y)) + c(-resx, +resx, -resy, +resy)/2 plot(r) writeFormats() writeRaster(r, "volcano3.nc", "CDF", over=T) print(nc_open("volcano3.nc")) ### {ncdf4} plot(raster("volcano3.nc")) r1 <- raster("volcano.nc") r2 <- raster("volcano2.nc") r3 <- raster("volcano3.nc") # extents cbind(sapply(paste0("r", 1:3), function(x) as.vector(extent(get(x))))) # r1 r2 r3 # [1,] 50 50 50 # [2,] 8750 8750 8750 # [3,] 50 50 50 # [4,] 6150 6150 6150 { windows() par(mfrow=c(2,2)) sapply(paste0("r", 1:3), function(x) plot(get(x), main=x)) par(mfrow=c(1,1)) } ### </code> ================================================ Éder Comunello Agronomist (UEM), MSc in Environ. Sciences (UEM) DSc in Agricultural Systems Engineering (USP/Esalq) Brazilian Agricultural Research Corporation (Embrapa) Dourados, MS, Brazil |<O>| ================================================ GEO, -22.2752, -54.8182, 408m UTC-04:00 / DST: UTC-03:00
participantes (3)
-
Mateus Dias Nunes
-
Thiago V. dos Santos
-
Éder Comunello