Atividade Pré-Curso – Produtos de Oceanografia por Satélite: Conceitos, Acesso e Processamento

Olá! Bem vindos à atividade pré-curso!

Contato: Em caso de dúvidas, enviar um e-mail para:

E-mail: [email protected]

Ao finalizar a atividade pré-curso, os participantes estarão capacitados a:

• Instalar e utilizar as ferramentas básicas para iniciar a manipular produtos de oceanografia por satélites com a linguagem de programação Python localmente.
• Acessar e baixar dados abertos dos servidores FTP da NOAA.
• Realizar de operações como:
  • Leitura de arquivos NetCDF.
  • Criação de um plot básico e visualização dos valores dos pixels.
  • Leitura de metadados, adição de títulos e legenda ao plot.
  • Adição de mapas e modificação da projeção.
  • Adição de shapefiles e textos.
  • Criação de paletas de cores personalizadas.
  • Plot de regiões específicas.
  • Leitura de dados de uma coordenada específica e sua visualização nos plots.

Durante o curso online dos dias 24, 26 e 31/05 vamos aprender também:

• Download de dados utilizando scripts.
• Criação de animações.
• Operações entre arquivos.
• Séries temporais.
• Acesso e plot de diversos produtos de oceanografia por satélite.
• Diversos tipos de plot (barbelas, vetores, etc).
• Sobreposição de produtos de oceanografia à dados do satélite GOES-16.
• Comparação de dados.
• Entre outras operações!

Duração estimada da atividade Pré-curso: 3:00 h

Abaixo, a estrutura das atividades práticas. Nesta atividade pré-curso veremos algumas operações mais básicas, que serão revisadas no dia 24/05. Já no tempo restante do dia 24/05 e nos dias 25 e 31/05, veremos novos conceitos.

Pré-requisito (desejável):

• Python 3.9 (linguagem de programação utilizada no treinamento)
• Um gerenciador de pacotes (para instalar bibliotecas)
• Um gerenciador de ambientes virtuais [envs] (nosso ambiente de programação)
• Um editor de texto ou IDE (para escrever nosso código)
• Amostras de produtos de oceanografia por satélite (dados a serem manipulados)

Nesta atividade, utilizaremos as seguintes ferramentas:

Para os primeiros três itens acima (“Python 3.9”, “Gerenciador de Pacotes” e “Gerenciador de Ambientes Virtuais”), a ferramenta “Miniconda” será suficiente.

Como editor de texto ou IDE (ambiente de desenvolvimento integrado), existem muitas opções disponíveis (Gedit, Visual Studio Code, Spyder, PyCharm, Atom, Jupyter Notebooks, etc.), mas pela simplicidade, neste procedimento usaremos o “Notepad++” (ou “Notepadqq” no Linux).

Para obter amostras de produtos de oceanografia por satélite, vamos baixar dados de servidores FTP abertos da NOAA utilizando o software Filezilla.

Observação: Execução local x Execução na nuvem:

É possível executar os scripts Python que veremos tanto localmente (em seus próprios computadores) como na nuvem (na plataforma Google Colab por exemplo). Durante o treinamento, vamos exemplificar ambos! Neste procedimento pré-curso, veremos como instalar as ferramentas necessárias e executar os scripts localmente. Já durante o curso, veremos como executar os scripts na nuvem, sem a necessidade de instalação local. Assim os alunos terão uma visão geral de ambas as possibilidades.

Obs.: O procedimento mostrado a seguir foi realizado no Windows 10.

Com estas informações preliminares, mãos a obra!

Procedimento de instalação dos software necessários:

1) Faça o download e instale como administrador a ferramenta Miniconda para Python 3.9 (aproximadamente 60 MB) no seguinte link:

docs.conda.io/en/latest/miniconda.html

Nos seguintes links, você encontrará as instruções para instalação do Miniconda: Windows e Linux.

Observações (instalador Windows):

• Durante a instalação, não é necessário selecionar “Add Anaconda to my Path environment variable”.
• Selecionar “Register Anaconda as my default Python 3.9”

Observações (instalador Linux):

• Na instalação Linux, é possível escolher o diretório de instalação com o parâmetro “-p”. Escolha o diretório adequado em sua máquina:

./Miniconda3-py39_4.11.0-Linux-x86_64.sh -p /meu_diretorio/

Observações (instalador Windows / Linux):

• A instalação durará aproximadamente 5 minutos.

Agora que temos a ferramenta Miniconda instalada, vamos criar nosso ambiente virtual Python, que neste procedimento será chamado chamado “ocean” (Obs.: você pode dar qualquer nome ao ambiente virtual).

Criando um Ambiente Virtual Python e instalando bibliotecas:

2) Vamos criar um ambiente virtual Python chamado “ocean” e instalar as seguintes bibliotecas:

• matplotlib: biblioteca para criação de plots
• netcdf4: leitura / escrita de arquivos NetCDF
• pyhdf: leitura / escrita de arquivos HDF
• cartopy: criação de mapas e outras análises de dados geoespaciais
• imageio: ler e gravar uma ampla variedade de dados de imagem
• boto3: integração com serviços na nuvem da Amazon
• gdal: reprojeção / manipulação de dados geoespaciais
• scipy: manipulação de arrays entre outros
• pandas: manipulação e análise de dados
• scikit-learn: biblioteca de aprendizado de máquina de código aberto. Vamos utilizá-la apenas para o cálculo de alguns parâmetros.

