A saga da água

Cai de pé e corre no chão

Tudo indica que em Titã, a lua de Saturno, chove metano líquido. O metano escorre em rios e forma mares de hidrocarbonetos. Não temos a sorte de ter aqui na Terra lagoas de gasolina. Mas temos água.

Aqui, a água evapora dos oceanos e é trazida por ventos para os continentes. Como que tudo que sobe tem que descer, chove. A água volta ao estado líquido e é drenada pela ação da gravidade de volta aos oceanos. O detalhe é que a água ao cair encontra a superfície da crosta terreste, que é um lugar bem esquisito.

Em contato com a crosta, a água da chuva passa a integrar o longo e contínuo processo de dissecação do terreno, de erosão das rochas e transporte de sedimentos para os oceanos. Na prática, isso quer dizer que a chuva forma bacias de drenagem, ou bacias hidrográficas. E um problema muito prático consiste em descobrir a área de uma bacia a montante de um ponto de saída, ou exutório.

Nesse post, vamos explorar formas de analisar essa questão, usando o QGIS e modelos digitais de elevação (MDE).

O divisor de águas

Se o seu terreno for bastante acidentado, é possível interpretar o padrão de drenagem por uma simples inspeção visual. E aqui entra em cena novamente o Modelo Digital de Elevação (MDE) como dado fundamental para a análise.

Olhe essa imagem de satélite, que eu importei como mapa de fundo no QGIS:

O lugar fica no município de Caará, no Rio Grande do Sul. É um lugar com o terreno bem acidentado, onde existem várias cachoeiras de águas limpas, não muito distante de Porto Alegre. Mas a vegetação e ausência de sombreamento na imagem não permitem que consigamos interpretar muito bem o padrão de drenagem.

Agora, no QGIS,  inseri uma camada de sombreamento semi-transparente. Ficou melhor:

Mas dá para melhorar mais. Removi a imagem de satélite e deixei apenas um MDE estilizado e o sombreamento. Aqui, já é possível ver o padrão de drenagem:

Mas considerando um ponto exutório, onde exatamente é o limite da bacia de drenagem para esse ponto?

A linha que limita a bacia é o divisor de águas, um caminho no terreno que funciona como uma cumeeira de um telhado. Da água que chove, parte escorre para um lado, parte para o outro.

Uma maneira de interpretar melhor o MDE é gerar curvas de níveis vetoriais (shapefile). Para isso, Vá em GDAL > Extraction > Contour e execute o algoritmo informando o MDE e ajustando os parâmetros opcionais. Eu gerei aqui as curvas de nível para intervalos de 10 metros:

Exatamente como o cume de um telhado, o divisor de águas é bem mais visível com o auxílio das curvas de nível:

O bom e velho método expedito

Visualizar o padrão de drenagem e identificar onde passa o divisor de águas não responde diretamente a nossa questão. Queremos saber a área, em forma e quantidade, da bacia de drenagem.

Com as curvas de nível, podemos resolver esse problema de forma muito prática no QGIS. A partir do ponto exutório, desenhe manualmente a área de bacia, seja pela ferramenta de mensuração de área, seja ao editar os vértices do polígono de um shapefile criado para esse propósito. Só precisa ter o cuidado de manter o divisor de águas sempre ortogonal às curvas de nível. Esse é o método expedito.

Para ser rápido, o método expedito precisa ser ruim. Para ser bom, você precisaria editar um polígono com vértices muito próximos uns dos outros, um processo muito lento para ser feito manualmente.

Para se ter uma ideia do problema, o MDE derivado da missão SRTM tem píxeis de 30m por 30m, o que implica que a resolução espacial desse dado é de 30m. Logo, se você quiser obter um polígono de bacia de drenagem com essa resolução, os vértices desse polígono precisam ter a densidade de pontos semelhante: a cada 30m, um vértice. Em bacias com centenas de quilômetros quadrados, a menos que você fique uma semana inteira clicando e criando vértices, esse método só irá produzir aproximações grosseiras, de resolução muito inferior ao dado original disponível no MDE.

Entrando na Matrix

