quarta-feira, 6 de junho de 2007

Métodos de interpolação

Introdução

Este trabalho tem por objectivo o desenvolvimento de uma carta de precipitação de uma região do Sahel no Ocidental africano (região identificada na figura 1), a partir de um mapa vectorial apenas com os valores de precipitação para algumas estações nesta mesma zona.

Figura 1 – Mapa do Sahel [1]

Em Sistemas de Informação Geográfica, para atingir este objectivo utiliza-se uma metodologia, denominada interpolação. Esta metodologia consiste numa estimativa dos valores do atributo dos locais não amostrados, partindo de outros pontos dos quais se tem valores obtidos no terreno, na mesma zona para a qual pretendemos fazer a estimativa. A base da interpolação é a tendência que os valores de um atributo têm para ser mais semelhantes em zonas próximas e mais diferentes em zonas afastadas.
É então possível passar um conjunto de dados de amostragens ou observações, no formato vectorial, para um mapa contínuo, no formato raster (que é o formato mais utilizado quando se pretende comparar mapas, fazer operações de álgebra de mapas ou adicionar vários layers à mesma imagem). O mapa obtido tem várias utilidades, entre elas e provavelmente a mais importante, a sua utilização como base para apoiar processos de decisão espacial.
O conjunto de dados a analisar neste trabalho estão organizados na figura 2, apresentando por cada estação o valor da precipitação que foi medido nesta.

Figura 2 – Dados recolhidos no terreno

Analisando esta figura é notória uma tendência para que a quantidade de precipitação diminua de cima para baixo, o que revela que a precipitação (variável em análise) é direccional.

Para uma melhor análise a nível quantitativo do conjunto de dados recolhidos fez-se um histograma destes dados, que é apresentado na figura 3.

Figura 3 – Histograma dos dados recolhidos

Este histograma foi construído com base em 252 amostras e cada classe tem uma amplitude de 5mm, a média é de 148,345mm, o desvio padrão 42,055mm, o valor máximo de precipitação 228mm e o mínimo zero.
Analisando a forma do histograma pode-se considerar que é na zona entre os 160 e 180mm que se verifica maiores frequências, isto é, a quantidade de dados obtidos neste intervalo é superior à quantidade de dados obtidos noutras classes.

Com esta ideia aproximada da distribuição dos dados tanto no espaço como a nível estatístico pode-se proceder à interpolação, tendo sempre em conta que esta não dará uma superfície que exprima a realidade tal qual como ela é, mas sim uma estatística e uma aproximação dos valores reais. Esta margem de erro que existe na interpolação pode ser avaliada de várias formas, através de processos de validação por exemplo, e de uma forma menos rigorosa, através da comparação dos histogramas dos mapas obtidos entre si e com o histograma dos dados.


Métodos que podem ser utilizados na interpolação

Existem vários métodos que podem ser utilizados para interpolar uma imagem, sendo que cada um deles tem as suas vantagens e desvantagens e estão melhor ou pior adaptados a um determinado tipo de dados.

A primeira forma de interpolação a ser utilizada é a análise de tendências.
Este método é aproximado, global e estatístico, isto é, as tendências gerais dos dados são respeitadas e traduzidas por uma função para toda a área em estudo, que tem em conta a aleatoriedade da variável e avalia a superfície criada estatisticamente.
A análise de tendências consiste na adaptação de um polinómio aos dados que foram recolhidos, através de uma regressão múltipla dos valores do atributo em função da localização geográfica.
O método apresenta algumas vantagens como por exemplo os erros aleatórios encontrados em cada ponto de recolha de dados ser independente e o ajustamento ser dado por dados estatísticos, no entanto, só costuma produzir bons resultados quando a variável em estudo apresenta uma relação conhecida com outras variáveis e é muito sensível à existência de outliers.
Existem outras limitações inerentes à análise de tendências, como por exemplo o facto de um modelo polinomial produzir uma superfície arredondada, que depois não se verifica na realidade.