Obs: Nem todas as bibliotecas acima serão utilizadas na atividade pré-curso, mas já vamos deixar todas as bibliotecas instaladas para o curso online, caso o aluno deseje executar os scripts localmente.

No Windows: Para criar o ambiente virtual, abra o “Anaconda Prompt” recentemente instalado, como administrador:

No Linux: Não é necessário abrir o Anaconda Prompt, os comandos podem ser usados diretamente no terminal

Criando o Ambiente Virtual (env) Python

Para a criação do ambiente virtual, insira o seguinte comando:

conda create --name ocean -c conda-forge matplotlib cartopy imageio boto3 netcdf4 pyhdf gdal scipy pandas scikit-learn

Onde:

ocean: é o nome do seu ambiente virtual Python. Poderia ser qualquer nome, como “oceanografia”, “curso_oceano”, etc.

matplotlib cartopy imageio boto3 netcdf4 pyhdf gdal scipy pandas scikit-learn: as bibliotecas que desejamos instalar. Existem muitas outras bibliotecas disponíveis, estas são as que utilizaremos no treinamento.

Observações sobre a instalação das bibliotecas no ambiente:

• Durante a instalação, insira “y” e Enter para continuar, quando solicitado.
• ATENÇÃO: A criação do ambiente virtual e instalação das bibliotecas levará um certo tempo. Dependendo da conexão, poderá levar um tempo considerável.

A mensagem abaixo confirma a correta criação do ambiente virtual.

Obs.: “base” (visto entre parênteses na imagem acima) é o nome do ambiente padrão do conda, e não é recomendada sua utilização para a instalação de bibliotecas devido a possíveis conflitos entre versões. É uma boa prática criar ambientes específicos para aplicações específicas (como este treinamento por exemplo), tendo assim o maior controle sobre suas instalações.

Após o procedimento, ative o ambiente virtual “ocean” recém criado com o seguinte comando:

conda activate ocean

Observação: No Linux, você pode utilizar também o comando:

source activate ocean

Ao ativá-lo, você verá o nome “ocean” entre parênteses (ao invés de “base”). Isso significa que o ambiente “ocean” foi ativado corretamente.

Informações Adicionais sobre Ambientes Virtuais:

Comandos úteis do conda

Ainda que os seguintes comandos não sejam necessários para este treinamento, são muito úteis no desenvolvimento em geral:

• Desativar um ambiente virtual:
• conda deactivate (no ambiente já ativado)
• Ver a lista de ambientes virtuais já criados:
• conda env list
• Ver a lista de pacotes (bibliotecas, etc.) instalados em um ambiente virtual:
• conda list (no ambiente já ativado)
• Deletar um ambiente virtual (exemplo com o ambiente “ocean”):
• conda remove –name ocean –all

Encontre mais informações sobre a gestão de ambientes virtuais neste link.

Estamos prontos para começar a programar com Python, no entanto, precisamos de um editor para escrever nosso código e também precisamos de produtos de oceanografia por satélite de amostra. Vamos ver como baixá-los!

Baixando um Editor de Textos:

3-) Existem diversas opções de editores ou IDE’s disponíveis para criação dos nossos scripts. Mas a título de simplicidade (e evitar que o aluno encontre problemas na configuração das opções de editores mais avançados), vamos utilizar o Notepad++. Falaremos mais sobre isso no treinamento online.

No Windows: Faça o download e instale o Notepad++ no seguinte link (apenas 4 MB):

notepad-plus-plus.org/downloads/

No Linux: Caso deseje, você pode instalar o Notepadqq com o seguinte comando:

sudo apt-get install notepadqq

Ou então usar o Gedit, ou outro editor Linux.

Observação:

• Você pode usar qualquer editor de textos ou IDE que desejar (Gedit, Visual Studio Code, Spyder, PyCharm, Atom, Jupyter Notebooks, etc.), salvando o arquivo de código usando a extensão “.py”. Apenas não exemplificaremos a configuração destas opções durante o treinamento.

Baixando Amostras de Produtos de Oceanografia por Satélite no formato NetCDF:

4-) Crie um diretório na sua máquina para os arquivos desta atividade. Neste procedimento exemplo utilizamos o diretório VLAB\Python\ (em: C:\VLAB\Python\). Este será nosso diretório de trabalho. Use o diretório que desejar.

Vamos baixar uma amostra de Temperatura da Superfície do Mar (TSM) diária no formato NetCDF, dos servidores FTP abertos da NOAA.

Para isso, vamos utilizar o software FileZilla, que pode ser baixado no seguinte link:

filezilla-project.org/download.php

FileZilla é uma solução FTP de código aberto gratuita.

Baixe e instale a versão gratuita mais recente do FileZilla. Neste procedimento, a versão 3.59.0 foi utilizada.

Vamos baixar um arquivo global diário de TSM (5 km) do seguinte servidor FTP da NOAA:

ftp.star.nesdis.noaa.gov

Nosso arquivo de interesse está localizado no seguinte diretório:

/pub/socd/mecb/crw/data/5km/v3.1_op/nc/v1.0/daily/sst/2022

Execute o FileZilla. No campo “Host”, insira ftp.star.nesdis.noaa.gov e aperte Enter.

Aparecerá uma janela informando que o servidor não suporta FTP sobre TLS. Clique “OK”.

Após a conexão bem sucedida, em “Endereço remoto”, insira o seguinte diretório:

/pub/socd/mecb/crw/data/5km/v3.1_op/nc/v1.0/daily/sst/2022

Ao confirmar com um Enter, a lista de arquivos neste diretório será mostrada na parte inferior da interface:

Agora vamos escolher onde o FileZilla deverá salvar os arquivos baixados. Em “Endereço local”, escolha seu diretório de trabalho (neste procedimento exemplo, estamos utilizando C:\VLAB\Python. Selecione o seu diretório de trabalho.

Neste procedimento, vamos baixar o arquivo “coraltemp_v3.1_20220101.nc”, que é um arquivo diário de TSM global (5km) do dia 01/01/2022.

Os arquivos diários de TSM global da NOAA possuem a seguinte nomenclatura:

coraltemp_v3.1_YYYYMMDD.nc

Onde: YYYY: ano , MM: mês e DD: dia

Para baixar o NetCDF desta data, clique com o botão direito sobre o arquivo e selecione “Baixar”:

O console do FileZilla irá mostrar o andamento do download (pasta de destino, porcentagem, entre outras informações):

Ao finalizar o download você verá uma mensagem no canto inferior direito indicando que o download foi completado:

E o arquivo baixado estará localizado na pasta escolhida como destino:

E assim baixamos nosso dado de amostra. Como vamos ver no curso, existem muitos links abertos de FTP da NOAA, com uma infinidade de dados disponíveis, que podem ser baixados de modo semelhante! No curso vamos ver também como fazer este mesmo procedimento, porém utilizando scripts, o que é muito mais simples e rápido.

Neste ponto já temos nosso ambiente de programação e nosso editor de texto instalados e uma amostra de dado disponível. Estamos prontos para iniciar! Vamos criar nosso primeiro script.

SCRIPT 1: PLOT BÁSICO

IMPORTANTE: Nesta atividade pré-curso não vamos explicar com muitos detalhes o que faz cada instrução. Os conceitos serão revisados e explicados com mais detalhes no curso online. O mais importante nessa etapa é conseguir executar os scripts sem maiores problemas.

Abra seu editor de texto ou IDE (neste procedimento, o “Notepad ++”), insira o código abaixo e salve-o como “script_01.py” no seu diretório de trabalho (C:\VLAB\Python\ , no nosso caso). Este é o mesmo diretório onde você baixou o arquivo NetCDF no passo anterior.

Importante: Na linha 12, insira o nome do arquivo NetCDF que você baixou no passo anterior. Caso tenha sido o arquivo do dia 01/01/2022, será o mesmo nome do exemplo a seguir.

A imagem abaixo mostra o script aberto no editor Notepad++:

Obs: Na imagem exemplo acima o Notepad++ foi configurado para um estilo de cores escuras, chamado “Obsidian”, em “Configurações” > “Configurador de estilos” > “Selecionar tema:” > “Obsidian”. Isto é completamente opcional. Voltemos ao script.

Basicamente, na linha 10 especificamos o arquivo NetCDF que desejamos manipular, na linha 13, lemos o conteúdo do conjunto de dados chamado “analysed_sst” (os dados de TSM), na linha 16 configuramos o tamanho do nosso plot (10×5 polegadas neste exemplo), na linha 19 configuramos a aparência do nosso plot (temperaturas mínimas, máximas, origem do dado e paleta de cores) e na linha 22 salvamos o plot no formato PNG. Já na linha 25, mostramos o plot numa janela interativa.

O mais importante para entendermos nesse primeiro script é a seguinte linha:

# Extract the Sea Surface Temperature
data = file.variables['analysed_sst'][0,:,:]

Nela, estamos indicando que desejamos ler o “dataset” de TSM disponível nesse NetCDF, chamado “analysed_sst”. Esse dado é uma matriz de 3 dimensões, sendo: data, linhas e colunas. Neste arquivo em específico, temos apenas uma data disponível, portanto nossa matriz 3d tem apenas uma camada. Por isso indicamos o índice “0” no primeiro campo. Já nos outros dois, utilizamos dois pontos “:”, indicando que desejamos ler toda a matriz de linhas e colunas, ou seja, todo o dado global. No Script 6 deste procedimento vamos verificar como ler apenas uma região específica.

De momento, vamos apenas executar o script. No “Anaconda Prompt”, com o ambiente virtual “ocean” ativado…

conda activate open

…acesse o diretório C:\VLAB\Python\ (ou o seu diretório de trabalho):

cd C:\VLAB\Python

Execute o script “script_01.py” com o seguinte comando:

python script_01.py

A imagem abaixo mostra o procedimento acima:

Ao executar o script_01.py, veremos o plot resultante em uma janela interativa:

• Visualizando os valores dos pixels:

Ao mover o ponteiro do mouse sobre o plot, no canto superior direito você pode visualizar a linha e coluna correspondente, além dos valores de TSM (°C) de um pixel específico.

Vamos adicionar algumas informações ao plot no nosso próximo script!

SCRIPT 2: LENDO METADADOS E ADICIONANDO LEGENDA, TÍTULO E DATA

Vamos criar nosso segundo script e plotar uma imagem com mais informações sobre o dado. Crie e salve um arquivo chamado “script_02.py”.

Neste exemplo vamos ler a data do arquivo diretamente de seus metadados. Porém, como vamos ver no curso, o arquivo possui a seguinte informação relacionada à data: “segundos desde 01/01/1981 00:00 UTC”. Precisaremos somar esta quantidade de segundos à esta data inicial para calcular a data de aquisição. O Python facilita a manipulação desse tipo de informação (datas e horários) através da biblioteca “datetime”, que vamos importar neste segundo script na seguinte linha:

from datetime import datetime, timedelta  # Library to convert julian day to dd-mm-yyyy

O fragmento de código abaixo mostra o comando para adicionar uma legenda, os comandos para calcular a data do arquivo e os comandos para adicionar um título ao plot (os comandos serão explicados durante o curso):

# Add a colorbar
plt.colorbar(label='Sea Surface Temperature (°C)', extend='both', orientation='horizontal', pad=0.05, fraction=0.05)

# Getting the file time and date
add_seconds = int(file.variables['time'][0])
date = datetime(1981,1,1,0) + timedelta(seconds=add_seconds)
date_formatted = date.strftime('%Y-%m-%d')
	
# Add a title
plt.title(f'NOAA Coral Reef Watch Daily 5 km SST - {date_formatted}', fontweight='bold', fontsize=10, loc='left')
plt.title('INPE / CGCT / DISSM', fontsize=10, loc='right')

Abaixo, o script completo até o momento:

Execute o script “script_02.py” com o seguinte comando:

python script_02.py

Ao executar o script, aparecerá uma nova janela com o nosso plot, já com os títulos e legendas adicionados:

Até o momento nosso dado está sendo plotado mostrando as informações em linhas e colunas (3600 linhas x 7200 colunas), isso porque não georreferenciamos nossa imagem. No terceiro exemplo, vamos utilizar a biblioteca Cartopy para adicionar mapas e outras informações ao nosso plot.

SCRIPT 3: ADICIONANDO MAPAS COM A CARTOPY

Crie e salve um arquivo chamado “script_03.py”.

Vamos importar a bilbioteca cartopy (para adicionar mapas entre outras características) e a biblioteca numpy (para facilitar algumas operações) com as seguintes linhas:

import cartopy, cartopy.crs as ccrs       # Plot maps
import cartopy.feature as cfeature        # Common drawing and filtering operations
import numpy as np                        # Import the Numpy package

Como vamos ver, precisamos informar à biblioteca Cartopy qual a extensão do nosso dado (latitudes e longitudes mínimas e máximas). Para isso vamos ler pela primeira vez as matrizes de latitudes e longitudes do nosso dado:

# Reading the lats and lons 
lats = file.variables['lat'][:]
lons = file.variables['lon'][:]

O primeiro passo para usar a Cartopy é informar a projeção desejada (no exemplo abaixo, “PlateCaree”, ou cilíndrica equidistante) e informar qual são os limites geográficos do nosso plot:

# Use the Cilindrical Equidistant projection in cartopy
ax = plt.axes(projection=ccrs.PlateCarree())

# Define the image extent
img_extent = [lons.min(), lons.max(), lats.min(), lats.max()]

Após isso, podemos adicionar as linhas de costa, países, linhas de referência de latitudes e longitudes, além de uma máscara de continentes (pois não há dados nessas regiões):

# Add coastlines, borders and gridlines
ax.coastlines(resolution='50m', color='black', linewidth=0.8)
ax.add_feature(cartopy.feature.BORDERS, edgecolor='black', linewidth=0.5)
gl = ax.gridlines(crs=ccrs.PlateCarree(), color='white', alpha=1.0, linestyle='--', linewidth=0.25, xlocs=np.arange(-180, 180, 30), ylocs=np.arange(-90, 90, 10), draw_labels=True)
gl.top_labels = False
gl.right_labels = False

# Add a land mask
ax.add_feature(cfeature.LAND)

Abaixo, o script completo até o momento:

Execute o script “script_03.py” com o seguinte comando:

python script_03.py

Ao executar o script, aparecerá uma nova janela com o nosso plot, já com os mapas adicionados:

Veja que agora, ao mover o ponteiro do mouse sobre o plot, além das temperaturas, vemos as coordenadas de cada ponto.

Uma das características da biblioteca Cartopy é a possibilidade de mudar a projeção do plot, o que pode ser útil em alguns casos. Muitos produtos encontrados nas páginas da NOAA são mostrados na projeção Mercator por exemplo. Vamos fazer alguns testes no nosso quarto script:

SCRIPT 4: MODIFICANDO A PROJEÇÃO DO PLOT

Crie e salve um arquivo chamado “script_04.py”.

No seguinte link, podemos ver as projeções disponíveis da Cartopy:

https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html 

Nesta atividade vamos testar três projeções: Mercator, Robinson e Mollweide.

As únicas modificações em relação ao script anterior, é a declaração da projeção utilizada e a declaração da região, conforme o exemplo abaixo:

# Use the Mercator projection in cartopy
ax = plt.axes(projection=ccrs.Mollweide(central_longitude=0.0))

# Method used for "global" plots
ax.set_global()

Na segunda instrução estamos dizendo que este é um dado global, o que simplifica um pouco o nosso script.

Também vamos adicionar uma imagem de fundo (uma das opções padrão da Cartopy), ao invés da máscara de continentes:

# Add a background map
ax.stock_img()

Abaixo, o script completo até o momento:

Execute o script “script_04.py” com o seguinte comando:

python script_04.py

Ao executar o script, aparecerá uma nova janela com o nosso plot, na projeção configurada (Mercator neste caso):

Agora vamos testar outras projeções. Por ser um dado global as modificações são bem simples. Mude a linha 23 de “Mercator” para “Robinson”:

# Use the Robinson projection in cartopy
ax = plt.axes(projection=ccrs.Robinson(central_longitude=0.0))

Mude também o subtítulo do plot:

plt.title('Robinson Projection', fontsize=10, loc='right')

E mude o tamanho do plot para 13 x 7 polegadas:

# Choose the plot size (width x height, in inches)
plt.figure(figsize=(13,7))

Execute o script 4 novamente:

python script_04.py

Você verá o mesmo plot, agora na projeção Robinson:

Modifique agora a projeção para “Mollweide”:

# Use the Mollweide projection in cartopy
ax = plt.axes(projection=ccrs.Mollweide(central_longitude=0.0))

E execute o script novamente. Você verá o plot na projeção Mollweide:

Com estas simples modificações, podemos alterar a projeção do nosso plot. Nossos plots já estão bem interessantes, porém alguns detalhes podem ser melhorados. No nosso próximo script, vamos ver como criar paletas de cores personalizadas.

SCRIPT 5: PALETAS DE CORES PERSONALIZADAS

Crie e salve um arquivo chamado “script_05.py”.

Até o momento, estamos utilizando a paleta de cores “jet”, umas das opções padrão da matplotlib, que podem ser encontradas no seguinte link:

matplotlib.org/stable/tutorials/colors/colormaps.html

Apesar de variadas, nem sempre atendem as necessidades do plot, sendo necessário por vezes criarmos paletas de cores personalizadas. É possível criar paletas de cores com qualquer número de cores, com qualquer cor ou combinação de cores. Vamos ver dois exemplo, um com uma escala contínua e outra discreta.

Vamos reproduzir a paleta de cores de TSM encontrada na página Worldview da NASA, que pode ser visualizada neste link. Esta é uma paleta de intervalos contínuos, que não consta nas paletas padrão da matplotlib.

Como esta paleta não está disponível na lista de colormaps padrão da matplotlib, para reproduzi-la, precisamos especificar via código a representação de cada intervalo de cores, porém em hexadecimal. Como podemos fazer isso?

Não é difícil. Vamos acessar a seguinte página:

https://imagecolorpicker.com/

Esta página possibilita extrairmos os valores das cores de uma imagem de referência.

Na página, clique no botão “Use Your Image”:

Na próxima janela, selecione a aba “Website Url” e cole o endereço da página de interesse, como o link do Worldview informado mais acima. Clique “OK”.

A página tirará um printscreen automaticamente e você pode clicar em qualquer ponto da imagem da página para verificar o valor hexadecimal das cores, podendo assim utilizar estes valores nos scripts.

Agora vamos inserir os valores hexadecimais no nosso script e criar uma paleta de cores personalizada. Primeiramente, vamos importar o módulo “colors” da matplotlib, necessário para essa operação:

import matplotlib.colors                  # Matplotlib colors  

No exemplo, foram obtidas através da página “imagecolorpicker” 18 amostras de cores da escala do Worldview (clicando em intervalos espaçados da mesma). Além disso, especificamos a cor que será atribuída aos valores acima do limite superior e uma cor que será atribuída aos valores abaixo do limite inferior (sendo as mesmos valores hexadecimais do 1° valor e 18° valor da lista, neste exemplo).

O código abaixo mostra as cores já em hexadecimal, e as demais instruções necessárias para criar a paleta de valores contínuos. As cores do intervalo entre os valores hexadecimais informados, serão criados automaticamente pela função “LinearSegmentColormap” da matplotlib.

# Create a custom color scale:
# Reference color scale from NASA Wordview: https://worldview.earthdata.nasa.gov/
# HEX values got from: https://imagecolorpicker.com/:
colors = ["#2d001c", "#5b0351", "#780777", "#480a5e", "#1e1552", 
          "#1f337d", "#214c9f", "#2776c6", "#2fa5f1", "#1bad1d", 
          "#8ad900", "#ffec00", "#ffab00", "#f46300", "#de3b00", 
          "#ab1900", "#6b0200", '#3c0000']
cmap = matplotlib.colors.LinearSegmentedColormap.from_list("", colors)
cmap.set_over('#3c0000')
cmap.set_under('#2d001c')
vmin = -2.0
vmax = 35.0

Também, no nosso script precisamos informar na instrução “imshow” a paleta de cores (que chamamos “cmap”) criada acima:

# Plot the image
img = ax.imshow(data, vmin=vmin, vmax=vmax, origin='lower', extent=img_extent, cmap=cmap) 

Abaixo, o script completo até o momento:

Execute o script “script_05.py” com o seguinte comando:

python script_05.py

Ao executar o script, aparecerá uma nova janela com o nosso plot, com a nova paleta de cores criada via código:

Ao compararmos com o NASA Worldview para a mesma data (01/01/2022 neste exemplo), obtivemos o mesmo resultado:

Agora vamos criar uma paleta de cores de intervalos discretos, utilizando como referência a seguinte paleta da NOAA, que possui 17 cores (com exceção do branco):

Vamos colocar o seguinte link na página “imagecolorpicker.com” para que a página faça um printscreen automático:

psl.noaa.gov/map/images/sst/sst.daily.gif

E verificar os 17 valores hexadecimal desta escala de cores:

# Create a custom color scale:
# Reference color scale from NASA Wordview: https://psl.noaa.gov/map/clim/sst.shtml
# HEX values got from: https://imagecolorpicker.com/:
colors = ["#0000a1", "#0000fe", "#0034ff", "#0356fc", "#0199fc", 
          "#01daff", "#238e25", "#00a001", "#4bcc4e", "#62fe64", 
          "#f5fe81", "#ffff02", "#ffb200", "#fd5106", "#ff2600", 
          "#fb6a94", "#b52c64"]
cmap = matplotlib.colors.ListedColormap(colors)
cmap.set_over('#b52c64')
cmap.set_under('#0000a1')
vmin = 0.0
vmax = 34.0

O que muda da escala contícua para a discreta é o seguinte comando que de “LinearSegmentedColormap”:

cmap = matplotlib.colors.LinearSegmentedColormap.from_list("", colors)

Muda para “ListedColormap”:

cmap = matplotlib.colors.ListedColormap(colors)

Além disso, vamos especificar na legenda os valores que desejamos mostrar. Para que esses valores batam com as divisões, precisamos ter o número de cores da escala + 2 valores. Neste exemplo, temos 19 valores (17 cores da escala criada + 2 valores). Também, adicionamos o parâmetro “ticks” à instrução “colorbar”:

# Define the ticks to be shown
ticks = [-2, 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34]

# Add a colorbar
plt.colorbar(img, label='Sea Surface Temperature (°C)', extend='both', orientation='horizontal', pad=0.05, fraction=0.05, ticks=ticks)

Abaixo, o script completo com a paleta com valores discretos:

Ao executar o script, vemos o mesmo plot, agora com a paleta de cores discreta:

Veja que os valores selecionados na lista “ticks” batem exatamente com os intervalos da paleta de cores criada via código.

Até o momento, estamos plotando a imagem em toda sua extensão global. No próximo exemplo, vamos verificar como plotar apenas uma região específica do dado.

SCRIPT 6: PLOTANDO UMA REGIÃO ESPECÍFICA

Crie e salve um arquivo chamado “script_06.py”.

Até o momento, estamos lendo o dado de TSM do seguinte modo:

# Extract the Sea Surface Temperature
data = file.variables['analysed_sst'][0,:,:]

Esta é uma matriz de 3 dimensões (data, linhas e colunas), sendo que para este dado em específico temos apenas uma data no arquivo, portanto sempre utilizamos o índice “0” no primeiro item entre colchetes (é uma matriz 3d de camada única). Já nos outros dois campos, utilizamos dois pontos “:”, indicando que desejamos ler a matriz completa de linhas e colunas (3600 linhas x 7200 colunas), ou seja, o dado global.

Para lermos apenas os dados de uma região específica do globo, precisamos saber exatamente quais são os índices iniciais e finais das linhas e colunas que correspondem as coordenadas desejadas. Felizmente, nosso dado NetCDF possui matrizes de latitudes e longitudes, das quais podemos extrair os índices de interesse. Precisamos de 4 índices: longitude mínima, longitude máxima, latitude mínima e latitude máxima. Com esses 4 índices, podemos ler o dados apenas para a região de interesse.

Vamos ver como fazer.

O primeiro passo é indicar no script qual a região de interesse. No exemplo abaixo estamos criando uma lista de coordenadas chamada “extent”, e declarando quatro valores (longitude mínima, latitude mínima, longitude máxima, latitude máxima). No exemplo abaixo temos os valores para a região aproximada da América do Sul.

# Select the extent [min. lon, min. lat, max. lon, max. lat]
extent = [-93.0, -60.00, -25.00, 18.00] # South America

Também, vimos em um script anterior a possibilidade de extrairmos os valores de latitudes e longitudes do nosso dado utilizando os seguintes comandos:

# Reading the lats and lons 
lats = file.variables['lat'][:]
lons = file.variables['lon'][:]

Com a região desejada e as matrizes de latitude e longitude armazenadas, podemos encontrar os índices de nosso dado, do seguinte modo:

# Latitude lower and upper index
latli = np.argmin( np.abs( lats - extent[1] ) )
latui = np.argmin( np.abs( lats - extent[3] ) )
 
# Longitude lower and upper index
lonli = np.argmin( np.abs( lons - extent[0] ) )
lonui = np.argmin( np.abs( lons - extent[2] ) )

Basicamente, estamos subtraindo as coordenadas desejadas (“extent”) das nossas matrizes de latitudes e longitudes (“lats” e “lons”). Após essa subtração, encontramos os índices desejados com a função numpy “argmin”, que retorna o índice do menor valor encontrado. Se por exemplo, declarei uma longitude mínima de -93°, ao subtrair 93 de minha matriz “lons”, a posição (ou índice) que estiver mais próximo de -93°, será o índice retornado pelo “np.argmin”, pois o valor da subtração será igual ou estará próximo de 0. Meio complicado, mas explicaremos com mais detalhes no treinamento.

Com os 4 índices calculados, vamos ler a nossa matriz de dados novamente, porém, apenas a região de interesse:

# Extract the Sea Surface Temperature
data = file.variables['analysed_sst'][ 0 , latli:latui , lonli:lonui ]

Estamos dizendo que desejamos apenas os dados que estejam entre os índices mínimos e máximos de latitude e mínimos e máximos de longitude.

Após isso, precisamos dizer para a Cartopy qual é nossa região de interesse:

# Use the Cilindrical Equidistant projection in cartopy
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([extent[0], extent[2], extent[1], extent[3]], ccrs.PlateCarree())

# Define the image extent
img_extent = [extent[0], extent[2], extent[1], extent[3]]

Abaixo, o script completo até o momento:

Veja que alteramos também o tamanho do plot e o subtítulo e mudamos a barra de cores para vertical.

Execute o script “script_06.py” com o seguinte comando:

python script_06.py

Ao executar o script, veremos apenas a região de interesse sendo plotada.

Vamos adicionar ao nosso plot alguns itens adicionais, como shapefiles e textos.

SCRIPT 7: ADICIONANDO SHAPEFILES E TEXTOS

Crie e salve um arquivo chamado “script_07.py”.

Agora que já sabemos plotar uma região específica, vamos adicionar algumas decorações aos nossos plots. Neste exemplo, vamos adicionar dois shapefiles, um de estados e províncias mundiais, e outro da METAREA V (área de responsabilidade de Brasil perante a OMM, para a qual o Serviço Meteorológico Marinho emite avisos de mau tempo e previsões meteorológicas aos navegantes e toda a comunidade marítima). Além disso, vamos adicionar alguns textos indicando as subregiões da METAREA V e um texto adicional, dentro de uma caixa, no canto inferior direito do plot.

Antes de criar os scripts, vamos baixar os shapefiles e colocá-los no nosso diretório de trabalho.

O shapefile de estados e províncias mundiais pode ser baixado no seguinte link:

https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_1_states_provinces.zip

O shapefile da METAREA V pode ser baixado no seguinte link:

https://github.com/diegormsouza/Oceanography_Python_May_2022/raw/main/Metareas.zip

O conteúdo do nosso diretório de trabalho nesse momento será:

Extraia os dois arquivos “.zip” dos shapefiles no diretório de trabalho, para que possamos utilizá-los no nosso script.

O primeiro passo para adicionar shapefiles ao nosso plot é importar o módulo “shapereader” da biblioteca Cartopy:

import cartopy.io.shapereader as shpreader # Import shapefiles

Com ele, podemos importar os dois shapefiles ao nosso plot:

# Add a shapefile
shapefile = list(shpreader.Reader('ne_10m_admin_1_states_provinces.shp').geometries())
ax.add_geometries(shapefile, ccrs.PlateCarree(), edgecolor='gray',facecolor='none', linewidth=0.3)

# Add a shapefile
shapefile = list(shpreader.Reader('Metareas.shp').geometries())
ax.add_geometries(shapefile, ccrs.PlateCarree(), edgecolor='black',facecolor='none', linewidth=1.0)

Na instrução “shpreader.Reader” indicamos o nome dos arquivos “.shp” que foram extraídos dos arquivos “.zip”, no nosso caso ‘ne_10m_admin_1_states_provinces.shp’ e ‘Metareas.shp’. Também indicamos a cor no parâmetro “edgecolor” e a espessura da linha no parâmetro “linewidth”.

De acordo com a OMM, na METAREA V nos temos algumas sub-regiões (chamadas de “Forecast Areas”, conforme imagem abaixo:

O nosso shapefile da METAREA V que estamos utilizando não contém as nomenclaturas de cada sub-região, então vamos adicioná-las via código. Será um bom exemplo de como adicionar anotações ao nosso plot. Vamos utilizar o seguinte código:

# Add texts
txt1 = ax.annotate("METAREA V", xy=(-28, 9), xycoords=ccrs.PlateCarree()._as_mpl_transform(ax), color='black', size=10, clip_on=True, annotation_clip=True, horizontalalignment='center', verticalalignment='center', transform=ccrs.PlateCarree())
txt2 = ax.annotate("A", xy=(-48, -32), xycoords=ccrs.PlateCarree()._as_mpl_transform(ax), color='black', size=10, clip_on=True, annotation_clip=True, horizontalalignment='center', verticalalignment='center', transform=ccrs.PlateCarree())
txt3 = ax.annotate("B", xy=(-43, -27), xycoords=ccrs.PlateCarree()._as_mpl_transform(ax), color='black', size=10, clip_on=True, annotation_clip=True, horizontalalignment='center', verticalalignment='center', transform=ccrs.PlateCarree())
txt4 = ax.annotate("C", xy=(-46, -25), xycoords=ccrs.PlateCarree()._as_mpl_transform(ax), color='black', size=10, clip_on=True, annotation_clip=True, horizontalalignment='center', verticalalignment='center', transform=ccrs.PlateCarree())
txt5 = ax.annotate("D", xy=(-38, -22), xycoords=ccrs.PlateCarree()._as_mpl_transform(ax), color='black', size=10, clip_on=True, annotation_clip=True, horizontalalignment='center', verticalalignment='center', transform=ccrs.PlateCarree())
txt6 = ax.annotate("E", xy=(-36, -17), xycoords=ccrs.PlateCarree()._as_mpl_transform(ax), color='black', size=10, clip_on=True, annotation_clip=True, horizontalalignment='center', verticalalignment='center', transform=ccrs.PlateCarree())
txt7 = ax.annotate("F", xy=(-33, -10), xycoords=ccrs.PlateCarree()._as_mpl_transform(ax), color='black', size=10, clip_on=True, annotation_clip=True, horizontalalignment='center', verticalalignment='center', transform=ccrs.PlateCarree())
txt8 = ax.annotate("G", xy=(-36, -2), xycoords=ccrs.PlateCarree()._as_mpl_transform(ax), color='black', size=10, clip_on=True, annotation_clip=True, horizontalalignment='center', verticalalignment='center', transform=ccrs.PlateCarree())
txt9 = ax.annotate("H", xy=(-46, 2), xycoords=ccrs.PlateCarree()._as_mpl_transform(ax), color='black', size=10, clip_on=True, annotation_clip=True, horizontalalignment='center', verticalalignment='center', transform=ccrs.PlateCarree())
txt10 = ax.annotate("N", xy=(-22, -13), xycoords=ccrs.PlateCarree()._as_mpl_transform(ax), color='black', size=10, clip_on=True, annotation_clip=True, horizontalalignment='center', verticalalignment='center', transform=ccrs.PlateCarree())
txt11 = ax.annotate("S", xy=(-22, -17), xycoords=ccrs.PlateCarree()._as_mpl_transform(ax), color='black', size=10, clip_on=True, annotation_clip=True, horizontalalignment='center', verticalalignment='center', transform=ccrs.PlateCarree())

No primeiro parâmetro do comando “annotate” especificamos o texto desejado, já no parâmetro “xy” especificamos a posição do texto, neste caso, em coordenadas de longitude e latitude. Especificamos também a cor e seu tamanho, entre outros parâmetros.

E uma última anotação que adicionaremos ao nosso plot é um pequeno texto dentro de uma caixa no canto inferior direito, indicando algumas siglas do INPE:

# Add a text inside the plot
from matplotlib.offsetbox import AnchoredText
text = AnchoredText("INPE / CGCT / DISSM", loc=4, prop={'size': 7}, frameon=True)
ax.add_artist(text)

Abaixo, o script completo até o momento:

Execute o script “script_07.py” com o seguinte comando:

python script_7.py

Ao executar o script, veremos os shapefiles e anotações adicionados ao nosso plot:

Para finalizar os exemplos do nosso procedimento pré-curso, vamos ver como adicionar o valor de TSM de um determinado pixel ao nosso plot, em uma caixa de texto estilizada.

SCRIPT 8: LENDO DADOS DE UMA COORDENADA E ADICIONANDO AO PLOT

Crie e salve um arquivo chamado “script_08.py”.

No script 6 vimos como plotar uma região específica, e para isso calculamos os índices de nossa região utilizando as matrizes de latitude e longitude. Neste exemplo vamos fazer algo semelhante, porém, ao invés de quatro índices como no script 6 (min lon, min lat, max lon e max lat), vamos calcular apenas dois, que serão nossos pontos de lat e lon desejados.

Faremos no seguinte modo:

# Reading the data from a coordinate
lat_point = -30
lon_point = -30
lat_idx = np.argmin(np.abs(lats - lat_point))
lon_idx = np.argmin(np.abs(lons - lon_point))
data_point = file.variables['analysed_sst'][ 0 , lat_idx , lon_idx ].round(2)

No exemplo acima, definimos uma coordenada específica (lat: -30° e lon: -30° neste caso), calculamos os índices mais próximos a esta coordenada e extraímos da nossa matriz de dados apenas este ponto. Além disso, arredondamos o resultado em duas casas decimais (.round(2)).

Após obter os dados de interesse na nossa matriz, adicionamos a informação a uma caixa de texto, do seguinte modo:

# Adding the data as an annotation
# Add a circle
ax.plot(lon_point, lat_point, 'o', color='red', markersize=5, transform=ccrs.Geodetic(), markeredgewidth=1.0, markeredgecolor=(0, 0, 0, 1))
# Add a text
txt_offset_x = 0.8
txt_offset_y = 0.8
plt.annotate("Lat: " + str(lat_point) + "\n" + "Lon: " + str(lon_point) + "\n" + "SST: \n" + str(data_point) + ' °C', xy=(lon_point + txt_offset_x, lat_point + txt_offset_y), xycoords=ccrs.PlateCarree()._as_mpl_transform(ax), fontsize=7, fontweight='bold', color='gold', bbox=dict(boxstyle="round",fc=(0.0, 0.0, 0.0, 0.5), ec=(1., 1., 1.)), alpha = 1.0)

Primeiro adicionamos um circulo na coordenada de interesse, definimos um “offset” (em graus) para a caixa de texto (para que ela não fique sobreposta ao círculo), e finalmente criamos a caixa de texto com o conteúdo desejado.

Assim como no script 7, no primeiro parâmetro da instrução “annotate” especificamos o texto, no segundo a posição, e nos demais a cor, tamanho entre outras características. A diferença em relação ao script 7, é que utilizamos o parâmetro “bbox”, onde definimos as características da caixa de texto.

Abaixo, o script completo até o momento:

Execute o script “script_08.py” com o seguinte comando:

python script_08.py

Ao executar o script, veremos o valor da TSM para a coordenada de interesse adicionada ao plot:

Com este script nós terminamos as atividades pré-curso! Nesta atividade vimos alguns dos scripts básicos para manipulação de dados de TSM. Durante o curso, nós faremos uma revisão destes scripts e veremos novos scripts com novas operações para uma variedade de dados!

ACESSO AOS SCRIPTS DO PRÉ-CURSO E DURANTE O CURSO ONLINE

Encontre os scripts vistos na atividade pré-curso no link abaixo:

https://github.com/diegormsouza/Oceanography_Python_May_2022

No curso veremos aproximadamente 30 scripts, que serão disponibilizados no link acima.