Para acessar a resolução espacial disponível no MDE, precisamos entrar na matrix.

Um MDE nada mais é que uma matriz de números que representam a altitude de uma amostra da superfície. Essa amostra é o píxel na imagem que vemos, que possui uma resolução espacial. Como já falamos, para os dados da missão STRM, o píxel é uma amostra aproximadamente quadrada de 30m por 30m.

Sendo uma matriz, podemos aplicar diversas técnicas de computação para identificar o gradiente do terreno, a direção de fluxo e para onde o fluxo aponta. O processo é análogo ao cálculo diferencial de duas variáveis, com x e y sendo a longitude e a latitude, e z sendo a altitude do terreno. A vantagem é que o QGIS faz as contas para nós.

A análise de bacia hidrográfica no QGIS pode ser feito de várias formas, usando pacotes de algoritmos diferentes, como o SAGA e o GRASS. Aqui vamos explorar o SAGA. Como reza a cartilha, vamos antes de tudo avaliar o fluxograma do processo como um todo.

Queremos obter:

1. a área de drenagem associada a um ponto exutório, e;
2. o polígono da bacia de drenagem associada a esse ponto exutório.

O fluxograma que nos levará para esses resultados, explicado a seguir, é o seguinte:

Passo 0: entre na Matrix

Obviamente, você precisa de um MDE para fazer essa análise. Eu fiz um post demonstrando como produzir um MDE a partir dos dados brutos da missão SRTM. Clique aqui para acessar o post. Siga os passos até o final, reprojetando o MDE para coordenadas planas (a não ser que seu MDE seja gigante, de abrangência nacional ou algo do gênero).

Passo 1: tapar buracos

A primeira etapa consiste em produzir um MDE sem depressões no terreno. Os “buracos” no terreno podem ser resultado direto da mensuração do sensor que produziu o dado ou produzidas durante os processos de confecção do mosaico e reprojeção do MDE. Na primeira alternativa, podem ser depressões legítimas do terreno ou depressões ilegítimas, como em lugares muito planos onde até mesmo talhões de eucalipto são detectados como pequenas elevações no terreno. Ainda também existe a possibilidade de serem resultados aleatórios do processo de coleta do dado, um ruído.

Independentemente da sua origem, não queremos as depressões! O motivo é o seguinte: os algoritmos que vamos usar “navegam” na matriz do MDE, seguindo a linha de fluxo morro abaixo ou morro acima. Se existe uma depressão no caminho, os algoritmos entendem que é o fim da linha e encerram suas contas.

Entre as alternativas disponíveis, recomendo o SAGA Fill Sinks. Também recomendo ajustar os parâmetros avançados para a interpolação ser bilinear –  vai faciltar a vida do seu computador. O algoritmo está em Processing > SAGA > Terrain Analysis – Hydrology > Fill Sinks.

Se você não acredita que existem depressões no terreno, observe as curvas de nível no seu MDE original e depois compare com as curvas de nível do MDE gerado pelo SAGA Fill Sinks. Se você encontrar curvas fechadas em lugares como vales ou planícies, ali existe uma depressão. No meu MDE, encontrei várias, de vários tamanhos.

Passo 2: morro abaixo

Agora que temos um MDE liso, sem depressões, podemos navegar na matriz do MDE e seguir as linhas de fluxo, dos lugares mais altos em direção aos lugares mais baixos. Morro abaixo. Precisamos fazer isso para descobrir exatamente onde a água se acumula e forma riachos e rios em nosso MDE.

Sabendo exatamente onde a água passa, podemos obter, com a ferramenta de identificação do QGIS, as coordenadas X e Y exatas do píxel do exutório que queremos avaliar a área de drenagem.

Entre algumas formas de fazer isso, a mais prática é o algoritmo SAGA Catchment Area. O algoritmo está em Processing > SAGA > Terrain Analysis – Hydrology > Catchment Area. Cuidado, pois você precisa informar para o algoritmo o MDE sem depressões, e não o original! Dos parâmetros avançados, eu costumo solicitar interpolação bilinear, para dar uma trégua ao meu computador.

