VLAB: Processamento de Dados de Satélites Geoestacionários (Pré-Curso)

brazil_640

Atividade Pré-Curso – Processamento de dados de Satélites Geoestacionários com Python

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

E-mail: diego.souza@inpe.br

Ao finalizar as atividades, os participantes estarão capacitados a:

  • Utilização de ferramentas básicas para iniciar a manipular imagens de satélite com a linguagem de programação Python
  • Realização de operações como:
    • Leitura de arquivos NetCDF da série GOES-R (GOES-16 ou 17)
    • Criar um plot básico do GOES-16 e visualizar os valores dos pixels (temperaturas de brilho / refletâncias)
    • Modificar as paletas de cores, ler data e hora, adicionar um título e uma legenda ao plot
    • Adicionar mapas e shapefiles
    • Criação de composições RGB
    • Manipular dados do sensor GLM
  • Além disso, durante o curso no dia 04/03 vamos aprender:
    • Baixar dados da nuvem via código
    • Realizar operações de recorte e reprojeção
    • Manipular produtos derivados do GOES-16
    • Entre outras operações úteis

Duração estimada da atividade pré-curso: 2:00 h

Abaixo, a estrutura do curso. Nesta atividade pré-curso veremos algumas operações básicas, que serão revisadas na manhã do dia 04/03. Já na tarde do dia 04/03, veremos novos conceitos mais avançados. No final deste procedimento, encontre um pacote com todos os arquivos que serão utilizados no curso.

Pré-requisitos (desejável):

  • Assistir os vídeos curtos introdutórios sobre o GOES-16: Link
  • Assistir o webinar sobre mecanismos de recepção de dados de satélite e introdução ao processamento de dados: Link
SHOWCast_Cloud_27
SHOWCast_Cloud_28

Ferramentas necessárias para o pré-curso:

  • Python 3.8 (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 imagens GOES-R (dados a serem manipuladors)

Nesta atividade, utilizaremos as seguintes ferramentas:

Para os primeiros três itens acima (“Python 3.8”, “Gerenciador de Pacotes” y “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 (Visual Studio Code, Spyder, PyCharm, Atom, Jupyter Notebooks, Gedit, etc.), mas pela simplicidade, neste procedimento usaremos o “Notepad ++” (ou “Notepadqq” no Linux).

Para obter amostras de imagens do GOES-16, vamos baixar dados da nuvem (Amazon, neste caso). Também podemos obter imagens de outros mecanismos de recepção (estações GRB, GNC-A, servidores UNIDATA, etc.).

Procedimento de instalação:

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

https://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.8”

Observações (instalador Linux):

  • Na instalação Linux, é possível escolher o diretório de instalação com o parâmetro “-p”:
./Miniconda3-latest-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 exemplo será chamado chamado “workshop” (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 “workshop” e instalar as seguintes bibliotecas:

  • matplotlib: biblioteca para criação de plots
  • netcdf4: leitura / escrita de arquivos NetCDF
  • cartopy: criação de mapas e outras análises de dados geoespaciais
  • 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

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

Executing Anaconda Prompt

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

Criando o ambiente virtual Python

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

conda create --name workshop -c conda-forge matplotlib netcdf4 cartopy boto3 gdal scipy pandas

Onde:

workshop: é o nome do seu ambiente virtual Python. Poderia ser qualquer nome, como “satellite”, “goes”, etc.

matplotlib netcdf4 cartopy boto3 gdal scipy pandas: as bibliotecas que queremos instalar

Observações:

  • Durante a instalação, insira “y” + Enter para continuar, quando solicitado.
  • A criação do ambiente virtual e instalação das bibliotecas levará aproximadamente 10 minutos.

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

conda activate workshop

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

source activate workshop

Informações adicionais: ainda que os seguintes comandos não sejam necessários para este treinamento, são muito úteis:

Desativar um ambiente virtual Python: conda deactivate (no ambiente já ativado)
Ver a lista de ambientes virtuais criados: conda env list
Ver a lista de pacotes (bibliotecas, etc.) instalados em um ambiente virtual: conda list

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 algumas imagens do GOES-16 de amostra.

Baixando um Editor de Textos:

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

https://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”.

Baixando imagens do GOES-16 no formato NetCDF:

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

5) Acesse a seguinte página:

http://home.chpc.utah.edu/~u0553130/Brian_Blaylock/cgi-bin/goes16_download.cgi

Na interface, escolha as seguintes opções:

  • Satellite: Escolha “GOES-16/East”
  • Domain: Full Disk
  • Product: ABI L2 Cloud and Moisture Imagery
  • Date and Hour (UTC): Qualquer data e hora que desejar!

Clique em “Submit” e baixe um arquivo da Banda 13, clicando nos quadros azuis com os minutos da hora selecionada (escolha o minuto desejado).

Observação: Para os prints deste procedimento, baixamos arquivos do GOES-16 para o dia 17 de julho de 2019às 12:00 UTC.

GOES-Download.jpg

O NetCDF será baixado!

Coloque o NetCDF no diretório recém criado (por exemplo: D:\VLAB\Python\)

Estamos prontos para iniciar!

SCRIPT 1: PRIMEIRO PLOT E VISUALIZAÇÃO DOS VALORES DE PIXEL

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.

6) Abra seu editor de texto (neste procedimento, o “Notepad ++”), insira o código abaixo e salve-o como “Script_01.py” no seu diretório de trabalho (D:\VLAB\Python\ , no nosso caso). Este é o mesmo diretório onde você tem sua imagem NetCDF):

Importante: Na linha 9, insira o nome do arquivo NetCDF que você baixou no passo anterior.

# Training: Python and GOES-R Imagery: Script 1 - Basic Plot / Extracting Pixel Values
#-----------------------------------------------------------------------------------------------------------
# Required modules
from netCDF4 import Dataset      # Read / Write NetCDF4 files
import matplotlib.pyplot as plt  # Plotting library
#-----------------------------------------------------------------------------------------------------------
# Open the GOES-R image
# Download files at this link: http://home.chpc.utah.edu/~u0553130/Brian_Blaylock/cgi-bin/goes16_download.cgi
file = Dataset("OR_ABI-L2-CMIPF-M6C13_G16_s20191981200396_e20191981210116_c20191981210189.nc")

# Get the pixel values
data = file.variables['CMI'][:]
#-----------------------------------------------------------------------------------------------------------
# Choose the plot size (width x height, in inches)
plt.figure(figsize=(7,7))

# Plot the image
plt.imshow(data, vmin=193, vmax=313, cmap='Greys')
#-----------------------------------------------------------------------------------------------------------
# Save the image
plt.savefig('Image_01.png')

# Show the image
plt.show()

VOCÊ PODE BAIXAR O SCRIPT ACIMA (Script_01.py) NESTE LINK

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

