ANÁLISE GEOESTATÍSTICA DE DADOS METEOROLÓGICOS DO ESTADO DO PARANÁ UTILIZANDO UM SOFTWARE LIVRE CLÁUDIO MARCHAND KRÜGER

Professor e Coordenador - Engenharia Civil - UnicenP/Centro Universitário Positivo [email protected]

dadaVinci Vinci, Curitiba, , Curitiba,v. v.2 2, n., n.1,1,p.p.79-86, 1-208,2005 2005

87

ANÁLISE GEOESTATÍSTICA DE DADOS METEOROLÓGICOS ...

RESUMO No presente trabalho são realizadas análises geoestatísticas de dados obtidos de estações meteorológicas distribuídas no território do Estado do Paraná. O software livre R foi utilizado em todas as etapas do trabalho. Para as análises geoestatísticas, foi utilizada a extensão geoR. A estatística espacial é o conjunto de métodos estatísticos nos quais a localização espacial desempenha um papel importante na análise de dados. O principal objetivo da geoestatística é modelar variações espaciais contínuas assumindo uma estrutura de correlação espacial da variável aleatória analisada. O software R e a extensão geoR demonstraram grande eficiência para este tipo de estudo, sendo uma alternativa viável em comparação com softwares comerciais especializados nesta área de conhecimento. Palavras-chave: Geoestatística, software livre, interpolação espacial.

ABSTRACT In this paper geostatistical analysis of meteorological data collected from stations distributed along the Paraná State are presented. The free software R was used in all phases of the work. For the geostatistical analysis, the geoR software extension was employed. Spatial Statistics is the set of statistical methods in which spatial locations play an important role in the analysis of data. The main aim of geostatistics is to model continuous spatial variation assuming a spatial correlation structure for the random variable under analysis. The free softwares R and geoR showed to be very efficient for this kind of study, and they are a viable alternative in comparison to commercial software dedicated to this area of knowledge. Key-words: Geostatistics, free software, spatial interpolation.

88

da Vinci , Curitiba, v. 2 , n. 1, p. 87-104, 2005

CLÁUDIO MARCHAND KRÜGER

ANÁLISE GEOESTATÍSTICA DE DADOS METEOROLÓGICOS DO ESTADO DO PARANÁ UTILIZANDO UM SOFTWARE LIVRE CLÁUDIO MARCHAND KRÜGER

1 INTRODUÇÃO No presente trabalho são realizadas análises geoestatísticas dos dados meteorológicos obtidos de estações meteorológicas distribuídas no território do Estado do Paraná. O software livre R foi utilizado em todas as etapas do trabalho. R é uma linguagem e ambiente computacional para cálculos estatísticos e elaboração de gráficos (www.r-project.org). Faz parte do projeto de softwares livres conhecido como GNU (www.gnu.org) e é similar ao software comercial S desenvolvido pelos Laboratórios Bell. O programa proporciona uma grande variedade de recursos de estatística e gráficos e é facilmente extensível através de módulos desenvolvidos por grupos de pesquisas em diferentes áreas da matemática e da computação distribuídos ao redor do mundo. No presente estudo, foi utilizada uma extensão do software R desenvolvida especialmente para análises geoestatísticas e previsões espaciais, chamada “geoR” (RIBEIRO e DIGGLE, 2001). Esta extensão foi desenvolvida no Departamento de Matemática e Estatística da Universidade de Lancaster, Reino Unido (www.est.ufpr.br/~paulojus/geoR) e tem recebido contribuições do Prof. Paulo Justiniano Ribeiro Jr., do Departamento de Estatística da UFPR. 2 O QUE É GEOESTATÍSTICA? A geoestatística tem suas origens em problemas ligados à estimativa de reservas minerais (CRESSIE, 1993). Seu desenvolvimento posterior acabou ocorrendo em grande parte de forma independente das correntes principais da estatística espacial. A estatística espacial é o conjunto de métodos estatísticos nos quais a localização espacial desempenha um papel importante na análise de dados. O principal objetivo da geoestatística é modelar variações espaciais contínuas assumindo uma estrutura de correlação espacial da variável analisada (DIGGLE e RIBEIRO, 2000). A geoestatística assume que a distribuição das diferenças de variáveis entre dois pontos amostrados é a mesma para toda a região em estudo, e que isto depende somente da distância entre eles e da orientação dos pontos. Ou seja, as diferenças existentes entre as diversas variáveis devem ser consistentes, porém não constantes.

