[R-br] Criar arquivo de saida a partir de dados NetCDF

[Thiago V. dos Santos]

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 em 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 em 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.