No “Anaconda Prompt”, com o ambiente virtual Python “workshop” ativado, acesse o diretório D:\VLAB\Python\ (ou seu diretório de trabalho):

cd D:\VLAB\Python

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

python Script_01.py

Aparecerá uma nova janela com seu primeiro plot GOES-16 feito com Python:

Python_Quickstart_Plot1.jpg
  • Visualizando as temperaturas de brilho:

Mova o ponteiro do mouse sobre o plot e você verá os valores dos pixels da Banda 13 em Temperaturas de brilho (K) sendo mostrados na janela. Na imagem de exemplo abaixo, a temperatura do topo de nuvem neste pixel exemplo é de 227 K:

Python_Quickstart_Plot1b.jpg
  • Zoom em uma região específica

Para fazer zoom em uma determinada região, simplesmente clique no ícone da lupa na parte superior da janela e selecione a região que deseja visualizar.

Python_Quickstart_Plot1c.jpg

Para voltar a visualização completa, clique no ícone “Home”.

Além da janela de visualização, foi salva uma imagem PNG chamada ‘Image_01.png’ no seu diretório de trabalho (D:\VLAB\Python\ neste exemplo).

O que pode ser mudado rapidamente neste script?

A imagem a ser lida: Tente baixar outras imagens da Amazon (outras horas, datas e bandas) e faça um novo plot.

O tamanho da imagem gerada: Tente gerar imagens com outras dimensões. O que você mudaria para escolher o tamanho da imagem em cm ao invés de polegadas?

Os valores mínimos e máximos: Tente mudar os valores mínimos e máximos em K para modificar o contraste e saturação.

SCRIPT 2: CONVERSÃO PARA CELSIUS, LER A DATA E HORA, ADICIONAR UMA LEGENDA E UM TÍTULO

Vamos criar nosso segundo script. Crie e salve um arquivo chamado “Script_02.py”.

Para converter os valores dos pixels para graus Celsius, adicione a seguinte subtração na linha 12:

# Get the pixel values
data = file.variables['CMI'][:] - 273.15

Ao fazer isto, temos que modificar os valores mínimos e máximos da nossa paleta de cores, de 193 ~ 313 K para -80 ~ 40 ° C:

# Plot the image
plt.imshow(data, vmin=-80, vmax=40, cmap='Greys')

Para adicionar uma legenda (colorbar), adicione a seguinte linha ao script:

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

Para ler a data e hora do arquivo GOES-16, adicione as seguintes linhas ao script:

from datetime import datetime # Basic Dates and time types
# Extract date
date = (datetime.strptime(file.time_coverage_start, '%Y-%m-%dT%H:%M:%S.%fZ'))

Para adicionar um título, adicione as seguintes linhas ao script:

# Add a title
plt.title('GOES-16 Band 13 ' + date.strftime('%Y-%m-%d %H:%M') + ' UTC', fontweight='bold', fontsize=10, loc='left')
plt.title('Full Disk', fontsize=10, loc='right')

Abaixo, o script completo até o momento:

#Training: Python and GOES-R Imagery: Script 2 - Basic Operation / Colorbar / Title / Date
#-----------------------------------------------------------------------------------------------------------
# Required modules
from netCDF4 import Dataset      # Read / Write NetCDF4 files
import matplotlib.pyplot as plt  # Plotting library
from datetime import datetime    # Basic Dates and time types
#-----------------------------------------------------------------------------------------------------------
# Open the GOES-R image
# Download files at this link: http://home.chpc.utah.edu/~u0553130/Brian_Blaylock/cgi-bin/goes16_download.cgi
file = Dataset("OR_ABI-L2-CMIPF-M6C13_G16_s20191981200396_e20191981210116_c20191981210189.nc")

#Get the pixel values
data = file.variables['CMI'][:] - 273.15
#-----------------------------------------------------------------------------------------------------------
# Choose the plot size (width x height, in inches)
plt.figure(figsize=(7,7))

# Plot the image
plt.imshow(data, vmin=-80, vmax=40, cmap='jet')

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

# Extract the date
date = (datetime.strptime(file.time_coverage_start, '%Y-%m-%dT%H:%M:%S.%fZ'))

# Add a title
plt.title('GOES-16 Band 13 ' + date.strftime('%Y-%m-%d %H:%M') + ' UTC', fontweight='bold', fontsize=10, loc='left')
plt.title('Full Disk', fontsize=10, loc='right')
#-----------------------------------------------------------------------------------------------------------
# Save the image
plt.savefig('Image_02.png')

# Show the image
plt.show()

VOCÊ PODE BAIXAR O SCRIPT ACIMA (Script_02.py) NESTE LINK

Execute o “Script_02.py” com o seguinte comando:

python Script_02.py

A seguinte imagem será visualizada, com título e legenda:

  • Modificando a paleta de cores (cmap):

Agora, na linha 19, mude o cmap de ‘Greys’ para ‘jet’.

# Plot the image
plt.imshow(data, vmin=-80, vmax=40, cmap='jet')

Verifique no seguinte link as paletas de cores padrão do Matplotlib / Python (novas paletas de cores podem ser criadas, como veremos mais adiante):

https://matplotlib.org/3.1.0/tutorials/colors/colormaps.html

Além da janela de visualização, foi salva uma imagem PNG chamada ‘Image_02.png’ no seu diretório de trabalho (D:\VLAB\Python\ neste exemplo).

O que pode ser mudado rapidamente neste script?

  • A paleta de cores. Tente criar plots com direferentes colormaps do link acima.
  • O título da imagem. Tente mudar o tamanho e o texto do título.
  • A aparência da barra de cores da legenda. Tente modificar a orientação para “vertical” e o tamanho da barra de cores modificando o valor da variável “fraction”

SCRIPT 3: ADICIONANDO MAPAS COM CARTOPY

Vamos criar nosso terceiro script. Crie e salve um arquivo chamado “Script_03.py”.

  • Parâmetros do GOES-16 (longitude central, altitude e extensão):

Para utilizar o cartopy, além de importar a biblioteca:

import cartopy, cartopy.crs as ccrs  # Plot maps

Precisamos especificar os parâmetros do GOES-16 para usar a projeção geoestacionária com o cartopy (longitude central do satélite, altitude e extensão). De momento vamos inserir os parâmetros manualmente (veremos como ler os parâmetros do header do arquivo mais adiante).

# Use the Geostationary projection in cartopy
ax = plt.axes(projection=ccrs.Geostationary(central_longitude=-75.0, satellite_height=35786023.0))
img_extent = (-5434894.67527,5434894.67527,-5434894.67527,5434894.67527)

Este é o nosso script completo no momento:

# Training: Python and GOES-R Imagery: Script 3 - Adding a Map with Cartopy (Readin parameters from header)
#-----------------------------------------------------------------------------------------------------------
# Required modules
from netCDF4 import Dataset          # Read / Write NetCDF4 files
import matplotlib.pyplot as plt      # Plotting library
from datetime import datetime        # Basic Dates and time types
import cartopy, cartopy.crs as ccrs  # Plot maps
#-----------------------------------------------------------------------------------------------------------
# Open the GOES-R image
# Download files at this link: http://home.chpc.utah.edu/~u0553130/Brian_Blaylock/cgi-bin/goes16_download.cgi
file = Dataset("OR_ABI-L2-CMIPF-M6C13_G16_s20191981200396_e20191981210116_c20191981210189.nc")

# Get the pixel values
data = file.variables['CMI'][:] - 273.15
#-----------------------------------------------------------------------------------------------------------
# Choose the plot size (width x height, in inches)
plt.figure(figsize=(7,7))

# Use the Geostationary projection in cartopy
ax = plt.axes(projection=ccrs.Geostationary(central_longitude=-75.0, satellite_height=35786023.0))
img_extent = (-5434894.67527,5434894.67527,-5434894.67527,5434894.67527)

# Add coastlines, borders and gridlines
ax.coastlines(resolution='10m', color='white', linewidth=0.8)
ax.add_feature(cartopy.feature.BORDERS, edgecolor='white', linewidth=0.5)
ax.gridlines(color='white', alpha=0.5, linestyle='--', linewidth=0.5)

# Plot the image
img = ax.imshow(data, vmin=-80, vmax=40, origin='upper', extent=img_extent, cmap='Greys')

# Add a colorbar
plt.colorbar(img, label='Brightness Temperatures (°C)', extend='both', orientation='horizontal', pad=0.05, fraction=0.05)

# Extract date
date = (datetime.strptime(file.time_coverage_start, '%Y-%m-%dT%H:%M:%S.%fZ'))

# Add a title
plt.title('GOES-16 Band 13 ' + date.strftime('%Y-%m-%d %H:%M') + ' UTC', fontweight='bold', fontsize=10, loc='left')
plt.title('Full Disk', fontsize=10, loc='right')
#-----------------------------------------------------------------------------------------------------------
# Save the image
plt.savefig('Image_03.png')

# Show the image
plt.show()

VOCÊ PODE BAIXAR O SCRIPT ACIMA (Script_03.py) NESTE LINK

Execute o “Script_03.py” com o seguinte comando:

python Script_03.py

E o seguinte plot será criado:

  • Visualizando as latitudes, longitudes e temperaturas de brilho:

Mova o ponteiro do mouse sobre o plot e, além dos valores dos pixels da Banda 13 em temperaturas de brilho (° C), você verá as coordenadas dos pixels:

Python_Quickstart_Plot3b.jpg

Além da janela de visualização, foi salva uma imagem PNG chamada ‘Image_03.png’ no seu diretório de trabalho (D:\VLAB\Python\ neste exemplo).

O que pode ser mudado rapidamente neste script?

  • As cores dos mapas (nas variáveis “color” y “edgecolor”).
  • As espessuras das linhas (nas variáveis “linewidth”).

Neste link, veja uma lista de cores com o nome de cada cor, que podem ser especificados no código.

SCRIPT 4: LENDO PARÂMETROS DO “HEADER” (CABEÇALHO) DOS ARQUIVOS

Vamos criar nosso quarto script. Crie e salve um arquivo chamado “Script_04.py”.

  • Lendo parâmetros do GOES-16 do header (longitude central, altitude e extensão):

No script anterior, inserimos os parâmetros do GOES-16 manualmente:

# Use the Geostationary projection in cartopy
ax = plt.axes(projection=ccrs.Geostationary(central_longitude=-75.0, satellite_height=35786023.0))
img_extent = (-5434894.67527,5434894.67527,-5434894.67527,5434894.67527)

Podemos calcular estes parâmetros utilizando informações obtidas do cabeçalho do arquivo NetCDF do GOES-16, do seguinte modo:

#Use the Geostationary projection in cartopy
longitude_of_projection_origin = file.variables['goes_imager_projection'].longitude_of_projection_origin
perspective_point_height = file.variables['goes_imager_projection'].perspective_point_height
ax = plt.axes(projection=ccrs.Geostationary(central_longitude=longitude_of_projection_origin, satellite_height=perspective_point_height))

# Extent of data in decimais (2712 * 0.000056 * 35786023.0)
xmin = file.variables['x'][:].min() * perspective_point_height 
xmax = file.variables['x'][:].max() * perspective_point_height
ymin = file.variables['y'][:].min() * perspective_point_height 
ymax = file.variables['y'][:].max() * perspective_point_height
img_extent = (xmin, xmax, ymin, ymax)

Isto é muito útil, pois, por exemplo, independente da posição do satélite (seja o GOES-16, ou o GOES-17), o código será o mesmo. Existem muitas outras informações no header dos arquivos!

Este é o nosso script completo no momento:

# Training: Python and GOES-R Imagery: Script 4 - Adding a Map with Cartopy (Reading parameters from header)
#-----------------------------------------------------------------------------------------------------------
# Required modules
from netCDF4 import Dataset          # Read / Write NetCDF4 files
import matplotlib.pyplot as plt      # Plotting library
from datetime import datetime        # Basic Dates and time types
import cartopy, cartopy.crs as ccrs  # Plot maps
#-----------------------------------------------------------------------------------------------------------
# Open the GOES-R image
# Download files at this link: http://home.chpc.utah.edu/~u0553130/Brian_Blaylock/cgi-bin/goes16_download.cgi
file = Dataset("OR_ABI-L2-CMIPF-M6C13_G16_s20191981200396_e20191981210116_c20191981210189.nc")

# Get the pixel values
data = file.variables['CMI'][:] - 273.15
#-----------------------------------------------------------------------------------------------------------
# Choose the plot size (width x height, in inches)
plt.figure(figsize=(7,7))

# Use the Geostationary projection in cartopy
longitude_of_projection_origin = file.variables['goes_imager_projection'].longitude_of_projection_origin
perspective_point_height = file.variables['goes_imager_projection'].perspective_point_height
ax = plt.axes(projection=ccrs.Geostationary(central_longitude=longitude_of_projection_origin, satellite_height=perspective_point_height))

# Extent of data in decimais (2712 * 0.000056 * 35786023.0)
xmin = file.variables['x'][:].min() * perspective_point_height 
xmax = file.variables['x'][:].max() * perspective_point_height
ymin = file.variables['y'][:].min() * perspective_point_height 
ymax = file.variables['y'][:].max() * perspective_point_height
img_extent = (xmin, xmax, ymin, ymax)

# Add coastlines, borders and gridlines
ax.coastlines(resolution='10m', color='white', linewidth=0.8)
ax.add_feature(cartopy.feature.BORDERS, edgecolor='white', linewidth=0.5)
ax.gridlines(color='white', alpha=0.5, linestyle='--', linewidth=0.5)

