Embed Size (px)
Citation preview
INPE-10436-NTC/361
MODELAGEM DE IRRADIÂNCIA SOLAR USANDO OCONJUNTO DE FERRAMENTAS IMAGE PROCESSING
WORKBENCH (IPW)
Brenner Stefan Gomes Silva
Nota Técnica
INPESão José dos Campos
2004
Índice
1. Introdução........................................................................................................... 12. Materiais e Dados Preliminares....................................................................... 12.1. Modelo Digital de Elevação (MDE) ................................................................ 12.2. Mapa de Albedo................................................................................................. 12.3. Obtenção de Parâmetros Atmosféricos (6S) .................................................. 23. Modelagem de Irradiância Solar ..................................................................... 34. Referências Bibliográficas............................................................................... 17
1
1. Introdução
Apresentam-se os procedimentos necessários para simulação da irradiância de
ondas curtas sobre a superfície e em condições de céu livre de nuvens. Tal
modelagem implica em calcular a influência da posição da terra sobre o valor da
irradiância exoatmosférica, a transferência radiativa na atmosfera e os efeitos
topográficos da área de estudo. O conjunto de ferramentas Image Processing
Workbench, ou IPW (Frew, 1990) disponibiliza funções para isso. Como exemplo
apresenta-se aqui a modelagem, ou simulação, da irradiância realizada em um
ambiente montanhoso.
2. Materiais e Dados Preliminares
É utilizado aqui o conjunto de ferramentas IPW em ambiente Unix. A versão 1.1 do
IPW pode ser obtida em http://www.icess.ucsb.edu/~ipw2/.
2.1. Modelo Digital de Elevação (MDE)
A interferência topográfica sobre a irradiância é calculada a partir do MDE. É
necessária a grade regular correspondente ao MDE em formato ASCII (texto) sem
cabeçalho. Para mapeamentos regionais, como em superfícies mamelonares, é
sugerida a interpolação por krigagem (Valeriano, 2002), por evitar artefatos do
MDE que formariam sombras inexistentes.
2.2. Mapa de Albedo
O mapa de albedo pode ser obtido a partir de imagens de satélite (Dugay e
LeDrew, 1992). Requer-se a priori correção radiométrica e posterior transformação
dos valores digitais em reflectância (Luiz et al, 2003). A correção atmosférica pode
ser realizada por diferentes métodos, p.ex. Dark Object Subtraction (Chavez, 1998)
ou pelo 6S, ou Second simulation of the satellite signal in the solar spectrum. A correção
topográfica pode ser realizada com o auxílio do MDE (Felicísimo e Garcia-
2
Manteca, 1989). O mapa final de albedo é uma grade regular em formato ASCII
(texto) sem cabeçalho.
2.3. Obtenção de Parâmetros Atmosféricos (6S)
A formulação da transferência radiativa incorporada no IPW requer os seguintes
parâmetros de entrada: espessura óptica da atmosfera (τ), albedo de espalhamento
simples (ω) e fator de assimetria do espalhamento (g). Tais parâmetros dependem
do comprimento de onda (λ) e devem ser obtidos para o respectivo intervalo
espectral da simulação. Utiliza-se aqui um modelo padronizado extraído do
algoritmo 6S, ou Second simulation of the satellite signal in the solar spectrum (Vermote
et al., 1997). Esta opção é útil quando não há dados coletados em campo para uma
determinada data, o que seria mais adequado (Zullo, 1994). Para o cálculo de
irradiância estes dados estão associados a uma altitude (t).
Nesta nota os parâmetros atmosféricos (τ, ω, g) foram divididos em intervalos de
comprimento de onda, para o cálculo particionado da irradiância. Esses valores são
apresentados na Tabela 1.
TABELA 1. Parâmetros atmosféricos utilizados para setembro e janeiro.
Comp.onda (λ) Esp.ópt (τ) alb.esp (ω) par.assim (g)
0,3 a 0,4 0,84 0,97 1,000,4 a 0,5 0,41 0,95 0,970,5 a 0,6 0,25 0,93 0,790,6 a 0,7 0,18 0,92 0,670,7 a 0,9 0,12 0,88 0,530,9 a 1,3 0,07 0,82 0,401,3 a 2,2 0,04 0,76 0,292,2 a 3 0,02 0,77 0,21
FONTE: 6S (Vermote et al., 1997).
3
3. Modelagem de Irradiância Solar
A função “toporad” do IPW é utilizada para calcular a distribuição espacial da
irradiância solar sobre uma área de interesse, dispondo inicialmente do Modelo
Digital de Elevação (MDE), de um mapa de albedo e parâmetros atmosféricos. São
descritos, em cinco passos, os arquivos de entrada e saída e comandos que acessam
as funções necessárias para simulação da irradiância. Em seguida é apresentado
um script escrito para o processamento em lote da irradiância em diferentes datas e
horários do dia.
Os arquivos ASCII são convertidos para o formato de imagem, tal como é utilizado
no IPW, através do comando import. Estes e os demais comandos do IPW estão
mais bem descritos em na ajuda fornecida pelo próprio programa. Como o MDE a
ser importado tem uma amplitude altimétrica de 2.365 metros, optou-se pela
importação de uma imagem de 16 bits para conter essas informações. A Figura 1
apresenta uma amostra do MDE utilizado como exemplo no procedimentos para a
obtenção da irradiância simulada no IPW. A expressão de importação é a seguinte:
% import -l 1848 -s 3435 -i mde.txt > mde.ipw
4
FIGURA 1 - Amostra extraída da altimetria importada em formato de imagem noIPW (mde.ipw) para ilustrar os passos de simulação de irradiância.
Após a importação é necessária a criação de um cabeçalho contendo a resolução da
imagem (mde.ipw), que corresponde à resolução da malha altimétrica (MDE). Este
cabeçalho deve conter também as coordenadas do primeiro pixel da imagem,
estabelecido como o pixel da extremidade superior esquerda. As coordenadas da
imagem são utilizadas para os cálculos da posição do sol em seu movimento
aparente nas diferentes horas do dia. Estes dados são inseridos no cabeçalho da
imagem e a partir daí podem ser realizados os cálculos de sombreamento e do
ângulo de incidência da luz solar. O comando que insere o cabeçalho de
georreferenciamento é o mkgeoh, como no exemplo:
5
%mkgeoh -o 7511796.036,482460.000 -d 30,30 -u meters -c UTM mde.ipw >
mde0.ipw
Onde as coordenadas planas (UTM) do primeiro pixel (p[0,0]) localizado no canto
superior esquerdo são: X= 482460.000, Y= 7511796.036; e a resolução é de 30X30
metros.
Para a execução do comando toporad, no 5o passo, a imagem de entrada deve ainda
conter um cabeçalho especificando a posição do sol que pode ser criado a partir do
comando mksunh:
%mksunh -z 0.672 -a 102.24 mde0.ipw > mde0.ipw
Onde, -z é o cosseno do ângulo zenital solar; e -a é o azimute solar (em radianos
partindo do sul, positivo a leste). Ambos são obtidos pelo comando sunang para 21
de junho de 2003. Esta data (solstício de inverno) é utilizada como exemplo aqui.
Os cabeçalhos inseridos no mde0.ipw são propagados nas matrizes resultantes de
cada passo.
10. Passo: Primeira aproximação da irradiância.
A obtenção dos valores de radiância solar direta e difusa para um dia de céu limpo
se dá a partir da função elevrad. Esta função recebe a altimetria da imagem de
entrada e escreve como saída uma imagem com duas bandas, uma para a
irradiância direta normal ao sol e outra para a irradiância difusa em uma superfície
horizontal. Esta execução geralmente precede a função toporad, que calcula a
irradiância global sobre o MDE.
Os parâmetros seguintes são necessários para a estimativa da transferência e
deposição de energia radiante na atmosfera a partir da radiância exoatmosférica:
6
-z: A altitude na qual foram obtidos os parâmetros atmosféricos (espessura óptica,
albedo de espalhamento simples e fator de assimetria do espalhamento).
Corresponde a altitude de medição destes parâmetros, caso realizada por um
fotômetro. Na utilização de parâmetros padronizados opta-se pela altitude média.
-t : Define um valor para a espessura óptica (τ);
-w: Define um valor para o albedo de espalhamento simples (ω);
-g: Define um valor para o fator de assimetria do espalhamento (g);
-r: Define um valor de albedo médio da superfície;
-s : A irradiância provinda do sol e que chega no topo da atmosfera. Seu valor pode
ser calculado pelo comando solar. Este comando permite a obtenção da radiância
exoatmosférica para um certo intervalo espectral.
-u : O cosseno do angulo zenital solar, obtido através do comando sunang.
O particionamento dos parâmetros atmosféricos (τ, ω e g) faz com que a função
elevrad seja executada uma vez para cada intervalo de comprimento de onda. As
coordenadas centrais são utilizadas na execução dos comandos solar e sunang,
predecessores ao elevrad. Essas funções fornecem a radiância exoatmosférica e os
ângulos de azimute e zenital solar. Como referência indica-se as coordenadas
centrais da área de estudo (Lat. 22o30’ e Long 45o00’). Um exemplo para o intervalo
0,3 a 0,4 no dia 21/06/2003 as 10:00 horas é dado a seguir:
%solar -d 2003,6,21 -w .3,.4
96.7298 W m^-2 (Radiância exoatmosférica)
%sunang –t 2003,6,21,10,00 –b -22,30 –l -45,00 –z 180
GMT Sat Jun 21 12:30:00 2003
7
-z 47.784 –u 0.671858 –a 102.24
onde
-z: Ângulo zenital; -u:Cosseno(z) e -a: Azimute solar.
%elevrad -z 1550 -t .31 -w .93 -g .75 -r .13 -s 13.4467 -u 0.671858 mde0.ipw >
mde1.ipw
A decomposição de mde1.ipw permite a obtenção das imagens de irradiância
direta normal ao vetor solar (mde11.ipw) e irradiância difusa normal ao horizonte
(mde12.ipw). As instruções do IPW tratam por bandas as imagens que formam
uma imagem composta. A separação destas bandas, ou imagens, é feita pelo
comando demux. Exemplo:
%demux -b 0 mde1.ipw > mde11.ipw
%demux -b 1 mde1.ipw > mde12.ipw
A Figura 2 apresenta a imagem mde11.ipw e seu respectivo histograma e a Figura
12 para mde12.ipw. O comando hist constrói o histograma para uma imagem de
entrada e prhist imprime as ocorrências estatísticas para cada valor digital. O
programa de visualização xgraph, disponível junto com o IPW serve para visualizar
histogramas como estes, que podem também ser analisados em planilhas.
FIGURA 2 - Irradiância direta calculada normal ao vetor solar pela função elevrad eseu respectivo histograma.
8
Figura 3 - Irradiância difusa calculada normal ao horizonte pela função elevrad e
seu respectivo histograma.
20.Passo: Azimute e declividade do terreno.
A partir da imagem de altimetria (MDE) são calculadas duas novas imagens
(bandas), uma contendo a exposição (azimute) e outra a declividade. A seguir estas
bandas são separadas em duas novas imagens:
%gradient mde0.ipw > mde2.ipw
%demux -b 0 mde2.ipw > mde21.ipw (declividade)
%demux -b 1 mde2.ipw > mde22.ipw (exposição)
Estas imagens são necessárias para o cálculo do fator de visibilidade do céu a partir
da função viewcalc. Esta função em conjunto com horizon é chamada pelo comando
viewf, posteriormente utilizado para o cálculo do fator de visibilidade do céu
usando diretamente o MDE (mdeg.ipw). A imagem composta mde2.ipw é
utilizada para calcular o ângulo de incidência (entre o vetor solar e a normal ao
terreno).
30.Passo: Ângulo de incidência.
O comando a seguir calcula a imagem contendo os valores de cosseno do ângulo
de incidência para cada pixel da imagem de altimetria (mdeg.ipw):
9
%shade -u 0.440168 -a 133.07 mde2.ipw > mde3.ipw
A transformação da irradiância direta normal ao vetor solar (mde11.ipw) para
irradiância direta normal ao relevo da superfície é feita multiplicando-se a imagem
mde11.ipw por mde3.ipw. O comando mux cria uma imagem composta destas
duas e o comando mult realiza a multiplicação (pixel a pixel) entre elas:
%mux mde3.ipw mde11.ipw | mult > mde3direta.ipw
Considerando que mde3.ipw contém somente os valores para os quais o sol está
diretamente incidindo, é considerado que este também assume a função de uma
mascara que exclui os pixels em sombra. Pode-se verificar este fato através do
histograma da imagem mde3.ipw.
%hist mde3.ipw | prhist > mde3.xg
%xgraph mde3.xg
A Figura 4 apresenta a imagem mde3.ipw calculada pela função shade e seu
respectivo histograma.
FIGURA 4 – Imagem contendo os valores de cosseno do ângulo de incidência
(mde3.ipw) obtidos pelo comando shade em relação à elevação eazimute solar para as coordenadas centrais da área.
É observado que 21394 pixels recebem o valor zero, ou seja, mais da metade dos
pixels não recebe iluminação direta. Logo esta imagem vai servir também como
uma mascara na multiplicação pelo valor da irradiância direta (mde11.ipw).
10
40. Passo: Fator de visibilidade do céu.
O comando viewf realiza o calculo do fator de visibilidade do céu (Vc) e de
configuração do terreno (Ct) a partir da matriz altimétrica (mde0.ipw) com
cabeçalho de coordenadas. O resultado é uma imagem composta (Vc e Ct) que pode
ser separada pelo comando demux. O cálculo dos horizontes é realizado pela
função horizon para 16 direções em cada célula da imagem.
%viewf mde0.ipw > mde4.ipw
%demux -b 0 mde4.ipw > mde41.ipw
%demux -b 1 mde4.ipw > mde42.ipw
O fator de visibilidade do céu é apresentado na Figura 5 junto a seu histograma
com valores de mínimo 0 e máximo 1. A imagem de configuração do terreno é
similar a do fator de visibilidade do céu, porém os valores mínimos e máximos no
histograma são invertidos (1 a 0):
11
FIGURA 5 - Imagem e histograma do fator de visibilidade do céu calculado sobre oMDE (mde0.ipw) pelo comando viewf.
A imagem da irradiância difusa é calculada pela multiplicação do fator de
visibilidade do céu com a irradiância difusa normal ao horizonte (mde12.ipw):
%mux mde12.ipw mde4.ipw| mult > mde4difusa.ipw
Até aqui são obtidas as porções direta (Figura 6) e difusa provinda do céu (Figura
7) da irradiância para cada ponto (30x30) da área de estudo.
FIGURA 6 - Imagem e histograma da irradiância direta normal a superfície(mde3direta.ipw).
12
FIGURA 7 - Imagem e histograma da irradiância difusa provinda do céu(mde4difusa.ipw).
50.passo: Irradiância global.
As imagens da Figura 6 e Figura 7 são complementares às calculadas pela função
toporad. Esta função calcula a distribuição topográfica da radiação solar para um
único horário, utilizando como entrada uma única imagem composta de 6 outras
imagens (bandas): irradiância direta e difusa (mde1.ipw), ângulo de iluminação
local (mde2.ipw), fator de visibilidade do céu e de configuração do terreno
(mde4.ipw) e a imagem referente ao albedo da superfície. A partir deste último
dado os valores de irradiância global incluem ainda a fração da radiação difusa
que é refletida pelo terreno. A matriz com os valores de albedo é importada tal
como o MDE (ASCII sem cabeçalho).
A imagem composta é formada pelo comando mux e em seguida a irradiância
global é calculada através do comando toporad. A Figura 18 apresenta o resultado
do comando toporad.
%mux marirradg.ipw marshade.ipw marviewf.ipw albedoiso.ipw | toporad >
mde5global.ipw
13
FIGURA 8 - Imagem de irradiância global obtida pelo comando toporad(mde5global.ipw) e seu respectivo histograma;
Comparando-se os valores tabulares dos pixels, verifica-se uma pequena diferença
nos valores da soma da irradiância direta e difusa com a irradiância global gerada
por toporad. Isto se dá devido ao acréscimo dado pela irradiância refletida pela
superfície.
%primg -a -i mde5global.ipw > mde5global.prt
%primg -a -i mde3direta.ipw > mde3direta.prt
%primg -a -i mde4difusa.ipw > mde4difusa.prt
mde5global.prt = mde3direta.prt + mde4difusa.prt + Albedo
Uma grade é considerada imagem no IPW e para transforma-la novamente em
uma grade regular (ASCII) utiliza-se o comando primg. Os valores de cada célula
estão quantizados linearmente de 0 a 65.536 (16 bits). Recuperam-se os dados de
irradiância (em W/m2) através de uma transformação linear, usando os valores
mínimo e máximo. Estes podem ser extraídos do cabeçalho da matriz através do
comando prhdr. As imagens associadas à algumas das grades de irradiância são
apresentadas na Figura 9.
14
Arquivos de execução em lote (Shell scripts) podem ser utilizados em ambiente
Unix para a modelagem de irradiância global em diferentes datas e momentos do
dia. Um exemplo destes arquivos, sem a repetição de comandos para cada horário,
é apresentado no a seguir.
#csh//22/12/2003//07:30 -u 0.377547 -a 72.777
gradient MDE16g.ipw | shade -u 0.377547 -a 72.777 > d070shad.ipw
elevrad -z 1500 -u 0.377547 -t 0.31 -w 0.93 -g 0.75 -r 0.115 -s 14.3514 MDE16g.ipw > d070a1.ipwelevrad -z 1500 -u 0.377547 -t 0.26 -w 0.96 -g 0.74 -r 0.115 -s 45.4167 MDE16g.ipw > d070a2.ipwelevrad -z 1500 -u 0.377547 -t 0.30 -w 0.96 -g 0.73 -r 0.115 -s 57.8211 MDE16g.ipw > d070a3.ipwelevrad -z 1500 -u 0.377547 -t 0.27 -w 0.96 -g 0.73 -r 0.115 -s 88.4045 MDE16g.ipw > d070a4.ipw
mux d070a1.ipw d070shad.ipw MDEviewf.ipw albedos.ipw | toporad > d070a1top.ipwmux d070a2.ipw d070shad.ipw MDEviewf.ipw albedos.ipw | toporad > d070a2top.ipwmux d070a3.ipw d070shad.ipw MDEviewf.ipw albedos.ipw | toporad > d070a3top.ipwmux d070a4.ipw d070shad.ipw MDEviewf.ipw albedos.ipw | toporad > d070a4top.ipw
\rm d070a1.ipw d070a2.ipw d070a3.ipw d070a4.ipw
elevrad -z 1500 -u 0.377547 -t 0.19 -w 0.96 -g 0.71 -r 0.115 -s 437.614 MDE16g.ipw > d070b.ipwelevrad -z 1500 -u 0.377547 -t 0.11 -w 0.94 -g 0.69 -r 0.115 -s 456.546 MDE16g.ipw > d070c.ipwelevrad -z 1500 -u 0.377547 -t 0.04 -w 0.79 -g 0.69 -r 0.115 -s 244.049 MDE16g.ipw > d070d.ipwelevrad -z 1500 -u 0.377547 -t 0.03 -w 0.68 -g 0.68 -r 0.115 -s 16.9692 MDE16g.ipw > d070e.ipw
mux d070b.ipw d070shad.ipw MDEviewf.ipw albedos.ipw | toporad > d070btop.ipwmux d070c.ipw d070shad.ipw MDEviewf.ipw albedos.ipw | toporad > d070ctop.ipwmux d070d.ipw d070shad.ipw MDEviewf.ipw albedos.ipw | toporad > d070dtop.ipwmux d070e.ipw d070shad.ipw MDEviewf.ipw albedos.ipw | toporad > d070etop.ipw
\rm d070b.ipw d070c.ipw d070d.ipw d070e.ipw d070shad.ipw
\\add070mux d070a1top.ipw d070a2top.ipw d070a3top.ipw d070a4top.ipw d070btop.ipw d070ctop.ipw d070dtop.ipwd070etop.ipw > d070top.ipwlincom -c 1 d070top.ipw > d070.ipw\rm d070top.ipw d070a1top.ipw d070a2top.ipw d070a3top.ipw d070a4top.ipw\rm d070btop.ipw d070ctop.ipw d070dtop.ipw d070etop.ipwprimg -a -r -i d070.ipw > d070.sprcat MNThdr d070.spr > d070c.sprcat d070c.spr MNTend > d070.spr\rm d070c.sprdos2unix d070.spr d070.spr
...
15
A figura 10 ilustra a irradiância obtida em 5 momentos (8, 10, 12, 14 e 16 horas, de
cima para baixo) dos dias de solstícios de verão (esquerda) e inverno (direita).
FIGURA 9 –Irradiância simulada para o solstício de verão e o de inverno
O fluxograma (modelo GEO-OMT) que esquematiza os passos para a simulação de
irradiância, tal como foi feito no IPW, é apresentado na Figura 10.
16
FIGURA 10 - Fluxograma detalhado dos passos executados no IPW para a
obtenção da irradiância de superfície.
17
4. Referências Bibliográficas
Felicíssimo, A. M.; Garcia-Manteca, P. Corrección del efecto topográfico en las
imagenes Landsat mediante el uso de un modelo digital de elevaciones. In:
Reunión Científica del Grupo de Trabajo en Teledetección, 3., 1989, Madrid,
España. Anais... Madrid: Asociación Española de Teledetección, 1989. p. 209-216.
Frew, J. The Image Processing Workbench. 1990. 382p. Phd. Thesis (Department
of Geography) - University of California, Santa Barbara, USA. 1990.
Luiz, A. J. Gürtler, S.; Gleriani, J. M.; Epiphânio, J. C. N.; Campos, R. C.
Reflectância a partir do número digital de imagens ETM+. In: Simpósio Brasileiro
de Sensoriamento Remoto, 11., 2003, Belo Horizonte. Anais… São José dos
Campos: INPE, 2003. p.2071-2078. CD ROM
Meador, W. E. Weaver, W. R. Two-stream approximations to radiative transfer in
planetary atmospheres: a unified description of existing methods and a new
improvement. Journal of Atmospheric Science, v. 37, p. 630-643, 1980.
Valeriano, M. M. Modelos digitais de elevação de microbacias elaborados com
krigagem . São José dos Campos: Instituto Nacional de Pesquisas Espaciais, 2002.
54 p. (INPE-9364-RPQ/736)
Vermote, E.; Tanré, D.; Deuzé, J.; Herman, M. e Morcette, J. Second simulation of
the satellite signal in the solar spectrum: An overview. IEEE Transactions on
Geoscience and Remote Sensing, v. 35, p. 675-686, 1997.