No caso em estudo realizaram-se análises de tendências com polinómios de graus diferentes, que geraram superfícies com variações ligeiramente diferentes como se pode observar nas figuras 4,5 e 6. Para obter estas figuras e o histograma no IDRISI utilizaram-se os seguintes comandos:

Figura 4 – Análise de tendências ajustada com polinómio de grau 1


Figura 5 – Análise de tendências ajustada com polinómio de grau 2


Figura 6 – Análise de tendências ajustada com polinómio de grau 3

Obtidas as superfícies, pode-se verificar que todas elas apresentam o mesmo comportamento direccional dos dados, mas que na superfície de grau 3 esse padrão já é um pouco mais complexo e não é tão contínuo como nas outras duas superfícies.


Outro método que pode ser utilizado é o inverso da distância pesada (IDW), em que o valor do ponto que se está a estimar é obtido fazendo uma média ponderada com os pontos mais próximos. Desta forma, quanto mais próximo estiver um ponto do ponto que se está a estimar maior será o seu peso e quanto mais afastado estiver, menor será o seu peso.
O IDW é um método local, pois aplica-se repetidamente a diversos subconjuntos de pontos do total de pontos amostrados; é um método determinístico uma vez que não utiliza nem a estatística nem a probabilidade; e finalmente é um método exacto visto que todos os pontos obtidos através de amostragem são respeitados e a superfície obtida passa por todos os pontos conhecidos.
Uma particularidade deste método é ser possível definir qual o peso da distância no cálculo da média ponderada através do expoente que é aplicado à distância (quanto maior for o expoente, maior é a influência da distância na média).

Para se implementar este método no IDRISI e o respectivo histograma utilizaram-se os seguintes comandos:

Neste caso específico foram calculadas as superfícies através deste método com a utilização de expoente 1 e expoente 2, como se pode verificar nas figuras 7 e 8.

Figura 7 – Inverso do peso da distância, expoente 1

Figura 8 – Inverso do peso da distância, expoente 2

Estas duas figuras apresentam algumas variações entre si, que no entanto não são muito acentuadas e como os mapas obtidos através da análise de tendências, apresentam o mesmo comportamento direccional revelado pelos dados, embora com formas mais acentuadas e menos suaves que as obtidas por análise de tendências.

Pode também ser utilizado como método de interpolação o chamado método dos polígonos de Thiessen. Este interpolador é local, determinístico e exacto, como o IDW e pelas mesmas razões.
Este método normalmente ajusta melhor dados nominais do que dados contínuos, como é o caso, uma vez que as alterações nas fronteiras dos polígonos são abruptas e todo o polígono fica com o valor que estava atribuído ao ponto que foi utilizado para calcular o polígono.
O cálculo dos polígonos depende da configuração dos pontos amostrais e é determinada a região mais próxima dum ponto, segundo um polígono de Thiessen ou um poliedro de Voronoi (como também são chamados).
Este método é robusto, mas não é sensível ao tipo de variável que está em estudo.
Para se utilizar este método no IDRISI é necessário seguir os seguintes passos:

Gerou-se uma superfície usando como interpolador os Polígonos de Thiessen que se encontra representada na figura 9.

Figura 9 – Polígonos de Thiessen

Como se pode observar na figura, o comportamento direccional mantém-se. Contudo para uma variável como a precipitação, que apresenta um comportamento muito mais continuo do que discreto, parece que os polígonos de Thiessen fornecem informação pouco detalhada e exacta, pois áreas significativas tomam um único valor de precipitação, o que não se adequa à realidade.