# Plot the image
img = ax.imshow(data, vmin=-80, vmax=40, origin='upper', extent=img_extent, cmap='Greys')

# Add a colorbar
plt.colorbar(img, label='Brightness Temperatures (°C)', extend='both', orientation='horizontal', pad=0.05, fraction=0.05)

# Extract the date
date = (datetime.strptime(file.time_coverage_start, '%Y-%m-%dT%H:%M:%S.%fZ'))

# Add a title
plt.title('GOES-16 Band 13 ' + date.strftime('%Y-%m-%d %H:%M') + ' UTC', fontweight='bold', fontsize=10, loc='left')
plt.title('Full Disk', fontsize=10, loc='right')
#-----------------------------------------------------------------------------------------------------------
# Save the image
plt.savefig('Image_04.png')

# Show the image
plt.show()

VOCÊ PODE BAIXAR O SCRIPT ACIMA (Script_04.py) NESTE LINK

Execute o “Script_04.py” com o seguinte comando:

python Script_04.py

E o mesmo plot do script anterior será criado:

Além da janela de visualização, foi salva uma imagem PNG chamada ‘Image_04.png’ no seu diretório de trabalho (D:\VLAB\Python\ neste exemplo).

SCRIPT 5: ADICIONANDO SHAPEFILES

Vamos criar nosso quinto script. Crie e salve um arquivo chamado “Script_05.py”.

Faça o download de um shapefile com os estados brasileiros neste link. Extraia o arquivo zip no seu diretório de trabalho (neste exemplo, D:\VLAB\Python\), aonde você seus scripts e imagens do GOES-16.

Esse deve ser o conteúdo do seu diretório de trabalho neste momento:

Para adicionar um shapefile, precisamos importar o “shapereader” do Cartopy:

import cartopy.io.shapereader as shpreader # Import shapefiles

E adicionar o shapefile a visualização com os seguintes comandos:

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

Este será nosso script completo nesse momento:

# Training: Python and GOES-R Imagery: Script 5 - Reading a Shapefile
#-----------------------------------------------------------------------------------------------------------
# Required modules
from netCDF4 import Dataset                # Read / Write NetCDF4 files import matplotlib.pyplot as plt            # Plotting library
from datetime import datetime              # Basic Dates and time types
import cartopy, cartopy.crs as ccrs        # Plot maps
import cartopy.io.shapereader as shpreader # Import shapefiles
#-----------------------------------------------------------------------------------------------------------
# Open the GOES-R image
# Download files at this link: http://home.chpc.utah.edu/~u0553130/Brian_Blaylock/cgi-bin/goes16_download.cgi
file = Dataset("OR_ABI-L2-CMIPF-M6C13_G16_s20191981200396_e20191981210116_c20191981210189.nc")

# Get the pixel values, and convert to Celsius
data = file.variables['CMI'][:] - 273.15
#-----------------------------------------------------------------------------------------------------------
# Choose the plot size (width x height, in inches)
plt.figure(figsize=(7,7))

# Use the Geostationary projection in cartopy
ax = plt.axes(projection=ccrs.Geostationary(central_longitude=-75.0, satellite_height=35786023.0))
img_extent = (-5434894.67527,5434894.67527,-5434894.67527,5434894.67527)

# Add a shapefile
# https://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/malhas_municipais/municipio_2019/Brasil/BR/br_unidades_da_federacao.zip
shapefile = list(shpreader.Reader('BR_UF_2019.shp').geometries())
ax.add_geometries(shapefile, ccrs.PlateCarree(), edgecolor='gold',facecolor='none', linewidth=0.3)

# Add coastlines, borders and gridlines
ax.coastlines(resolution='10m', color='white', linewidth=0.8)
ax.add_feature(cartopy.feature.BORDERS, edgecolor='white', linewidth=0.8)
ax.gridlines(color='white', alpha=0.5, linestyle='--', linewidth=0.5)

# Plot the image
img = ax.imshow(data, vmin=-80, vmax=40, origin='upper', extent=img_extent, cmap='Greys')

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

# Extract the date
date = (datetime.strptime(file.time_coverage_start, '%Y-%m-%dT%H:%M:%S.%fZ'))

# Add a title
plt.title('GOES-16 Band 13 ' + date.strftime('%Y-%m-%d %H:%M') + ' UTC', fontweight='bold', fontsize=10, loc='left')
plt.title('Full Disk', fontsize=10, loc='right')
#-----------------------------------------------------------------------------------------------------------
# Save the image
plt.savefig('Image_05.png')

# Show the image
plt.show()

VOCÊ PODE BAIXAR O SCRIPT ACIMA (Script_05.py) NESTE LINK

Execute o “Script_05.py” com o seguinte comando:

python Script_05.py

Importante: Dependendo do tamanho do shapefile, o script pode demorar um certo tempo para ser executado.

O seguinte plot será criado:

Além da janela de visualização, foi salva uma imagem PNG chamada ‘Image_05.png’ no seu diretório de trabalho (D:\VLAB\Python\ neste exemplo).

O que pode ser mudado rapidamente neste script?

  • O shapefile a ser lido. Tente adicionar outros shapefiles de sua preferência (por exemplo: reios ou cidades de sua região).

Veja o exemplo abaixo, uma concessionária de energia elétrica adicionou um shapefile de suas linhas de transmissão. Muito interessante!

SCRIPT 6: GLM + ABI

Neste exemplo da atividade prática, vamos fazer o download de um arquivo do sensor GLM na página de download de dados da nuvem mostrada no início deste procedimento, como no exemplo abaixo (baixe um arquivo no mesmo range de horário da imagem ABI que você baixou anteriormente):

Vamos criar nosso sexto script. Crie e salve um arquivo chamado “Script_06.py”.

Para ler os eventos, grupos y flashes de um arquivo GLM, podemos utilizar o código abaixo:

#Open the GLM file
fileGLM = Dataset("OR_GLM-L2-LCFA_G16_s20191981200000_e20191981200200_c20191981200224.nc")

e_lats = fileGLM.variables['event_lat'][:]
e_lons = fileGLM.variables['event_lon'][:]
g_lats = fileGLM.variables['group_lat'][:]
g_lons = fileGLM.variables['group_lon'][:]
f_lats = fileGLM.variables['flash_lat'][:]
f_lons = fileGLM.variables['flash_lon'][:]

img2 = ax.plot(e_lons,e_lats,'.r', markersize=10, transform=ccrs.PlateCarree(), alpha=0.01)
img3 = ax.plot(g_lons,g_lats,'.y', markersize=5, transform=ccrs.PlateCarree(), alpha=0.5)
img4 = ax.plot(f_lons,f_lats,'.g', markersize=2.5, transform=ccrs.PlateCarree(), alpha=1)

