Interpolação IDW e Krigagem de Dados TRMM
Elivaldo Rocha
10 de março de 2021
1 Descrição da Área de Estudo
Belém (inicialmente chamada Santa Maria de Belém do Pará ou Nossa Senhora de Belém do Grão Pará e frequentemente chamada de Belém do Pará) é um município brasileiro e capital do estado do Pará, situado na região Norte do país, a uma latitude 01º27’21" sul e longitude 48º30’16" oeste, fundado em 12 de janeiro de 1616 pelos portugueses, às margens da baía do Guajará.
É uma cidade histórica e portuária, localizada na Amazônia Oriental, ao extremo nordeste da maior floresta tropical do mundo, sendo a capital mais chuvosa do Brasil devido a seu clima equatorial, influenciada diretamente pela Amazônia. Belém possui uma área de 1 059,458 km² e uma altitude de dez metros ao nível médio do mar em média, estando a cerca de 2 140 km da capital federal, Brasília.
Alta pluviosidade, altitude próximo ao nivel médio do mar e praticamente plana são caracteristicas que favorecem alagamentos.
Mapa de localização.
2 Sensor TRMM 3B43
A precipitação é uma das variáveis meteorológicas mais importantes na região tropical, sendo sua distribuição espacial relevante dentre outras aplicações para o manejo de atividades agrícolas e tomadores de decisões de algumas atividades econômicas, tais como: a construção civil, empresas de energia elétrica, modelagem hidrológica e entre outras.
A Tropical Rainfall Measuring Mission (TRMM) é uma missão espacial conjunta entre a NASA e a Agência Nacional de Desenvolvimento Espacial do Japão, projetada para monitorar e estudar a precipitação tropical e subtropical e a liberação de energia associada.
A missão usa 5 instrumentos: Radar de precipitação (PR), TRMM Microwave Imager (TMI), Visible Infrared Scanner (VIRS), Clouds & Earths Radiant Energy System (CERES) e Lightning Imaging Sensor (LSI). O TMI e o PR são os principais instrumentos utilizados para a precipitação. Esses instrumentos são usados em um algoritmo que forma o conjunto de dados de calibração do Instrumento Combinado TRMM (TCI) (TRMM 2B31) para a Análise de Precipitação Multi-satélite TRMM (TMPA), cujas médias de precipitação mensal TMPA 3B43 e médias diárias e sub-diárias (3h) TMPA 3B42 são provavelmente os produtos mais relevantes relacionados ao TRMM para pesquisas climáticas.
3B42 e 3B43 estão disponíveis em resolução espacial de 0,25°, aproximadamente 27,8km, cobrindo 50°N a 50°S de 1998 até o 2019.
3 Visualizaçao dos dados TRMM 3B43
Os dados originalmente tem extensão nas coordenadas 10ºS 60ºW e 02ºN 038ºW, com resolução espacial de 0.25º, o que equivale aproximadamente 27,8 Km. Nos valores dos pixeis encontram-se a média mensal em milímetros por mês desde 01 de Janeiro de 1998 até 31 de dezembro de 2019. Sistema de Coordenadas de Referência : Geográfico - WGS 84 - EPSG:4326.Site para obtenção dos dados:
Os dados vetorias são oriundos do Institudo Brasileiro de Geografia e Estatística - IBGE-2019
O objetivo é criar um arquivo raster cuja a resolusão espacial seja de 400 metros.
Visualização dos Dados.
4 Descrição Dos Métodos de Interpolação Espacial
O processo de interpolação espacial possibilita a criação de uma superfície estatiística. Essa superfície contínua é gerada apartir de pontos com valores conhecidos do fenômeno de interesse. Existe vários métodos de interpolação,como por exemplo, Kriging, polígonos de Thiessen, IDW. Estes podem ser classificados pela natureza determinista e geoestatística.
5 Método da Ponderação pelo Inverso da Distância (IDW)
A interpolação IDW é classificada como determinística e produz uma superfície contínua apartir de pontos conhecidos. A ponderação é atribuida conforme a distância entre ponto no qual se quer prever e o ponto amostrado. Sendo assim, pontos que estão mais distantes tem peso menor e influencia menos na previsão, pontos que estão mais próximos tem peso maior e, também, maior influencia na previsão.
Esse cáculo é realizado através da equação descrita por Mello et al. (2003) abaixo :
\[Z_{p} = \frac{\sum_{i=1}^{n}(\frac{1}{d_{i}^{p}}*z_{i})}{\sum_{i=1}^{n}(\frac{1}{d_{i}^{p}})}\] onde:
\(\mathit{Z_{p}}\) é o atributo interpolado,
\(\mathit{z_{i}}\) valor do atributo do i-ésimo ponto de amostragem,
\(\mathit{d_{i}}\) distância euclidiana entre o i-ésimo ponto de vizinhança e o ponto amostrado,
\(\mathit{p}\) é o expoente utilizado.
\(\mathit{n}\) número de amostras.
5.1 Exemplo
Faremos a interpolação para encontrar a precipitação em mm/mês do ponto z:
Representação para Cálculo de IDW.
Solução: \[Z_{p} = \frac{\sum_{i=1}^{n}(\frac{1}{d_{i}^{p}}*z_{i})}{\sum_{i=1}^{n}(\frac{1}{d_{i}^{p}})} = \frac{(\frac{4}{3^{2}}+\frac{5}{2^2}+\frac{8}{6^2})}{(\frac{1}{3^2}+\frac{1}{2^2}+\frac{1}{6^2})} = \frac{(\frac{23}{12})}{\frac{7}{18}}\approx 4.92_{mm/mes}\]
6 Efeito do Expoente
O expoente \(\mathit{p}\) atribui alguns efeitos na interpolação, e sua escolha deve ser tomada em virtude de métodos estatísticos como, por exemplo, a validação cruzada.
Vejamos :
Baixo, menor que 2 : destacam anomalias locais;
2: inverso do quadrado da distâcia, o mais usado;
Altos, entre 3 e 5: suavizam anomalias locais; e
maior ou igial 10 : estimam poligonos planos, muito semelhante ao método conhecido como polígos de Varonoi.
6.1 Ponderação em função da distância
Efeito do Expoente.
6.2 Animação
7 Avaliar a eficiência de modelos por validação cruzada (CROSS-VALIDATION)
A validação cruzada é uma técnica para avaliar a capacidade de generalização de um modelo, a partir de um conjunto de dados. Esta técnica é amplamente empregada em problemas onde o objetivo da modelagem é a predição. Busca-se então estimar o quão preciso é este modelo na prática, ou seja, o seu desempenho para um novo conjunto de dados.
O conceito central das técnicas de validação cruzada é o particionamento do conjunto de dados em subconjuntos mutuamente exclusivos, e posteriormente, o uso de alguns destes subconjuntos para a estimação dos parâmetros do modelo (dados de treinamento), sendo os subconjuntos restantes (dados de validação ou de teste) empregados na validação do modelo.
\[Ac_{f} = \frac{1}{v}\sum_{i=1}^{v}\epsilon _{y_{i},\hat{y}_{i}} = \frac{\sum_{i=1}^{v}(y_{i}-\hat{y_{i}})}{v}\] onde:
\(\mathit{v}\) é o número de dados de validação, e
\(\mathit{\epsilon _{y_{i},\hat{y}_{i}}}\) é o resíduo dado pela diferença entre o valor real \(\mathit{y_{i}}\) e o valor predito \(\mathit{\hat{y}_{i}}\).
Desta forma é possível inferir quantitativamente a capacidade de generalizações dos modelos.
7.1 Raiz do Erro Médio Quadrático (Root Mean square Error - RMSE)
\[RMSE = \sqrt{\frac{\sum_{i=1}^{v}(y_{i}-\hat{y_{i}})^2}{v}}\]
O desvio médio quadrático ( RMSD ) ou erro quadrático médio ( RMSE ) é uma medida frequentemente usada das diferenças entre os valores (valores de amostra ou população) previstos por um modelo ou um estimador e os valores observados. O RMSD representa a raiz quadrada do segundo momento amostral das diferenças entre os valores previstos e os valores observados ou a média quadrática dessas diferenças. Esses desvios são chamados de resíduos quando os cálculos são realizados sobre a amostra de dados que foi usada para a estimativa e são chamados de erros(ou erros de predição) quando calculado fora da amostra. O RMSD serve para agregar as magnitudes dos erros nas previsões para vários pontos de dados em uma única medida de poder preditivo. O RMSD é uma medida de precisão para comparar erros de previsão de diferentes modelos para um conjunto de dados específico e não entre conjuntos de dados, pois depende da escala.
O RMSD é sempre não negativo e um valor 0 (quase nunca alcançado na prática) indicaria um ajuste perfeito aos dados. Em geral, um RMSD menor é melhor do que um maior. No entanto, as comparações entre diferentes tipos de dados seriam inválidas porque a medida depende da escala dos números usados.
RMSD é a raiz quadrada da média dos erros quadrados. O efeito de cada erro no RMSD é proporcional ao tamanho do quadrado do erro; portanto, erros maiores têm um efeito desproporcionalmente grande no RMSD. Consequentemente, o RMSD é sensível a outliers.
Menor RMSE dado pelo Melhor Expoente.
fonte: https://ascelibrary.org/cms/asset/d0bd0f01-979c-4e7a-8a81-2f7b96c81f5d/figure6.gif
7.2 Estimativa de explicação da variância em Porcentagem
Análogo ao R² - O coeficiente de determinação
\[1 - \frac{V_{ariação}D_{os}R_{esíduos}}{V_{ariação}D_{ados}O_{bservados}}\]
É uma medida de ajuste de um modelo estatístico generalizado aos valores observados de uma variável aleatória. Está entre 0 e 1, por vezes sendo expresso em termos percentuais. desta forma, expressa a quantidade da variância dos dados que é explicada pelo modelo.
No caso dos dados pontuais do sensor SRTMM 3B43, quanto maior, mais explicativo é o modelo, ou seja, melhor ele se ajusta aos dados pontuais. Por exemplo, um R² = 0,82341 significa que o modelo explica 82,34% da variância dos dados interpolados a partir dos dados pontuais incluídas no modelo.
Lembrando que:
\[V_{ariação}D_{os}R_{esíduos} = \frac{\sum_{i=1}^{v}(\epsilon _{i}-\mu _{\epsilon})^2}{v}\] \(\mathit{\epsilon_{i}}\) é o resíduo obtido apartir de \(\mathit{y_{i}-\hat{y}_{i}}\) ;
\(\mathit{\mu_{\epsilon}}\) é a média dos resíduos.
\[V_{ariação}D_{ados}O_{bservados} = \frac{\sum_{i=1}^{v}(y _{i}-\mu _{y})^2}{v}\] \(\mathit{y_{i}}\) é o dado observado no ponto \(\mathit{i}\) ;
\(\mathit{\mu_{y}}\) é a média dos dados observados.
8 Resulta Obtido pelo Método IDW
A chuva não cai uniformemente distribuída por toda a bacia, principalmente quando se trata de grandes bacias. A precipitação tem como característica a grande variabilidade espacial e temporal. Note que ao aumentarmos o expoente p, a predição tende a se assemelhar com os dados originais do TRMM com resolução espacial de 0,25º. No entando, o interessante é visualizar a variação dentro do minicípio de Belém. Sendo assim, elegeu-se p igual a 3.75 , ponto no qual apartir dele o gráfico tende a se manter constante.
8.1 Escolha do Melhor Expoente
Análise do Expoente.
8.2 Tabela do Modelo
| Predição | observado | Resíduo | Longitude em UTM | latitude em UTM |
|---|---|---|---|---|
| 147.3248 | 137.7717 | -9.553116 | 4346819 | 10319583 |
| 146.5778 | 139.7026 | -6.875172 | 4374613 | 10319443 |
| 152.7525 | 142.7469 | -10.005581 | 4402408 | 10319310 |
| 164.3419 | 155.8157 | -8.526134 | 4430203 | 10319182 |
| 176.6915 | 178.5787 | 1.887253 | 4457998 | 10319060 |
| 186.8877 | 189.6336 | 2.745857 | 4485792 | 10318945 |
8.3 Modelo de Regressão Entre os Dados Preditos e os Dados Observados
Modelo de regressão entre os dados preditos e observados.
8.4 Dados Interpolados com Resolução de 400 metros
Mapa de localização.
8.5 Visualização com Open Street MaP
9 Krigagem
Poderá ser entendido como uma predição linear ou uma forma da Inferência bayesiana. Parte do princípio que pontos próximos no espaço tendem a ter valores mais parecidos do que pontos mais afastados. A técnica de Krigagem assume que os dados recolhidos de uma determinada população se encontram correlacionados no espaço. Isto é, se a precipitação média mensal de um ponto p numa cidade é x mm/mês , é muito provável que se encontrem resultados muito próximos de x quanto mais próximos se estiver do ponto p (princípio da geoestatística). Porém, a partir de determinada distância de p, certamente não se encontrarão valores aproximados de x porque a correlação espacial pode deixar de existir.
Considera-se o método de Krigagem do tipo BLUE (Best Linear Unbiased Estimator - Melhor Estimador Linear não-Viciado): é linear porque as suas estimativas são combinações lineares ponderadas dos dados existentes; é não enviesada pois procura que a média dos erros (desvios entre o valor real e o valor estimado) seja nula; é a melhor porque os erros de estimação apresentam uma variância (variância de estimação) mínima. O termo Krigagem abrange um conjunto de métodos, sendo os mais usuais os seguintes:
9.1 Tipos de Krigagem
9.1.1 Krigagem Simples
Assume que as médias locais são relativamente constantes e de valor muito semelhante à média da população que é conhecida. A média da população é utilizada para cada estimação local, em conjunto com os pontos vizinhos estabelecidos como necessários para a estimação.
9.1.2 Krigagem Normal
As médias locais não são necessariamente próximas da média da população usando-se apenas os pontos vizinhos para a estimação. É o método mais usado em problemas ambientais.
9.1.3 Co-krigagem
É uma extensão do anterior a situações em que duas ou mais variáveis são espacialmente dependentes e a variável que se quer estimar não está amostrada com a intensidade com que estão as outras variáveis dependentes, utilizando-se os valores destas e as suas dependências para estimar a variável requerida.
9.2 Conceitos Matemáticos
O Método de Krigagem serve-se de diversas teorias explanadas na estatística. No entanto, para deixarmos mais claras as teorias de estatística usadas e mais direcionadas ao escopo deste texto, explicar-se-á alguns conceitos.
9.2.1 Semi-Variância e Semi-Variograma
A semi-variância é a medida do grau de dependência espacial entre duas amostras. A magnitude da semi-variância entre dois pontos depende da distância entre eles, implicando em semi-variâncias menores para distâncias menores e semi-variâncias maiores para distâncias maiores. O gráfico das semi-variâncias em função da distância a um ponto é chamado de Semi-variograma. A partir de uma certa distância a semi-variância não mais aumentará com a distância e estabilizar-se-á num valor igual à variância média, dando a esta região o nome de silo ou patamar (sill). A distância entre o início do semi-variograma e o começo do silo recebe o nome de range ou amplitude ou alcance. Ao extrapolarmos a curva do semi-variograma para a distância zero, podemos chegar a um valor não-nulo de semi-variância. Este valor recebe o nome de Efeito Pepita (Nugget Effect).
9.2.2 Modelos de Variograma
No Método de Krigagem normalmente são usados quatro tipos de variogramas. Neles, são usadas as seguintes variáveis:
Características do Semi-variograma.
\(\displaystyle{\gamma_{(h)}}\): semivariância;
\(\displaystyle{c_{0}}\): efeito pepita (nugget)-variação ao acaso-Fatores não relacionados ao espaço e Erros de amostragem;
\(\displaystyle{c}\): contribuição estruturada(Structured Variability);
\(\displaystyle{\alpha}\): alcance (range)-até onde ocorre a autocorrelação espacial;
\(\displaystyle{c_{0}+c}\): patamar ou variância assintótica máxima (Total Variability);
\(\displaystyle{h}\): distância de separação;
Para distância \(h \leq \alpha\), situa-se o domínio estruturado - onde a autocorrelação espacial influência; e
Para distância \(h \geq \alpha\), situa-se o domínio aleatório, onde a diferença dos valores dos pontos amostrais não dependem mais da distância.
modelos matematicos de variogramas.
9.2.2.1 Linear
Este modelo não apresenta silo e é muito simples. Sua curva pode ser representada por:
\[\displaystyle{ \gamma_{(h)} = c_{0} + ch}\]
9.2.2.2 Esférico
A forma esférica é a mais utilizada e possui silo. Sua forma é definida por:
\[\displaystyle{\gamma_{(h)} = \begin{cases} c_{0} + c[1.5(\frac{h}{\alpha }) - 0.5(\frac{h}{\alpha })^3]~~, & ~Se~h < \alpha \\ c_{0} + c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~, & ~Se~h> \alpha \end{cases}}\]
9.2.2.3 Exponencial
A curva do variograma exponencial respeita a seguinte equação:
\[\displaystyle{\gamma_{(h)}=c_{0}+c(1-e^\frac{-h}{\alpha})}\]
9.2.2.4 Gaussiano
A forma gaussiana é dada por:
\[\displaystyle{\gamma_{(h)} = \begin{cases} c_{0}+c(1-e^{\frac{-h^2}{\alpha^2}}), & ~Se~h < \alpha \\ c_{0} + c~~~~~~~~~~~~~, & ~Se~h> \alpha \end{cases}}\]
10 O Método de Krigagem
A geoestatística foi utilizada para avaliar a variabilidade espacial dos atributos estudados, sendo a dependência espacial expressa através dos semivariogramas, estimado pela Equação:
\[\mathit{\gamma _{(h)} = \frac{1}{2N_{(h)}}\sum_{i=1}^{N_{(h)}}\left [ Z_{(X_{i})} - Z_{(X_{i} + h)} \right ]^2}\]
onde \(\mathit{N_{(h)}}\) é o número de pares de valores medidos \(\mathit{Z_{(x_{i})}}\), Z(xi+h) \(\mathit{Z_{(x_{i}+h)}}\), separado pela distância \(\mathit{h}\). Se a variável for escalar O gráfico de \(\mathit{\gamma _{(h)}}\) versus os valores correspondentes de \(\mathit{h}\) é chamado de semivariograma.
Para estimar valores de precipitação em locais não amostrados e confeccionar os mapas, utilizou-se do método de interpolação krigagem ordinária (KRIG), que é um método geoestatístico, que se utiliza de um estimador linear não-viciado com mínima variância e leva em consideração a estrutura de variabilidade encontrada para o atributo e é definido pela seguinte equação:
\[\mathit{Z_{(X_{i},X_{i} + h)}^{*} = \sum_{i=1}^{n}\lambda_{i}Z_{(X_{i},X_{i} + h)}}\]
sujeito a:
\[\mathit{\sum_{i=1}^{n}\lambda_{i} = 1}\]
em que \(\mathit{Z_{(X_{i},X_{i} + h)}^{*}}\) é o estimador para um ponto \(\mathit{(X_{i},X_{i} + h)}\) da região e \(\mathit{\lambda_{i}}\) são os pesos usados na estimativa.
11 Modelgem do Semivariograma dos Dados TRMM
11.1 Análise Descritiva
Ao se utilizar os método de Krigagem é recomendável que os dados se assemelhe com uma distribuição normal. Caso não seja, deve-se realizar uma transformação afim de que possa torná-la normal. Igualmente importante, que seja uma transformação em que pode-se realizar o procedimento inverso, de volta para os dados originais.
Verificação de Normalidade dos Dados TRMM.
Análise descritiva dos Dados TRMM.
Como mostrado na figura acima e, para maior comprovação, foi submetido ao teste de normalidade de shapiro Wilk, que mostrou o \(p-value < 2.2*10^{-16}\). Portanto muito menor que p-value igual 0,05(Para p-values maior ou igual a 0,05, distribuição normal).
11.2 Mapa Variográfico
Seguindo na investigação dos dados TRMM, também é necessário verificar se os dados posui uma tentência em torno das coordenadas de longitude e latitute.
Os efeitos direcionais (tendência e anisotropia) são visualizados na forma de Mapa Variográfico. E caso seja presente, é necessário remover esses efeitos para seguir na análise variográfica.
Análise de Efeitos Direcionais dos Dados TRMM.
Verifica-se que as linhas que representam variogramas para diferentes direções não se cruzam diversas vezes, havendo clara distinção de tendências. Deste modo, consideramos falsa a premissa de isotropia.
Análise de Efeitos Direcionais dos Dados TRMM.
11.3 Semivariograma
Os dados foram interpolados para uma grade de pixels de tamanho 400x400 m, por krigagem ordinária. Utilizou-se o estimador clássico de Matheron (1963) para o cálculo das semivariâncias, e o alcance foi limitado à metade da distância máxima entre dois pontos (cutoff = 50% = 1416783m ) e os lags elegido como sendo a distancia minima entre os pares de pontos (width = 27333,34m). O ajuste do semivariograma experimental foi realizado pelo melhor modelo identificado na análise geoestatística e os parâmetros obtidos foram utilizados no processo de interpolação por krigagem. Após a interpolação ser processada foi realizada a criação e geração dos mapas temáticos para cada modelo teórico (Exponencial,Esférico e Gaussiano)
Variograma Experimental.
12 Resultados Obtidos Pelo Método de Krigagem Ordinaria
12.1 Modelo Exponencial
Modelo Exponencial.
12.2 Visualização do Modelo Exponencial com Open Street Map
12.3 Modelo Esférico
Modelo Esférico.
12.4 Visualização do Modelo Esférico com Open Street Map
12.5 Modelo Gaussiano
Modelo Gaussiano.
12.6 Visualização do Modelo Gaussiano com Open Street Map
13 Conclusão
Os modelos geoestatísticos, Exponencial e Esférico, apresentaram menor raiz do erro médio quadrático. E o modelo determinístico, IDW, apresentou o maior RMSE.
conforme tabela abaixo:
| Modelos | R-Quadrado | Raiz do Erro Médio Quadrático |
|---|---|---|
| Exponencial | 0.9922804 | 0.0879277 |
| Gaussiano | 0.9732477 | 0.1635797 |
| Esférico | 0.9921694 | 0.0885192 |
| IDW | 0.9020000 | 4.5300000 |