Operações básicas com arquivos NetCDF usando o R

Introdução

Saudações leitores,

É com prazer que eu apresento o meu blog. Nesse espaço vocês encontrarão dicas que podem ser úteis para estudantes ou pesquisadores que precisem lidar de alguma forma com dados espaciais em suas pesquisas.

O objetivo desse primeiro post é mostrar como ler, plotar e fazer operações algébricas simples em um arquivo climático relativamente comum, no format NetCDF, usando o R. A mesma lógica pode ser aplicada não apenas a dados atmosféricos, mas a qualquer dado em formato de ponto de grade: por exemplo imagens de satélite ou dados oceânicos.

Para ilustrar esse post eu vou usar as observações globais mensais de temperatura do ar em alta resolução espacial (0.5° x 0.5°) da Climate Research Unit (CRU). A versão atual desses dados (versão 4.03 no momento da redação desse texto) cobre os meses de janeiro de 1901 a dezembro 2018 e pode ser encontrada no link: crudata.uea.ac.uk/cru/data/hrg.

Entretanto, o arquivo da série completa é muito grande (~170 Mb de download e 3 Gb após descompactação). Para fins de simplificação, vamos usar o arquivo que cobre apenas janeiro de 2011 a dezembro de 2018. Portanto, baixe e descompacte o arquivo chamado cru_ts4.03.2011.2018.tmp.dat.nc.gz da página do CRU mencionada acima.

Antes de começarmos

Eu estou assumindo nesses artigos que você tenha um habilidades básicas de R, tais como instalação de pacotes, estruturas mais comuns de dados, acessar diretórios por linha de comando etc. Para obter ajuda sobre esses tópicos, recursos abundantes estão disponíveis na internet. Além de buscar no google, uma boa fonte de pesquisa é procurar por ou postar perguntas no stackoverflow.com. Além de instalar o R, eu recomendo usar o editor RStudio (www.rstudio.com). O RStudio é uma interface que torna o R muito mais simples de usar. Ele inclui editor de códigos e ferramentas de visualização de variáveis e figuras.

Nesse primeiro tutorial nós vamos usar o pacote raster. Pacotes, no R, são uma combinação de funções que realizam tarefas especializadas que não estão disponíveis nas funções padrão do R. Um desses pacotes, o raster (que eu vou mencionar frequentemente daqui em diante) é especializado em trabalhar com dados geográficos espacializados (chamados de “raster”).

Particularmente para abrir arquivos netCDF, o pacote raster precisa das funções do pacote ncdf4, que é responsável por acessar as bibliotecas netCDF que estão instaladas no seu computador. Portanto, para instalar o pacote ncdf4 você vai precisar instalar as bibliotecas netCDF no seu computador. Esse passo não é necessário caso você use o R no Windows ou no macOS. Porém, caso você esteja usando linux, por exemplo o Ubuntu, você pode instalar as bibliotecas netCDF com o seguinte comando (no seu sistema, e não no R):

sudo apt-get install libudunits2-dev libnetcdf-dev udunits-bin libudunits2-dev

Em seguida, instale os pacotes ncdf4 e raster:

# instala pacotes utilizados nesse artigo
install.packages("ncdf4")
install.packages("raster")

Vamos à prática

Vamos começar a nossa sessão no R carregando os pacotes que precisaremos para as tarefas de hoje:

# carrega os pacotes utilizados nesse artigo
library(ncdf4)
library(raster)

Lendo arquivos

Ler um arquivo NetCDF usando o pacote raster é bastante simples. Se o arquivo tivesse uma só camada de tempo (por exemplo um arquivo de topografia), bastaria usar o comando raster (você leu certo: o comando tem o mesmo nome do pacote).

Entretanto, como a maioria dos arquivos climáticos possui várias camadas de tempo - ou seja, consecutivas medidas da variável de estudo ao longo do tempo - nós precisamos usar a função brick. O comando a seguir cria uma variável chamada b, contendo todos os intervalos de tempo do arquivo NetCDF:

# lê o arquivo netcdf sem especificar variável
b <- brick("/Volumes/Storage/Research/data/cru_ts4.03.2011.2018.tmp.dat.nc")

Nota: Caso o arquivo NetCDF tenha mais de uma variável, você pode escolher qual variável abrir usando o argumento varname. No caso do arquivo que estamos analisando, o nome da variável é tmp. Dependendo do arquivo que você estiver analisando, pode ser precip ou algo parecida. Veja um exemplo de como especificar o nome da variável:

# lê o arquivo netcdf especificando variável
b <- brick("/Volumes/Storage/Research/data/cru_ts4.03.2011.2018.tmp.dat.nc", varname="tmp")