da Vinci , Curitiba, v. 2 , n. 1, p. 87-104, 2005

89

ANÁLISE GEOESTATÍSTICA DE DADOS METEOROLÓGICOS ...

De maneira resumida, os passos de um estudo de geoestatística consistem em: - análise exploratória dos dados: coleta das amostras a serem analisadas, cálculo de estatísticas básicas e análises de tendências; - análise estrutural dos dados: análise e inferência de valores correlacionados de uma variável no espaço ou no tempo, chamada de modelagem do variograma; - realização de inferências: aplicação das técnicas de Krigagem (nome genérico dado ao processo de estimativa baseado na teoria das variáveis regionalizadas) simples e ordinária. Essa etapa é comumente chamada de Krigagem ou simulação. 2.1 Modelagem do Variograma A natureza estrutural de um conjunto de dados é definida a partir da comparação de valores tomados simultaneamente em dois pontos, segundo uma determinada direção. A função variograma é uma medida da variância das diferenças nos valores da variável regionalizada entre pontos separados por uma certa distância. Pontos mais próximos, por estarem correlacionados, terão essa variância pequena, aumentando à medida que os pontos se distanciam. Pode-se dizer que as técnicas geoestatísticas de estimativa são superiores por permitirem o cálculo do erro associado, chamado variância de krigagem. 3 DADOS UTILIZADOS 3.1 Dados meteorológicos Foram utilizados os dados de temperatura média (graus Celsius), umidade relativa do ar média (%) e total de horas de insolação média mensal de 55 estações meteorológicas do Estado do Paraná e próximas em São Paulo e Santa Catarina. Os dados meteorológicos foram fornecidos pela COPEL (Sistema MET) e também obtidos de publicações oficiais (Boletim Meteorológico, 1988; Brasil, 1969 e 1992; Santa Catarina, 1986). A figura 1 mostra o mapa do Estado do Paraná e a localização das estações meteorológicas utilizadas. Figura 1 – Localização das estações meteorológicas

90

da Vinci , Curitiba, v. 2 , n. 1, p. 87-104, 2005

CLÁUDIO MARCHAND KRÜGER

4 ANÁLISE GEOESTATÍSTICA 4.1 Análise de variogramas Uma função de fundamental importância em geoestatística é a função denominada de variograma. O variograma de um processo estocástico Y(x) é a função (1) onde x e x´ representam as coordenadas genéricas de duas observações quaisquer Y(x) e Y(x´) em um espaço n-dimensional. Para um certo conjunto de dados observados (Y i , Xi ): i=1, ... n, os valores esperados teóricos na equação acima são substituídos pelas médias das diferenças ao quadrado entre pares de observações Yi e Yj cujas localizações são dadas pelas coordenadas genéricas xi e xj . Na análise, supõe-se inicialmente que Y (x), é estacionário. O variograma empírico de um conjunto de dados (Y i , xi ): i=1, ... n é o conjunto de pontos (uij , vij ): j > i onde u ij =|| xi - xj|| é a distância euclidiana entre xi e x j e (2) A aplicação da equação acima ao conjunto de todos os dados (Y i , Xi ): i=1, ... n dá origem a uma nuvem de pontos com variabilidade muito grande para a modelagem teórica do processo. Por esta razão, é mais útil suavizar a nuvem de pontos do variograma empírico através do cálculo das médias de v ij para grupos de pontos que estejam a uma certa distância de separação uij .

da Vinci , Curitiba, v. 2 , n. 1, p. 87-104, 2005

91

ANÁLISE GEOESTATÍSTICA DE DADOS METEOROLÓGICOS ...

Figura 2 – Variogramas empíricos e suavizados

92

da Vinci , Curitiba, v. 2 , n. 1, p. 87-104, 2005

CLÁUDIO MARCHAND KRÜGER

semivariance

semivariance

semivariance