Para sobrepor os dados do sensor GLM à imagem do sensorABI, podemos usar o seguinte script:

# Training: Python and GOES-R Imagery: Script 6 - ABI and GLM
#-----------------------------------------------------------------------------------------------------------
# Required modules
from netCDF4 import Dataset          # Read / Write NetCDF4 files
import matplotlib.pyplot as plt      # Plotting library
from datetime import datetime        # Basic Dates and time types
import cartopy, cartopy.crs as ccrs  # Plot maps
#-----------------------------------------------------------------------------------------------------------
# Open the GOES-R image
# Download files at this link: http://home.chpc.utah.edu/~u0553130/Brian_Blaylock/cgi-bin/goes16_download.cgi
file = Dataset("OR_ABI-L2-CMIPF-M6C13_G16_s20191981200396_e20191981210116_c20191981210189.nc")

# Get the pixel values
data = file.variables['CMI'][:] - 273.15

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

# Use the Geostationary projection in cartopy
ax = plt.axes(projection=ccrs.Geostationary(central_longitude=-75.0, satellite_height=35786023.0))
img_extent = (-5434894.67527,5434894.67527,-5434894.67527,5434894.67527)

# Add coastlines, borders and gridlines
ax.coastlines(resolution='10m', color='white', linewidth=0.8)
ax.add_feature(cartopy.feature.BORDERS, edgecolor='white', linewidth=0.5)
ax.gridlines(color='white', alpha=0.5, linestyle='--', linewidth=0.5)

# Plot the image
img = ax.imshow(data, vmin=-80, vmax=40, origin='upper', extent=img_extent, cmap='Greys')
#-----------------------------------------------------------------------------------------------------------
# Open the GLM file
# Download files at this link: http://home.chpc.utah.edu/~u0553130/Brian_Blaylock/cgi-bin/goes16_download.cgi
fileGLM = Dataset("OR_GLM-L2-LCFA_G16_s20191981200000_e20191981200200_c20191981200224.nc") 
e_lats = fileGLM.variables['event_lat'][:]
e_lons = fileGLM.variables['event_lon'][:]
g_lats = fileGLM.variables['group_lat'][:]
g_lons = fileGLM.variables['group_lon'][:]
f_lats = fileGLM.variables['flash_lat'][:]
f_lons = fileGLM.variables['flash_lon'][:] 
img2 = ax.plot(e_lons,e_lats,'.r', markersize=10, transform=ccrs.PlateCarree(), alpha=0.01)
img3 = ax.plot(g_lons,g_lats,'.y', markersize=5, transform=ccrs.PlateCarree(), alpha=0.5)
img4 = ax.plot(f_lons,f_lats,'.g', markersize=2.5, transform=ccrs.PlateCarree(), alpha=1)
#-----------------------------------------------------------------------------------------------------------
# Add a colorbar
plt.colorbar(img, label='Brightness Temperatures (°C)', extend='both', orientation='horizontal', pad=0.05, fraction=0.05)

# Extract the date
date = (datetime.strptime(file.time_coverage_start, '%Y-%m-%dT%H:%M:%S.%fZ'))

# Add a title
plt.title('GOES-16 Band 13 and GLM ' + date.strftime('%Y-%m-%d %H:%M') + ' UTC', fontweight='bold', fontsize=10, loc='left')
plt.title('Full Disk', fontsize=10, loc='right')
#-----------------------------------------------------------------------------------------------------------
# Save the image
plt.savefig('Image_06.png')

# Show the image
plt.show()

VOCÊ PODE BAIXAR O SCRIPT ACIMA (Script_06.py) NESTE LINK

Execute o “Script_06.py” com o seguinte comando:

python Script_06.py

O seguinte plot será criado:

Dê um zoom em uma região específica:

O que pode ser mudado rapidamente neste script?

  • Aparência dos pontos. Tente modificar as cores e tamanho dos pontos onde as descargas elétricas foram detectadas pelo sensor GLM.

SCRIPT 7: CRIANDO UMA COMPOSIÇÃO RGB (VERMELHO, VERDE, AZUL)

Vamos criar uma composição RGB Airmass com Python. Veja o procedimento de criação da composição Airmass encontrado no Guia Rápido do RAMMB:

ARM_Quick_Guide

Essa é a “receita” da composição Airmass:

Airmass_Recipe.jpg

Passo 1) Primeiramente, de acordo com a receita, precisamos ler quatro Bandas GOES-16:

  • Banda 08 (6.2 um)
  • Banda 10 (7.3 um)
  • Banda 12 (9.6 um)
  • Banda 13 (10.3 um)

Acesse a página de download da Amazon mostrada no início deste procedimento e baixe os NetCDF’s para a mesma data e horário para os 4 canais acima.

http://home.chpc.utah.edu/~u0553130/Brian_Blaylock/cgi-bin/goes16_download.cgi

Este deve ser o conteúdo do seu diretório de trabalho neste momento:

Vamos criar nosso sétimo script. Crie e salve um arquivo chamado “Script_07.py”.

Já sabemos como ler um arquivo do GOES-16, não é mesmo? Vamos ler 4 arquivos agora:

# Open the GOES-R image
file1 = Dataset("OR_ABI-L2-CMIPF-M6C08_G16_s20191981200396_e20191981210104_c20191981210182.nc")
# Get the pixel values, and convert to Celsius
data1 = file1.variables['CMI'][:] - 273.15

# Open the GOES-R image
file2 = Dataset("OR_ABI-L2-CMIPF-M6C10_G16_s20191981200396_e20191981210116_c20191981210188.nc")
# Get the pixel values, and convert to Celsius
data2 = file2.variables['CMI'][:] - 273.15

# Open the GOES-R image
file3 = Dataset("OR_ABI-L2-CMIPF-M6C12_G16_s20191981200396_e20191981210111_c20191981210185.nc")
# Get the pixel values, and convert to Celsius
data3 = file3.variables['CMI'][:] - 273.15

# Open the GOES-R image
file4 = Dataset("OR_ABI-L2-CMIPF-M6C13_G16_s20191981200396_e20191981210116_c20191981210189.nc")
# Get the pixel values, and convert to Celsius
data4 = file4.variables['CMI'][:] - 273.15

Passo 2) Após isso, devemos realizar as operações que vimos na receita (coluna central):

# RGB Components
R = data1 - data2
G = data3 - data4
B = data1

Passo 3) Também precisamos estabelecer os mínimos, máximos e valores de gamma de cada componente, de acordo com a receita:

# Minimuns, Maximuns and Gamma
Rmin = -26.2
Rmax = 0.6

Gmin = -43.2
Gmax = 6.7

Bmin = -29.25
Bmax = -64.65

R[R > Rmax] = Rmax
R[R < Rmin] = Rmin

G[G > Gmax] = Gmax
G[G < Gmin] = Gmin