Existe ainda outro interpolador bastante utilizado denominado kriging. Este faz parte da modelação geoestatística que tem como base a premissa que pontos mais próximos no espaço tendem a ter valores mais parecidos que pontos mais afastados.
Este método assume a variabilidade aleatória que existe entre os vários pontos amostrais, estudando variáveis regionalizadas, isto é, variáveis que apresentam continuidade ponto a ponto mas que as mudanças entre eles são tão complexas que não pode ser descrita por nenhuma função determinista.
Para se verificar se a variável é regionalizada faz-se um variograma, a partir da semivariância. O variograma indica o nugget (quando a semivariância é muito reduzida), o sill (patamar a que a semivariância estabiliza) e o range (distância a partir da qual se atinge o sill e já não se pode considerar o modelo válido). O variograma também pode indicar se a variável é anisotrópica (quando a variável tem comportamentos distintos segundo diferentes direcções).
O kriging baseia-se numa média ponderada que atribui pesos aos pontos conhecidos, através dos valores de semivariância que estes apresentam.
Este interpolador tem várias vantagens porque os erros de estimação apresentam a variância mínima, é linear (as suas estimativas são combinações lineares dos dados conhecidos) e não enviezada (tenta que a média dos erros seja nula).

Para avaliar se a variável em estudo é regionalizada ou não fez-se o variograma da precipitação na região do Sahel. Para construir o variograma ajustado aos dados da precipitação seguiram-se os seguintes comandos no IDRISI:

O variograma da precipitação na região em estudo é apresentado na figura 10, juntamente com a nuvem dos resíduos.

Figura 10 – Variograma e nuvem dos resíduos

Pela nuvem dos resíduos pode-se constatar que os dados apresentam realmente um comportamento direccional logo o variograma que vai ser calculado deve ter em conta esse aspecto.
Do variograma obtido apenas interessam os pontos e devem ser retirados 3 valores que serão essenciais para o ajuste de um modelo ao variograma. Esses 3 valores correspondem ao nugget (ou efeito pepita) que é de 270, ao sill que toma o valor de 550 e ao range que é de 200. Estes são os valores que serão utilizados para ajustar um modelo aos pontos. A figura 11 traduz esse ajuste.

Figura 11 – Variograma ajustado

O modelo que foi utilizado foi o Gaussiano e como se pode verificar este modelo apenas ajusta o conjunto dos dados até ao valor do range, logo para pontos a distâncias superiores a 200 este modelo não é efectivo e não tem um bom ajuste.

Depois de se confirmar que a variável é regionalizada e anisotrópica, procede-se com o kriging, seguinto os passos abaixo descritos no IDRISI e deste procedimento resultam dois mapas.

O primeiro mapa (figura 12) representa a distribuição da precipitação na área pretendida, o segundo (figura 13) representa o erro associado a cada zona do mapa de precipitação, portanto as zonas em que o modelo gerado através de kriging está mais ajustado.

Figura 12 – Kriging


Figura 13 – Mapa dos erros associados ao kriging

Analisando o mapa de precipitação mantém-se o comportamento direccional dos dados, como seria de esperar e a superfície é suave. O mapa dos erros permite verificar que no canto superior esquerdo o erro é muito elevado, logo aquela zona do mapa não estará muito ajustada. Também é possível verificar que as zonas com maior erro são zonas em que a densidade de pontos com valores conhecidos é menor, logo existe menos informação por onde estimar a superfície de distribuição da precipitação.


Análise comparativa dos resultados de cada método

Tendo obtido quatro mapas de interpolação diferentes para os mesmos dados, é necessário agora compará-los, bem como aos histogramas resultantes de cada mapa, com os dados recolhidos para determinar a qualidade de cada método neste caso.