Nas três figuras anteriores, as linhas tracejadas são “envelopes” dos variogramas, obtidos através da simulação por troca de posições dos pontos na amostra. Desta forma, a faixa mostrada no envelope mostra uma larga variação possível do variograma amostral. Esta grande variabilidade pode ser parcialmente explicada pelo número relativamente pequeno de pontos amostrais (55). Ainda assim, evidencia-se a existência de uma estrutura de correlação espacial nas três variáveis analisadas. Uma forma alternativa de análise de variogramas é a consideração de uma direção preferencial de variação da variável analisada. Por exemplo, ao se estudar a dispersão espacial da poluição produzida por uma indústria, vale a pena concentrar a análise na direção do vento dominante no local e dar menos importância às variações em outras direções. Pela análise dos variogramas direcionais, pode-se notar uma anisotropia razoavelmente bem definida havendo maior variação na variabilidade espacial quando se consideram ângulos de orientação de 45o e 135o. Estas orientações devem-se às características do relevo do Estado (mais marcante na direção Leste-Oeste) e da variação climática (com grande variabilidade no sentido NorteSul). Devido ao número reduzido de pontos amostrais, tornando difícil a avaliação do ângulo de anisotropia mais adequado baseando-se nos variogramas direcionais, as análises seguintes não consideraram o aspecto de anisotropia.

Figura 3 – Variogramas direcionais

da Vinci , Curitiba, v. 2 , n. 1, p. 87-104, 2005

93

ANÁLISE GEOESTATÍSTICA DE DADOS METEOROLÓGICOS ...

4.2 Ajuste de parâmetros do modelo geoestatístico O modelo geoestatístico adotado assume que as observações representam uma versão “com ruído” de um sinal que descreve a verdadeira variabilidade da variável aleatória na área de estudo. Um dos objetivos da análise geoestatística é prever o sinal na área ou alguma quantidade que possa ser descrita como uma função do sinal. Um procedimento bastante comum (mas não o mais preciso) para estimativa de parâmetros do modelo geoestatístico é ajustar uma curva ao variograma empírico (figura 4). Desta forma, os parâmetros do modelo podem ser lidos diretamente do modelo teórico ajustado. O jargão normalmente utilizado neste caso é chamar o intercepto da curva de nugget ou “pepita”, fazendo referência às origens da geoestatística na área de mineração. A diferença entre a assíntota horizontal do variograma e o nugget é chamado de sill (ou “patamar”) e a distância na qual o variograma teórico atinge seu máximo é conhecido como range (ou “alcance”). Estes nomes correspondem aos parâmetros (nugget), (sill) e (range), respectivamente. O parâmetro normalmente é expresso multiplicado por uma constante, dependendo do modelo utilizado. Para modelos onde o alcance é infinito, utiliza-se um “alcance prático”, que seria o ponto onde o variograma atinge 95% da assíntota. É interessante notar que, mesmo para uma distância igual a zero, normalmente os dados amostrais apresentam um resíduo de variabilidade, chamado pepita. Este resíduo pode ser atribuído a erros de observação ou a descontinuidades locais da variável analisada, por exemplo, a ocorrência de uma pepita de ouro em um campo de mineração. semivariância

Figura 4 – Componentes do variograma

No presente estudo, inicialmente, foram tentados ajustes “a sentimento”, para tentar encontrar valores iniciais coerentes dos parâmetros , e . Esta análise inicial revela-se importante, pois as análises subseqüentes, como ajuste por mínimos quadrados e máxima verossimilhança, usam os parâmetros iniciais nos procedimentos de otimização.

94

da Vinci , Curitiba, v. 2 , n. 1, p. 87-104, 2005

CLÁUDIO MARCHAND KRÜGER