B[B < Bmax] = Bmax (inverted!)
B[B > Bmin] = Bmin (inverted!)

gamma_R = 1
gamma_G = 1
gamma_B = 1

Passo 4) E vamos normalizar os dados entre 0 y 1:

# Normalize the data
R = ((R - Rmin) / (Rmax - Rmin)) ** (1/gamma_R)
G = ((G - Gmin) / (Gmax - Gmin)) ** (1/gamma_G)
B = ((B - Bmin) / (Bmax - Bmin)) ** (1/gamma_B)

Passo 5) Finalmente, vamos criar um “stack” dos componentes R, G y B para criar a composição final:

# Create the RGB
RGB = np.stack([R, G, B], axis=2)

# Eliminate values outside the globe
mask = (RGB == [R[0,0],G[0,0],B[0,0]]).all(axis=2)
RGB[mask] = np.nan

Passo 6) Vamos criar o plot!

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

# Use the Geostationary projection in cartopy
ax = plt.axes(projection=ccrs.Geostationary(central_longitude=-75.0, satellite_height=35786023.0))
img_extent = (-5434894.67527,5434894.67527,-5434894.67527,5434894.67527)

# Add coastlines, borders and gridlines
ax.coastlines(resolution='10m', color='white', linewidth=0.8)
ax.add_feature(cartopy.feature.BORDERS, edgecolor='white', linewidth=0.8)
ax.gridlines(color='white', alpha=0.5, linestyle='--', linewidth=0.5)

# Plot the image
img = ax.imshow(RGB, origin='upper', extent=img_extent)

# Extract date
date = (datetime.strptime(file1.time_coverage_start, '%Y-%m-%dT%H:%M:%S.%fZ'))

# Add a title
plt.title('GOES-16 Airmass RGB ' + date.strftime('%Y-%m-%d %H:%M') + ' UTC', fontweight='bold', fontsize=10, loc='left')
plt.title('Full Disk', fontsize=10, loc='right')

# Save the image
plt.savefig('Image_07.png')

# Show the image
plt.show()

Abaixo, nosso script completo:

# Training: Python and GOES-R Imagery: Script 7 - RGB Creation
#-----------------------------------------------------------------------------------------------------------
# Required modules
from netCDF4 import Dataset                # Read / Write NetCDF4 files
import matplotlib.pyplot as plt            # Plotting library
from datetime import datetime              # Basic Dates and time types
import cartopy, cartopy.crs as ccrs        # Plot maps
import numpy as np                         # Import the Numpy packag
#-----------------------------------------------------------------------------------------------------------
# Open the GOES-R image
# Download files at this link: http://home.chpc.utah.edu/~u0553130/Brian_Blaylock/cgi-bin/goes16_download.cgi
file1 = Dataset("OR_ABI-L2-CMIPF-M6C08_G16_s20191981200396_e20191981210104_c20191981210182.nc")
# Get the pixel values, and convert to Celsius
data1 = file1.variables['CMI'][:,:][::4 ,::4] - 273.15
# Open the GOES-R image
file2 = Dataset("OR_ABI-L2-CMIPF-M6C10_G16_s20191981200396_e20191981210116_c20191981210188.nc")
# Get the pixel values, and convert to Celsius
data2 = file2.variables['CMI'][:,:][::4 ,::4] - 273.15
# Open the GOES-R image
file3 = Dataset("OR_ABI-L2-CMIPF-M6C12_G16_s20191981200396_e20191981210111_c20191981210185.nc")
# Get the pixel values, and convert to Celsius
data3 = file3.variables['CMI'][:,:][::4 ,::4] - 273.15
# Open the GOES-R image
file4 = Dataset("OR_ABI-L2-CMIPF-M6C13_G16_s20191981200396_e20191981210116_c20191981210189.nc")
# Get the pixel values, and convert to Celsius
data4 = file4.variables['CMI'][:,:][::4 ,::4] - 273.15
#-----------------------------------------------------------------------------------------------------------
# RGB Quick Guide: http://rammb.cira.colostate.edu/training/visit/quick_guides/QuickGuide_GOESR_AirMassRGB_final.pdf

# RGB Components
R = data1 - data2
G = data3 - data4
B = data1

# Minimuns, Maximuns and Gamma
Rmin = -26.2
Rmax = 0.6
Gmin = -43.2
Gmax = 6.7
Bmin = -29.25
Bmax = -64.65

R[R > Rmax] = Rmax
R[R < Rmin] = Rmin
G[G > Gmax] = Gmax
G[G < Gmin] = Gmin
B[B < Bmax] = Bmax # Inverted! B[B > Bmin] = Bmin # Inverted!

gamma_R = 1
gamma_G = 1
gamma_B = 1

# Normalize the data
R = ((R - Rmin) / (Rmax - Rmin)) ** (1/gamma_R)
G = ((G - Gmin) / (Gmax - Gmin)) ** (1/gamma_G)
B = ((B - Bmin) / (Bmax - Bmin)) ** (1/gamma_B) 

# Create the RGB
RGB = np.stack([R, G, B], axis=2)

# Eliminate values outside the globe
mask = (RGB == [R[0,0],G[0,0],B[0,0]]).all(axis=2)
RGB[mask] = np.nan
#-----------------------------------------------------------------------------------------------------------
# Choose the plot size (width x height, in inches)
plt.figure(figsize=(7,7)) 

# Use the Geostationary projection in cartopy
ax = plt.axes(projection=ccrs.Geostationary(central_longitude=-75.0, satellite_height=35786023.0))
img_extent = (-5434894.67527,5434894.67527,-5434894.67527,5434894.67527)

# Add coastlines, borders and gridlines
ax.coastlines(resolution='10m', color='white', linewidth=0.8)
ax.add_feature(cartopy.feature.BORDERS, edgecolor='white', linewidth=0.8)
ax.gridlines(color='white', alpha=0.5, linestyle='--', linewidth=0.5)

# Plot the image
img = ax.imshow(RGB, origin='upper', extent=img_extent)

# Extract the date
date = (datetime.strptime(file1.time_coverage_start, '%Y-%m-%dT%H:%M:%S.%fZ'))

# Add a title
plt.title('GOES-16 Airmass RGB ' + date.strftime('%Y-%m-%d %H:%M') + ' UTC', fontweight='bold', fontsize=10, loc='left')
plt.title('Full Disk', fontsize=10, loc='right')
#-----------------------------------------------------------------------------------------------------------
# Save the image
plt.savefig('Image_07.png')

# Show the image
plt.show()

VOCÊ PODE BAIXAR O SCRIPT ACIMA (Script_07.py) NESTE LINK

Execute o “Script_07.py” com o seguinte comando:

python Script_07.py

O seguinte plot será criado:

O que pode ser mudado rapidamente neste script?

  • Adicionar shapefiles. Tente adicionar um shapefile de estados brasileiros, como fizemos no exercício 5.
  • Criação de outros RGB’s. Consulte os Guías Rápidos do RAMMB e tente criar outros RGB’s. O método de criação do RGB “Dust” é muito similar! Observação: O canal azul não está invertido na composição “Dust”.

SCRIPT 8: REALCES E SUBPLOTS

Neste último script do nosso pré-curso, vamos aprender dois conceitos interessantes: Realces de imagens e criação de subplots. Sobre os realces, vamos aprender a realçar uma imagem utilizando 4 métodos:

  • Juntando duas colormpas padrão da Matplotlib.
  • Juntando uma colormap customizada e uma colormap padrão.
  • Apenas cores customizadas.
  • Usando arquivos “CPT”

Neste script também vamos introduzir um arquivo com diversas funções úteis, chamado “utilities.py”. Vamos usá-lo bastante no curso do dia 04/03, e neste exemplo número 8, vamos utilizar a função que converte arquivos CPT para uma colormap Python.

VOCÊ PODE BAIXAR O SCRIPT “utilities.py” NESTE LINK

Baxe também o arquivo cpt “IR4AVHRR6.cpt”, que será utilizado em um dos exemplos de realce.

VOCÊ PODE BAIXAR O ARQUIVO “IR4AVHRR6.cpt” NESTE LINK

Este deve ser o conteúdo do seu diretório de trabalho neste momento:

Passo 1) Primeiramente, vamos importar as bilbiotecas necessárias. Vejam que estamos importando agora a biblioteca “cm” com algumas funções para manipulação de colormaps, além da função “loadCPT”, que está dentro do nosso arquivo “utilities.py”:

# Required modules
from netCDF4 import Dataset          # Read / Write NetCDF4 files
import matplotlib.pyplot as plt      # Plotting library
from datetime import datetime        # Basic Dates and time types
from matplotlib import cm            # Colormap handling utilities
import numpy as np                   # Scientific computing with Python
from utilities import loadCPT        # Import the CPT convert function

Passo 2) No nosso primeiro exemplo, vamos utilizar o range de -80 a 40 °C (ou seja, 120 valores). Então lemos a escala de cinza invertida (120 valores) e a convertemos em um array. Para o realce, lemos a paleta ‘jet’ invertida, com 40 valores (pois desejamos realçar entre -80 e -40 °C). Após isso, juntamos ambas as paletas de cores em um novo colormap, customizado:

# COLORMAP EXAMPLE 1
# Custom colormap joining matplotlib colormaps
# Available matplotlib colormaps: https://matplotlib.org/stable/tutorials/colors/colormaps.html
vmin1 = -80                                                       # Min. value
vmax1 = 40                                                        # Max. value
gray_cmap = cm.get_cmap('gray_r', 120)                            # Read the reversed 'gray' cmap
gray_cmap = gray_cmap(np.linspace(0, 1, 120))                     # Create the array
jet_cmap  = cm.get_cmap('jet_r', 40)                              # Read the reversed 'jet' cmap 
jet_cmap  = jet_cmap(np.linspace(0, 1, 40))                       # Create the array
gray_cmap[:40, :] = jet_cmap                                      # Join both cmaps arrays
my_cmap1 = cm.colors.ListedColormap(gray_cmap)                    # Create the custom colormap

Passo 3) No nosso segundo exemplo, ao invés de usar a colormap ‘jet’ para o realce, vamos criar uma paleta customizada, mas desta vez para realçar entre -30 e -80 °C (50 valores):

# COLORMAP EXAMPLE 2
# Creating the INPE DISSM IR colormap
# Online color picker: https://imagecolorpicker.com/
vmin2 = -80                                                       # Min. value
vmax2 = 40                                                        # Max. value
gray_cmap = cm.get_cmap('gray_r', 120)                            # Read the reversed 'gray' cmap
gray_cmap = gray_cmap(np.linspace(0, 1, 120))                     # Create the array
colors = ["#ffa0ff", "#0806ff", "#3bcfff", "#feff65", "#ff7516"]  # Custom colors
my_colors = cm.colors.ListedColormap(colors)                      # Create a custom colormap
my_colors = my_colors(np.linspace(0, 1, 50))                      # Create the array
gray_cmap[:50, :] = my_colors                                     # Join both cmaps arrays
my_cmap2 = cm.colors.ListedColormap(gray_cmap)                    # Create the custom colormap

Passo 4) No nosso terceiro exemplo, vamos utilizar uma paleta de cores apenas com cores customizadas:

# COLORMAP EXAMPLE 3
# Creating a linear colormap
# Online color picker: https://imagecolorpicker.com/
colors = ["#bc8462", "#ae656f", "#a44a79", "#962e97", "#6158c5", "#2b8ffb", "#5fcdff", "#94fff0", "#a5ff94", "#fff88c", "#ffbf52", "#ec7b27", "#b84827", "#a1333d", "#bd5478", "#cc6a99", "#d982b8"]
my_cmap3 = cm.colors.LinearSegmentedColormap.from_list("", colors)# Create a custom linear colormap
vmin3 = -80                                                       # Min. value
vmax3 = 40 

Passo 5) No nosso quarto exemplo, vamos utilizar uma paleta de cores convertida de um arquivo cpt:

# COLORMAP EXAMPLE 4
# Converts a CPT file to be used in Python
CPT archive: http://soliton.vm.bytemark.co.uk/pub/cpt-city/
cpt = loadCPT('IR4AVHRR6.cpt')                                    # Load the CPT file   
my_cmap4 = cm.colors.LinearSegmentedColormap('cpt', cpt)          # Create a custom linear colormap
vmin4 = -103.0                                                    # Min. value
vmax4 = 84.0 

Passo 6) Para criar subplots, primeiramente declaramos quantas linhas e colunas terá nosso plot completo (no exemplo abaixo, 2 linhas x 2 colunas, ou seja, 4 plots):

# Choose the plot size (width x height, in inches)
fig, axs = plt.subplots(2,2, figsize=(10,10)) # 2 rows x 2 columns 

Passo 7) Para cada subplot, especificamos em qual linha e em qual coluna desejamos adicionar o plot, na variável “axs[x,y]”:

IMPORTANTE: Existem muitos modos de criar subplots, este é apenas um dos exemplos.

# Plot 1 (first row, first column)
# Plot the image
img1 = axs[0,0].imshow(data, vmin=vmin1, vmax=vmax1, origin='upper', cmap=my_cmap1)
# Add a colorbar
plt.colorbar(img1, extend='both', orientation='vertical', pad=0.05, fraction=0.05, ax=axs[0,0])
# Add a title
axs[0,0].set_title('Example 1 - Joining 2 cmaps', fontweight='bold', fontsize=10, loc='center')

Abaixo, o script completo até o momento:

# Training: Python and GOES-R Imagery: Script 8 - Custom Colormaps - Enhancing IR Channels
#-----------------------------------------------------------------------------------------------------------
# Required modules
from netCDF4 import Dataset          # Read / Write NetCDF4 files
import matplotlib.pyplot as plt      # Plotting library
from datetime import datetime        # Basic Dates and time types
from matplotlib import cm            # Colormap handling utilities
import numpy as np                   # Scientific computing with Python
from utilities import loadCPT        # Import the CPT convert function
#-----------------------------------------------------------------------------------------------------------
# COLORMAP EXAMPLE 1
# Custom colormap joining matplotlib colormaps
# Available matplotlib colormaps: https://matplotlib.org/stable/tutorials/colors/colormaps.html
vmin1 = -80                                                       # Min. value
vmax1 = 40                                                        # Max. value
gray_cmap = cm.get_cmap('gray_r', 120)                            # Read the reversed 'gray' cmap
gray_cmap = gray_cmap(np.linspace(0, 1, 120))                     # Create the array
jet_cmap  = cm.get_cmap('jet_r', 40)                              # Read the reversed 'jet' cmap 
jet_cmap  = jet_cmap(np.linspace(0, 1, 40))                       # Create the array
gray_cmap[:40, :] = jet_cmap                                      # Join both cmaps arrays
my_cmap1 = cm.colors.ListedColormap(gray_cmap)                    # Create the custom colormap
#-----------------------------------------------------------------------------------------------------------
# COLORMAP EXAMPLE 2
# Creating the INPE DISSM IR colormap
# Online color picker: https://imagecolorpicker.com/
vmin2 = -80                                                       # Min. value
vmax2 = 40                                                        # Max. value
gray_cmap = cm.get_cmap('gray_r', 120)                            # Read the reversed 'gray' cmap
gray_cmap = gray_cmap(np.linspace(0, 1, 120))                     # Create the array
colors = ["#ffa0ff", "#0806ff", "#3bcfff", "#feff65", "#ff7516"]  # Custom colors
my_colors = cm.colors.ListedColormap(colors)                      # Create a custom colormap
my_colors = my_colors(np.linspace(0, 1, 50))                      # Create the array
gray_cmap[:50, :] = my_colors                                     # Join both cmaps arrays
my_cmap2 = cm.colors.ListedColormap(gray_cmap)                    # Create the custom colormap 
#-----------------------------------------------------------------------------------------------------------
# COLORMAP EXAMPLE 3
# Creating a linear colormap
# Online color picker: https://imagecolorpicker.com/
colors = ["#bc8462", "#ae656f", "#a44a79", "#962e97", "#6158c5", "#2b8ffb", "#5fcdff", "#94fff0", "#a5ff94", "#fff88c", "#ffbf52", "#ec7b27", "#b84827", "#a1333d", "#bd5478", "#cc6a99", "#d982b8"]
my_cmap3 = cm.colors.LinearSegmentedColormap.from_list("", colors)# Create a custom linear colormap
vmin3 = -80                                                       # Min. value
vmax3 = 40                                                        # Max. value
#-----------------------------------------------------------------------------------------------------------
# COLORMAP EXAMPLE 4
# Converts a CPT file to be used in Python
# CPT archive: http://soliton.vm.bytemark.co.uk/pub/cpt-city/
cpt = loadCPT('IR4AVHRR6.cpt')                                    # Load the CPT file   
my_cmap4 = cm.colors.LinearSegmentedColormap('cpt', cpt)          # Create a custom linear colormap
vmin4 = -103.0                                                    # Min. value
vmax4 = 84.0                                                      # Max. value
#-----------------------------------------------------------------------------------------------------------
# Open the GOES-R image
# Download files at this link: http://home.chpc.utah.edu/~u0553130/Brian_Blaylock/cgi-bin/goes16_download.cgi
file = Dataset("OR_ABI-L2-CMIPF-M6C13_G16_s20191981200396_e20191981210116_c20191981210189.nc")
# Get the pixel values
data = file.variables['CMI'][:] - 273.15
#-----------------------------------------------------------------------------------------------------------
# Choose the plot size (width x height, in inches)
fig, axs = plt.subplots(2,2, figsize=(10,10)) # 2 rows x 2 columns
#-----------------------------------------------------------------------------------------------------------
# Plot 1 (first row, first column)
# Plot the image
img1 = axs[0,0].imshow(data, vmin=vmin1, vmax=vmax1, origin='upper', cmap=my_cmap1)
# Add a colorbar
plt.colorbar(img1, extend='both', orientation='vertical', pad=0.05, fraction=0.05, ax=axs[0,0])
# Add a title
axs[0,0].set_title('Example 1 - Joining 2 cmaps', fontweight='bold', fontsize=10, loc='center')
#-----------------------------------------------------------------------------------------------------------
# Plot 2 (first row, second column)
# Plot the image
img2 = axs[0,1].imshow(data, vmin=vmin2, vmax=vmax2, origin='upper', cmap=my_cmap2)
# Add a colorbar
plt.colorbar(img2, label='Brightness Temperatures (°C)', extend='both', orientation='vertical', pad=0.05, fraction=0.05, ax=axs[0,1])
# Add a title
axs[0,1].set_title('Example 2 - Custom + cmap', fontweight='bold', fontsize=10, loc='center')
#-----------------------------------------------------------------------------------------------------------
# Plot 3 (second row, first column)
# Plot the image
img3 = axs[1,0].imshow(data, vmin=vmin3, vmax=vmax3, origin='upper', cmap=my_cmap3)
# Add a colorbar
plt.colorbar(img3, extend='both', orientation='vertical', pad=0.05, fraction=0.05, ax=axs[1,0])
# Add a title
axs[1,0].set_title('Example 3 - Only custom colors', fontweight='bold', fontsize=10, loc='center')
#-----------------------------------------------------------------------------------------------------------
# Plot 4 (second row, second column)
# Plot the image
img4 = axs[1,1].imshow(data, vmin=vmin4, vmax=vmax4, origin='upper', cmap=my_cmap4)
# Add a colorbar
plt.colorbar(img4, label='Brightness Temperatures (°C)', extend='both', orientation='vertical', pad=0.05, fraction=0.05, ax=axs[1,1])
# Add a title
axs[1,1].set_title('Example 4 - CPT file', fontweight='bold', fontsize=10, loc='center')
#-----------------------------------------------------------------------------------------------------------
# Save the image
plt.savefig('Image_08.png')
# Show the image
plt.show()

VOCÊ PODE BAIXAR O SCRIPT ACIMA (Script_08.py) NESTE LINK

Execute o “Script_08.py” com o seguinte comando:

python Script_08.py

O seguinte plot será criado:

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

Clique na imagem abaixo para fazer o download do pack de arquivos do curso (scripts vistos na atividade pré-curso + todos os demais scripts e arquivos auxiliares.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s