O resultado é um raster em que cada píxel tem um valor de área, em metros quadrados (se você estiver com um MDE não projetado em coordenadas planas, a área informada será em graus quadrados – e boa sorte na conversão de unidades). Portanto, resolvemos dois problemas de uma só vez: não só conseguimos extrair as coordenadas X e Y do píxel do exutório, mas também o valor do píxel em si já nos diz a área acumulada, em metros quadrados.

Com a ferramenta de indentificação (o cursor com um “izinho”) e a camada raster produzida pelo SAGA Catchment Area selecionada, clique no píxel que você quer que seja o exutório da bacia e avalie o valor de área de drenagem e as coordenadas X e Y. No meu caso aqui a área de drenagem deu 452 hectares. Copie as coordenadas em algum lugar, pois vamos usar no próximo passo.

Dica: se você obter um raster completamente preto, diferente do que eu apresento aqui em baixo, antes de assumir que algo deu errado confira a estilização do raster, alterando os limiares de máximo e mínimo.

Passo 3: morro acima

Obtendo as coordenadas X e Y do píxel que queremos avaliar a bacia de drenagem, o ponto exutório, podemos agora subir o morro. Isso quer dizer que vamos informar as coordenadas e um algoritmo vai navegar na matriz do MDE sem depressões. Enquanto o algoritmo estiver apenas subindo em termos de altitude, quer dizer que ele está dentro da bacia de drenagem. Assim que ele começar a descer, o algoritmo vai saber que saiu da bacia – cruzou o divisor de águas.

Com o SAGA, isso é feito pelo algortimo SAGA Upslope Area, também localizado em Processing > SAGA > Terrain Analysis – Hydrology. Não esqueça que o MDE a ser informado é o MDE sem depressões. Recomendo manter os parâmetros padrão, com exceção para a interpolação, que eu costumo mudar para bilinear.

O resultado é um raster pseudo-booleano. O valor de píxel é 100 para os píxeis que estão dentro da bacia e 0 para os que estão fora.

Passo 4: saindo da Matrix

Existem várias formas de converter o resultado do SAGA Upslope Area em um shapefile da bacia de drenagem. Uma delas é mantendo a fidelidade ao SAGA. Existe um algoritmo nesse pacote que converte classes específicas de um raster em shapefile. O algoritmo é o Vectorising Grid Classes, localizado em Processing > SAGA > Vector <-> raster. No formulário de parâmetros, especificique que você quer conveter a classe de valor de píxel igual a 100.

O resultado é um shapefile com o formato exato da bacia de drenagem de acordo com nosso MDE. As bordas do shapefile são serrilhadas porque o algortimo não faz suavizações dos píxeis. Você pode ficar tentado a querer suavizar a borda, aplicando outras ferramentas, mas isso irá apenas degradar a qualidade da resolução do nosso dado. Não é bonito, mas é honesto. E aqui estamos com nosso resultado final.

Bônus: automatizando (quase) tudo

A partir de um MDE já projetado, não seria possível automatizar todo o fluxograma?

Sim, é possível. Mas exite uma parte muito sensível ao processo: quando o usuário precisa inspecionar o raster de “Catchment Area” e obtêm ali as coordenadas X e Y do píxel que se deseja analisar como ponto exutório. Dá para automatizar esse procedimento, pedindo para o usuário um shapefile de ponto e fazendo as contas para aproximar o píxel exutório a partir desse ponto. Mas essa parte eu vou deixar para um outro post.

Sendo assim, vou aqui propor dois procedimentos automatizados, segmentados por um procedimento humano.

O primeiro procedimento automático é a geração do MDE sem depressões e a consequente geração do raster de área acumulada (Catchment Area). Usando o Modelador de processos, o processo ficou assim:

Clique aqui para baixar a primeira parte do modelo para o QGIS 2.

Após obter manualmente as coordenadas do ponto de interesse, começa o segundo procedimento automático. Nesse procedimento, nós podemos ignorar o raster produzido pelo SAGA Upslope Area e focar somente no resultado final, o shapefile da bacia.

Clique aqui para baixar a segunda parte do modelo para o QGIS 2.