Note que o comando brick carrega todas as camadas do arquivo NetCDF, lendo automaticamente suas informações sobre resolução, data e detalhes da cobertura espacial (lat x lon). Vamos dar uma olhada no sumário da variável que acabamos de criar:

# sumariza as informações do arquivo
b
class       : RasterBrick 
dimensions  : 360, 720, 259200, 60  (nrow, ncol, ncell, nlayers)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : /Volumes/Storage/Research/data/cru_ts3.24.2011.2015.tmp.dat.nc 
names       : X2011.01.16, X2011.02.15, X2011.03.16, X2011.04.16, X2011.05.16, X2011.06.16, X2011.07.16, X2011.08.16, X2011.09.16, X2011.10.16, X2011.11.16, X2011.12.16, X2012.01.16, X2012.02.15, X2012.03.16, ... 
Date        : 2011-01-16, 2018-12-16 (min, max)
varname     : tmp

Algumas informações úteis são fornecidas por esse sumário. Por exemplo:

• dimensions: o número de linhas (latitudes), colunas (longitudes), total de pixels (linhas multiplicadas por colunas), e de camadas de tempo (datas) do arquivo.
• resolution: mostra a resolução (ou “tamanho de pixel”) dos dados (em graus nesse caso). Um pixel de 0.5°x0.5° significa que cada pixel representa uma área de aproximadamente 55km x 55km da superfície terrestre.
• extent: mostra a cobertura espacial (ou seja, os limites máximos e mínimos de latitude e longitude) do arquivo NetCDF.
• coord. ref.: o sistema de coordenadas em que o arquivo está projetado. Este arquivo está em sistema de coordenadas lat/lon com o datum WGS84.
• date: a cobertura temporal dos dados, que nesse caso vai de 16 de janeiro de 2011 (a primeira camada) até 16 de dezembro de 2018 (a última camada).

Plotar os dados

Além dessas informações básicas que extraímos do arquivo, nós podemos também fazer uma inspeção visual exploratória dos dados plotando a variável b:

# plota a variável
plot(b)

Note que, por padrão, todas as 16 primeiras camadas do raster são plotadas. Caso você queira, por exemplo, plotar apenas a primeira camada, use o sistema de indexação do R:

# plota a primeira camada da variável
plot(b[[1]])

E você também pode plotar um intervalo de camadas, por exemplo 13 a 24 - que no nosso arquivo representaria as datas de 16/jan/2012 de a 16/dez/2012.

# plota um intervalo de camadas da variável
plot(b[[13:24]])

Apesar da função plot do R fornecer figuras muito básicas para objetos raster, ela serve bem para uma inspeção visual inicial nos dados. Em um tutorial futuro, eu vou explorar maneiras de plotar objetos raster com uma aparência muito mais “profissional”. Caso você queira se antecipar, veja alguns exemplos de mapas produzidos com o pacote rasterVis: oscarperpinan.github.io/rastervis.

Operações algébricas

Operações algébricas simples em objetos raster no R não exigem nenhuma sintaxe especial. Podemos tratar a nossa variável raster b como uma variável normal do R. Para ilustrar essa ideia, vamos super que queiramos converter os nossos dados, que estão expressos em °C, para uma nova variável que estará expressa em °F. Um código bem simples para fazer isso seria:

# Converte temperatura em graus C para graus F
# a formula usada é:
# t_F é igual a t_C multiplicada por 9/5 e adicionada a 32
b_em_F <- (b * 9/5) + 32

E agora vamos conferir o resultado:

# plota a 1a data da temperatura em F
plot(b_em_F[[1]])

O mapa parece consistente, com temperaturas variando de aproximadamente -50°F (-45°C) a 105°F (40°C).

Opearações mais complicadas, por exemplo que exijam mais de uma variável, são um pouco mais elaboradas e precisariam usar as funções calc ou overlay do pacote raster. Eu abordarei operações desse tipo em um tutorial futuro.

Palavras finais

Esse tutorial mostrou algumas funções básicas para abrir e plotar arquivos NetCDF no R. No próximo post eu mostrarei como fazer recorte espacial nesses dados do CRU. Fique ligado! Se você achar o conteúdo desse post útil, fique à vontade para divulgá-lo entre outras pessoas que podem estar interessadas.

Caso você tenha alguma dúvida, comentário ou sugestão, não hesite em entrar em contato. Se você tiver alguma sugestão específica de tema de artigo, me dispare um e-mail também. Se eu souber como ajudar, escrevo o artigo assim que possível.

Cheers, Thiago.