O primeiro método foi a análise de tendências, em que se obterem três mapas distintos ajustados com grau 1, 2 e 3 do polinómio.
No caso da análise de tendências ajustada com polinómio de grau pode ver-se que os valores (Fig.4) seguem a mesma distribuição direccional dos dados recolhidos (Fig.5), aumentando da parte inferior do mapa para o cimo. Analisando os dois histogramas (Figs.3 e 14) verificamos que aqui existem muitas discrepâncias, nos dados recolhidos ocorrem variações acentuadas em relação á quantidade de valores que se tem para cada nível de precipitação, sendo os valores entre 160 e 180mm aqueles em que a frequência de dados é maior; no caso do histograma da analise de tendências com polinómio de grau 1 apresenta mais ou menos a mesma frequência para todos os valores de precipitação entre 43 e 203mm aproximadamente, sendo esta progressivamente menor para precipitações acima e abaixo destes valores. Isto seria de esperar tendo em conta que no mapa os valores se distribuem de forma quase regular, para os valores intermédios de precipitação, sendo as áreas de valores extremos mais reduzidas.

Figura 14 – Histograma da análise de tendências ajustada com polinómio de grau 1


No mapa da analise de tendência com polinómio de grau 2 (Fig.5) tem-se uma distribuição um pouco mais distinta da do mapa dos dados (Fig.2), sendo que aos valores de precipitação em cada um dos mapas não correspondem bem à mesma área, porém mantêm a mesma direcção,
aumento dos valores de baixo para cima. Na análise dos histogramas (Figs.3 e 15) podemos ver que embora em ambos as frequências dos dados estejam a aumentar, os valores que não correspondem um ao outro, neste caso os valores de precipitação com maior frequência estão entre 170 e 190mm, ou seja um pouco superiores aos valores medidos.

Figura 15 – Histograma da análise de tendências ajustada com polinómio de grau 2

Por ultimo na analise de tendências com ajustamento a polinómio de grau 3 temos um mapa (Fig.6) com um padrão mais complexo do que os anteriores, mas que segue a mesma direcção porém de forma menos continua; este é dos três mapas aquele cuja distribuição mais se assemelha com a do mapa de dados (Fig.2), pois este também não tem uma distribuição tão contínua nem uniforme como os outros mapas fazem crer. Nos histogramas as semelhanças já não são assim tantas, enquanto que nos dados a frequência dos valores vai aumentando á medida que a precipitação aumenta, neste (Fig.16) existem dois picos de frequência, um entre os 144 e 164mm aproximadamente e o outro entre os 164 e os 184mm aproximadamente; é de notar que este ultimo corresponde a valores de precipitação muito próximos daqueles em que a frequência é maior no dados (160 e 180mm).

Figura 16 – Histograma da análise de tendências ajustada com polinómio de grau 3

Com o método do inverso da distância pesada obtiveram-se dois mapas sem diferenças muito significativas, um para exponente 1 (Fig.7) e outro para exponente 2 (Fig.8), com superfícies que seguem o mesmo comportamento direccional dos dados, embora de forma mais acentuada e menos suave. Os padrões são muito semelhantes aos das análises de tendências. Em comparação com os dados recolhidos o seu padrão é mesmo muito idêntico, revelando até valores de precipitação isolados no meio de outros diferentes.

Os seus histogramas (Figs.17 e 18) também são bastante semelhantes entre si, e com o histograma dos dados recolhidos (Fig.3). Em todos eles a maior frequência de dados obtidos está entre valores de precipitação de 160 e 180mm, sendo menores para valores mais baixos e mais altos, mas com algumas oscilações.

Figura 17 – Histograma do inverso do peso da distância com exponente 1


Figura 18 – Histograma do inverso do peso da distância com exponente 2

Outro método usado foi o dos polígonos de Thiessen, com o qual obtemos um mapa (Fig.9) com uma superfície bastante descontínua e discreta que atribui valores todos iguais a uma área á volta dum ponto com valor conhecido. Porem comparando-o com o mapa dos dados recolhidos é bem visível que segue a mesma direcção e um padrão algo semelhante.

Comparando o histograma referente a este método (Fig.19) com o dos dados recolhidos, as semelhanças são muito poucas, desta forma o mais sensato será assumi-las insignificantes.

Figura 19 – Histograma dos polígonos de Thiessen

