Como trabalhar com dados grib2 no Python

Colegas pesquisadores e entusiastas de SIG de código-livre,

Bem-vindos ao meu blog!

Gostaria de começar com um aviso - posso ser uma pesquisadora desta área, mas isso não significa que tudo o que faço ou escrevo aqui funcionará para você, em suas próprias configurações de desktop e versões de packages. Não tenho nenhuma responsabilidade se você perder dados ou estragar sua instalação. Eu também não autorizo nenhuma cópia do meu conteúdo.

Esta semana vamos falar sobre como abrir e trabalhar com dados raster georreferenciados do tipo GRIB, extensão .grib2.

Neste post, irei abrir o arquivo grib2, ler os dados da banda desejada e ler seu CRS e geotransformação.

Para abrir o GRIB, usarei pygrib. Para o CRS e a transformação, usarei, pelo menos, o rasterio e o GDAL. Neste post, mostrei como abrir, editar e salvar arquivos .tif com o rasterio, e neste, mostrei um pouco sobre a poderosa ferramenta GDAL Translate.

O que é GRIB?

Gridded Binary (GRIB) é um formato de dados comumente usado em meteorologia para armazenar dados históricos e de previsão do tempo. É o formato sugerido pela Organização Meteorológica Mundial (WMO).

Os dados GRIB podem ser GRIB 1 ou GRIB 2, uma das principais diferenças é a compressão usada para os rasters, por isso os arquivos GRIB2 são mais leves que os GRIB1.

Onde são encontrados dados GRIB?

Já vi dados GRIB sendo usados em dois ou três conjuntos de dados, mas o que mais me lembro é o MERGE, produto de precipitação (Sensoriamento Remoto) do INPE.

Como posso abrir arquivos GRIB fora do Python?

Normalmente eu abro os arquivos grib2 primeiro no QGIS (estou usando a versão 3.16.9), que abre sem problemas.

E como posso abrir arquivos GRIB ou GRIB2 usando Python?

Eu estou usando o pacote pygrib para abrir os arquivos. Em sua documentação, eles descrevem para duas maneiras de instalar o pygrib, usando pip ou conda. No entanto, só consegui fazer funcionar usando conda. Execute isso em seu prompt do Anaconda para instalar o pygrib:

conda install -c conda-forge pygrib

As dependências também serão instaladas pelo conda. Agora, em seu script python, importe o pacote pygrib.

import pygrib

A primeira coisa que você precisa fazer é abrir o arquivo

openpath=<path to file>
gribfile = pygrib.open(openpath)

Se você não abriu antes o arquivo com outro software como o QGIS, você pode querer saber o que há em cada banda do raster. A documentação do Pygrib apresenta uma maneira simples de mostrar os títulos de todas as bandas:

grbs.seek(0)
for grb in grbs:
    print(grb)

Isso irá gerar uma saída como esta:

1:Precipitation:kg m**-2 (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 200006021200
2:Pressure reduced to MSL:Pa (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 200006021200

Com base nessa saída, posso escolher qual banda gostaria de ler. No caso, a de Precipitação. Com base nisso, eu digito o seguinte:

selected_grb = grbs.select(name='Precipitation')[0]

Para selecionar os dados e ler as informações de latitude e longitude do arquivo:

data_actual, lats, lons = selected_grb.data()

Se eu chamar o tipo da variável data_actual

type(data_actual)

aparece isto:

numpy.ndarray

mostrando que é um ndarray. Posso usar esse array imediatamente, fazer operações com ele, etc.

Alguns arquivos grib são abertos no Python como uma matriz mascarada “Masked Array”. Uma matriz mascarada possui os dados em uma matriz salva em uma banda, e uma matriz com uma máscara binária em outra banda. A máscara geralmente é usada para mostrar quais pontos da malha do raster são válidos. Se isso acontecer, carregue o módulo para trabalhar com matrizes mascaradas:

import numpy.ma as ma

e uma etapa extra para extrair os dados é necessária:

data_extract=ma.getdata(data_actual, subok=True)

Para extrair a máscara da matriz mascarada:

mask_extract=ma.getmask(data_actual)

Como posso saber o CRS (projeção cartográfica) de um arquivo grib?

Existe uma maneira de fazer isso dentro do seu script usando pacotes Python, mas eu não consegui devido a incompatibilidades entre alguns pacotes que já estavam instalados no meu ambiente base do conda. Tive que fazer isso em um ambiente Python separado. Vou descerver as duas opções de qualquer maneira.

Também existe uma forma de saber o CRS sem usar Python, usando softwares de geoprocessamento.

Para fazer isso no QGIS (pode ser o seu primeiro ou último recurso, dependendo de quão familiarizado você está com recursos de SIG), vá para “Opção 3: abrir o raster no QGIS (ou no seu software SIG preferido)”.

Usar GDAL no Python para transformar grib2 em GeoTIFF

Para fazer isso no Python, continue lendo, listei a solução padrão-ouro e a solução alternativa.

Opção 1: se você consegue instalar o pacote OSGEO em seu ambiente padrão

Se você conseguir instalar o pacote OSGEO Python no mesmo ambiente em que está executando o restante do seu script, poderá ler a projeção dentro deste mesmo código. Isto seria o ideal.

Para instalar o pacote:

conda install -c conda-forge gdal

Depois de instalar, rode:

from osgeo import gdal
arq1 = gdal.Open (‘gribfile.grib2’)
GT_entrada = arq1.GetGeoTransform()
print(GT_entrada)

a geotransformação será impressa na tela (também é salva na variável “GT_entrada”).

Além disso, salve uma versão tif do arquivo. Para fazer isso, basta executar

save_arq = gdal.Translate('tif_file.tif’,arq1)

use rasterio para ler o tif e para printar a projeção.

with rasterio.open(‘tif_file.tif’) as tiffile:
    print(tiffile.crs)

Opção 2: caso haja incompatibilidade ao instalar o pacote GDAL

Crie um ambiente diferente apenas para instalar o pacote GDAL. Fiz isso diretamente no Anaconda Powershell Prompt.

conda create --name envgdal

Entre no ambiente

conda activate envgdal

Instale o pacote do GDAL

conda install -c conda-forge gdal

Entre no console Python simplesmente digitando

python

E rode a parte relativa ao GDAL diretamente deste local.

from osgeo import gdal
arq1 = gdal.Open (‘gribfile.grib2’)
GT_entrada = arq1.GetGeoTransform()
print(GT_entrada)
save_arq = gdal.Translate('tif_file.tif’,arq1)

Eu sugiro não rodar o rasterio junto com o GDAL neste ambiente. Execute a parte do rasterio em outro ambiente conda.

Opção 3: abrir o raster no QGIS (ou no seu software SIG preferido)

Certifique-se que tem QGIS instalado.

Abra o raster no QGIS, em seguida, abra o painel de informações e verifique o “CRS” da camada. Vou fazer isso no QGIS 3.16.9.

Abra o arquivo raster. Para fazer isso, clique em Data Source Manager.

Clique na opção Raster na barra esquerda. Selecione seu arquivo grib2, clique em Add e em Close.

No painel de camadas, clique com o botão direito na camada e clique em Properties.

Clique em Information na barra esquerda.

Deve estar disponível o sistema de coordenadas do arquivo. Caso não encontre as informações que deseja, você pode:

• Projetar o raster em um sistema de coordenadas diferente. Neste caso, o QGIS realizará a projeção adequada. Faça isso usando a ferramenta Warp (reproject).

• Salvar o arquivo grib2 como um arquivo tiff de dentro QGIS e abrir usando o rasterio. Para isso, clique com o botão direito do mouse na camada, selecione a opção Export, Save As.

Extras:

• E se eu quiser uma maneira alternativa de traduzir meu arquivo GRIB para .tif?

Outra opção seria um script bash simples usando os arquivos executáveis do GDAL. Mostrei como programar seus próprios scripts shell para o GDAL neste post.