4.2.1 Famílias de modelos de funções de correlação Os variogramas teóricos usados no presente estudo incluem o “efeito pepita”. No caso estacionário, o variograma com este efeito pode ser representado por (DIGGLE e RIBEIRO, 2000): (3) Na expressão acima é uma função de correlação espacial. Nem toda função teórica pode ser candidata a representar . Além de alguns requisitos não-triviais do ponto de vista matemático não citados aqui, usualmente exige-se que o modelo usado para a função de correlação incorpore os seguintes atributos (DIGGLE e RIBEIRO, 2000): seja monotonamente não-crescente com u (a correlação entre duas observações 1. decresce com o aumento da distância entre dois pontos amostrais); 2. quando (a correlação tende para zero para grandes distâncias entre pontos); tende para zero 3. pelo menos um parâmetro no modelo deve controlar a taxa na qual (visto que a distância de separação na qual a correlação torna-se desprezível não é conhecida de antemão). Adicionalmente, é útil incluir no modelo alguma flexibilidade na forma geral da função de correlação. Desta forma, um modelo paramétrico para a função de correlação deve ter um ou dois parâmetros, e um modelo para o variograma três ou quatro (dois parâmetros de correlação mais os dois componentes da variância ( e ) (DIGGLE e RIBEIRO, 2000). A seguir, algumas famílias de funções que atendem a estes requisitos e que foram utilizadas na presente análise: Família Esférica (4)

Como esta família depende apenas do parâmetro , não há flexibilidade para a forma. Adicionalmente, segundo DIGGLE e RIBEIRO (2000), a lógica de sua aplicação física é questionável em problemas no espaço bidimensional. Família Exponencial (5) onde > 0 e . A função de correlação exponencial corresponde ao caso em que k=1. Para k=2 , esta função chama-se função de correlação gaussiana. Esta família freqüentemente produz um ajuste qualitativamente razoável à estrutura de correlação dos da Vinci , Curitiba, v. 2 , n. 1, p. 87-104, 2005

95

ANÁLISE GEOESTATÍSTICA DE DADOS METEOROLÓGICOS ...

dados espaciais, mas as previsões baseadas neste modelo tendem a ser pouco robustas a pequenas variações sobre o modelo escolhido (DIGGLE e RIBEIRO, 2000). Família de Matérn (6) onde ( , k ) são parâmetros, G representa a função gama e Kk representa uma função de bessel modificada de ordem k. A família é válida para > 0 e k > 0. O caso k=0,5 corresponde exatamente à função de correlação exponencial, . A função gaussiana corresponde ao caso limite no qual . Segundo RIBEIRO et al. (2003), A família de Matérn é particularmente atrativa, pois o valor do parâmetro k controla a suavidade do processo do sinal subjacente. Desta forma, conclui-se que a família de Matérn é provavelmente a melhor escolha para função de correlação para uso geral, por sua flexibilidade e simplicidade (apenas dois parâmetros). As figuras a seguir mostram os variogramas empíricos de temperatura, umidade e insolação, juntamente com os ajustes (a sentimento) de funções de correlação exponencial e esférica e também os ajustes por mínimos quadrados ponderado (WLS) e ordinário (OLS). São mostrados também os ajustes da função de Matérn e a estimativa de parâmetros diretamente (sem auxílio do variograma) por máxima verossimilhança.

96

da Vinci , Curitiba, v. 2 , n. 1, p. 87-104, 2005

CLÁUDIO MARCHAND KRÜGER

Figura 5 – Semivariogramas teóricos

4.3 Uso do modelo para previsão e interpolação espacial (Kriging) Na maioria das aplicações de geoestatística, o maior interesse está na realização de inferências sobre o sinal S(x) , tais como prever o valor de S(x) sobre a área de estudo, ou estimar a probabilidade de que S(x) esteja acima de um certo valor limite. O problema da previsão então reduz-se ao estudo da distribuição condicionada de S(x) dados os valores observados y, e para tal considera-se uma malha de previsão sobre a região. Considerando um valor isolado S0 em um local genérico x0 , (S0 , Y ) possui, sob a hipótese de um modelo gaussiano, uma distribuição condicionada de previsão [ S0 | Y ] com média e variância dadas por: (7) (8) onde: r é um vetor (n x 1) com os elementos entre x e xj ;

onde || xi - xj || é a distância euclidiana

R é uma matriz (n x n) com os (i,j)-ésimos elementos

;

I é uma matriz identidade (n x n). Com todos os parâmetros do modelo considerados conhecidos, as equações acima são chamadas de equações de “kriging simples”. As fórmulas acima estão implementadas na função do R chamada krige.conv, cujo nome significa “kriging convencional”. As entradas para a função são o objeto geodata e coordenadas dos pontos amostrais e informações do modelo, incluindo valores dos parâmetros. A função calcula estimativas das médias e variâncias na malha de previsão.

da Vinci , Curitiba, v. 2 , n. 1, p. 87-104, 2005

97

ANÁLISE GEOESTATÍSTICA DE DADOS METEOROLÓGICOS ...

As figuras a seguir mostram os mapas com as previsões de temperatura, umidade e insolação, usando o modelo de Matérn. Por motivo de falta de espaço, serão omitidos os resultados do modelo de máxima verossimilhança. Em seguida, são mostrados os mapas com os erros de previsão e a seqüência de comandos do R. As coordenadas estão representadas em graus e as variáveis em graus Celsius (temperatura), porcentagem (umidade relativa) e horas (insolação).

Figura 6 – Mapas de previsão e erros de previsão

98

da Vinci , Curitiba, v. 2 , n. 1, p. 87-104, 2005

CLÁUDIO MARCHAND KRÜGER

5 CONCLUSÕES No presente trabalho foram analisados os dados meteorológicos de 55 estações distribuídas pelo Estado do Paraná e perto das fronteiras nos estados vizinhos de São Paulo e Santa Catarina. As seguintes conclusões foram obtidas da análise geoestatística dos dados de temperatura média, umidade relativa do ar média e insolação média mensal: - Os envelopes dos variogramas empíricos demonstraram a existência de estrutura de correlação espacial para as três variáveis meteorológicas analisadas. - A análise da anisotropia foi dificultada pela grande variabilidade dos variogramas direcionais causada pelo pequeno número de pontos amostrais, sendo desconsiderada nas análises subseqüentes. Ainda assim, evidenciou-se a existência de anisotropia nos dados das três variáveis analisadas. - Os parâmetros do modelo geoestatístico foram obtidos através do ajuste das funções de correlação esférica, exponencial e Matérn e também pela inferência direta aos dados através de estimativas por máxima verossimilhança. - As estimativas por máxima verossimilhança nem sempre se mostram visualmente as melhores, comparando-se os pontos do variograma empírico com a função teórica ajustada, mas os erros de previsão mostraram-se menores que os outros métodos, em alguns casos. - O modelo de Matérn foi o que apresentou os menores erros de previsão, entre as funções de correlação testadas. Pela literatura pesquisada, é o modelo mais vantajoso para uso geral entre as funções candidatas, por ser simples (apenas dois parâmetros) e flexível, adaptando-se melhor aos dados através da escolha do parâmetro k. - Os mapas com os erros de previsão demonstraram um resultado melhor para as variáveis temperatura e insolação média. A variável umidade relativa apresentou maiores desvios na previsão. - O software R e a extensão geoR demonstraram grande eficiência para este tipo de estudo, sendo uma alternativa viável em comparação com softwares comerciais especializados nesta área de conhecimento.

da Vinci , Curitiba, v. 2 , n. 1, p. 87-104, 2005

99

ANÁLISE GEOESTATÍSTICA DE DADOS METEOROLÓGICOS ...

REFERÊNCIAS BIBLIOGRÁFICAS BOLETIM METEOROLÓGICO. Dados meteorológicos do Estado de São Paulo. DAEE, 1988. BRASIL. Ministério da Agricultura. Escritório de Meteorologia. Normas Climatológicas. Rio de Janeiro, 1969. v.4 BRASIL. Ministério da Agricultura e Reforma Agrária. Normas Climatológicas: 19611990. Brasília, 1992. CRESSIE, N. Statistics for spatial data. New York: Wiley, 1993. DIGGLE, P. J.; RIBEIRO JÚNIOR, P. J. Model Based Geostatistics. In: SINAPE, 14., 2000, Caxambu. Anais... São Paulo: Associação Brasileira de Estatística, 2000. RIBEIRO JÚNIOR, P. J.; CHRISTENSEN, O. F.; DIGGLE, P. J. Geostatistical software : geoR and geoRglm. In:INTERNATIONAL WORKSHOP ON DISTRIBUTED STATISTICAL COMPUTING, 3. Anais... Viena: Áustria, 2003. RIBEIRO JÚNIOR, P. J.; DIGGLE, P. J. GeoR: a package for geostatistical analysis. R. News, v.1-2, June 2001. SANTA CATARINA. Gabinete de Planejamento e Coordenação Geral. Subchefia de Estatística, Geografia e Informática. Atlas do Estado de Santa Catarina. Rio de Janeiro: Aerofoto Cruzeiro, 1986.

100

da Vinci , Curitiba, v. 2 , n. 1, p. 87-104, 2005

CLÁUDIO MARCHAND KRÜGER

ANEXOS COMANDOS DA LINGUAGEM R ## ## Variogramas ## ## Temperatura temp.vario