O último método utilizado foi o Kriging, que fornece a informação em dois mapas, um do modelo (Fig.12) e outro dos erros a ele associados (Fig.13). Porém para a análise comparativa apenas nos interessa o modelo, que é dado como uma superfície suave, bastante semelhante ao modelo de análise de tendências, mas neste caso descontínua. Comparando-o com o mapa dos valores recolhido verifica-se que segue o mesmo comportamento direccional e tem uma distribuição muito semelhante.

Olhando para os histogramas, o dos dados recolhidos (Fig.3) e o do modelo (Fig.20), verifica-se que em termos de aspecto têm algumas semelhanças, contudo com uma análise pormenorizada verifica-se que isto não é verdade. No modelo a maior frequência dos dados ocorre para valores de precipitação da ordem dos 170 / 190mm, ou seja um pouco superiores aos valores em que as frequências são maiores nos dados recolhidos.

Figura 20 – Histograma do Kriging

Para validar a interpolação recorre-se a um conjunto de pontos de coordenadas e valores de precipitação conhecidos no mapa dos dados recolhidos e, a partir deles determinamos o seu valor de precipitação em cada mapa de cada método. Através desses valores criamos uma matriz de validação feira (Fig.21) que nos permite analisar os erros de cada modelo.

Tabela 1 – Matriz de validação

Pela análise desta tabela verifica-se que o método que apresenta maior erro médio, 643.666667, é o método dos polígonos de Thiessen, isto deve-se ao facto de se atribuírem a uma determinada área, em redor dum ponto de valor conhecido, o mesmo valor desse ponto. Assim ao verificar um valor de um ponto, por este modele, ele será igual ao valor do ponto que lhe estiver mais próximo, e não tem em conta que possa ser ligeiramente diferente.

Por outro lado o método do inverso do peso da distância de exponente 2 é aquele em que o erro médio é menor, 106.333333, isto porque este método atribui a um ponto cujo valor não se conhece um valor médio ponderado de todos os pontos á sua volta, desta forma continua a persistir a maior influencia dos pontos mais próximos, mas o ponto que se pretende determinar apresentará um valor diferente deles, o que minimiza as discrepâncias entre os valores estimados e os valores medidos.


Escolha do modelo mais apropriado

Perante todos estes modelos é necessários escolher aquele, que pelas suas características, melhor se adequa ao nosso estudo. Assim é necessário que o modelo cumpra com alguns requisitos.

Sendo a precipitação uma variável com uma distribuição relativamente continua é de esperar que o mapa de interpolação, feito através dos dados das estações de monitorização, também o seja. É também de esperar que este apresente um padrão de distribuição e uma direcção semelhante aos dados, neste caso valores de precipitação que em geral aumentam de cima para baixo.

Outro aspecto importante a considerar, agora em temos de histograma, é frequência dos valores, que devem ser semelhantes, embora em ordens de grandeza diferentes, tanto no histograma dos dados recolhidos como no histograma do modelo de interpolação.

Por fim é necessário verificar o valor dos erros médios, que devem ser o mais pequeno possível, isto para que os valores dos pontos no mapa interpolado sejam o mais parecido possível com os valores recolhidos nesses pontos.

Desta forma, embora muitos possam satisfazer uma ou mais condições, apenas num dos modelos podem verificar-se todas estas condições. Reconhece-se assim como mais apropriado o método do inverso da distância pesada com exponente 2, uma vez que nele se verificam todas as condições necessárias para se considerar um bom modelo de interpolação.

Todos os outros métodos apresentam uma ou outra condicionante que não lhe confere uma boa adaptação a este caso. A análise de tendências dá uma superfície muito suave e não tem em conta ou outliers, os polígonos de Thiessen não são apropriados para variáveis contínuas e por isso levavam a um erro médio bastante elevado, e o Kriging gera erros associados na zonas em que a densidade de pontos é menor.