Guia para o cálculo da Distância Vertical para a Rede de Drenagem (Vertical Distance to Channel Network, VDCN), tutorial passo-a-passo usando apenas os softwares de SIG de código-livre QGIS e SAGA GIS

Colegas pesquisadores e entusiastas do GIS de código-livre,

Primeiramente, bem-vindos ao meu blog!

Gostaria de começar com um aviso - posso ser uma pesquisadora desta área, mas isso não quer dizer que tudo o que faço ou escrevo aqui funcionará para você, em suas próprias configurações de desktop e versões de packages. Não assumo responsabilidade nenhuma se você perder dados ou bagunçar sua instalação. Eu também não autorizo nenhum tipo de cópia do meu conteúdo.

Hoje, o post é sobre como gerar o raster para a Distância Vertical para a Rede de Drenagem (Vertical Distance to Channel Network ou VDCN) usando QGIS Raster Calculator, e o SAGA cmd (Command Line Interface – CLI) ou o SAGA no QGIS, passo a passo.

A VDCN é a distância vertical até a rede de drenagem interpolada ou até a rede fluvial ou rede de drenagem do local.

Existe uma ferramenta no SAGA chamada “Vertical Distance to Channel Network” que pode ser usada para isso. Porém, o post atual é sobre a possibilidade de você mesmo realizar os cálculos e ter mais controle sobre esse processo. Ou, em alguns casos, a ferramenta SAGA leva muito tempo para big data espacial, então é melhor dividir o processamento em pequenas etapas.

Por que calcular a distância vertical até a rede de drenagem (Vertical Distance to Channel Network)?

O VDCN é um importante atributo do terreno que fornece informações sobre a profundidade aproximada do fluxo de base de uma bacia hidrográfica em um ponto no espaço. Isso é útil para saber a que distância vertical a superfície do terreno está do lençol freático, o que tem utilidade para várias aplicações em estudos de Hidrologia, Geomorfologia e Desastres Naturais.

O nível do fluxo de base é a elevação (absoluta) na que a água subterrânea do lençol freático pode ser encontrada. A menos que esse valor seja medido em campo, o nível de base tem de ser estimado por medições indiretas. Se você tiver medições de campo da profundidade do lençol freático em seu local, consistentes e espacialmente distribuídas de maneira correta, você pode simplesmente interpolar suas medições para um raster e esta é a distância vertical para a rede de drenagem.

Se você não tem essa medição in situ, você pode seguir as etapas abaixo, ou empregar outra ferramenta de sua escolha, à exemplo da “Vertical Distance to Channel Network” do SAGA GIS, que exige a rede de drenagem como uma das entradas.

Para este tutorial, você precisará:

• Seu Modelo Digital de Elevação (MDE) / Modelo Digital de Terreno (MDT)
• QGIS 3 ou superior (faça download em http://qgis.org)
• SAGA GIS incluído no QGIS ou a versão de Linha de Comando saga_cmd.exe (Baixe em http://saga-gis.org - Eu ensinei como usar neste post)

Bem, vamos lá. Pronto para começar?

1. (opcional) Reprojetar o MDE (raster)

Reprojete o raster para coordenadas projetadas (como por exemplo UTM), em vez de geográficas (latitude e longitude). Meu algoritmo favorito para fazer isso é o gdal Warp. Ele pode ser encontrado no QGIS 3 na caixa de ferramentas de processamento (Processing Toolbox).

2. Preencher depressões em seu modelo digital de elevação (MDE)

Use sua ferramenta preferida para preencher as depressões em seu MDE. Minha preferida é a SAGA Sink Removal. No Prompt, use saga_cmd (verifique este post para obter detalhes sobre como usar saga_cmd):

saga_cmd ta_preprocessor 2 -DEM input_DEM.tif -DEM_PREPROC DEM_filledsinks.sdat

Se o seu raster for pesado, isso pode demorar algum tempo. No QGIS 3.16, essa ferramenta encontra-se na Processing Toolbox com o nome “Fill Sinks”.

3. Calcule o Fluxo Acumulado (Flow Accumulation) usando uma abordagem Determinística do tipo Deterministic 8 (D8)

Eu gosto mais de usar a abordagem determinística para esse tipo de algoritmo, porque ela gera um único canal em vez de vários canais paralelos.

3.1 O que uma acumulação de fluxo do tipo D8 faz, e por que é importante preencher os sumidouros/depressões antes de rodar o algoritmo para acumular o fluxo?

Vamos fingir que este é o seu MDE (os números são os valores da elevação do terreno, em metros)

Olhando para ele, é fácil perceber que o canal está na terceira linha.

Para realizar a acumulação de fluxo, primeiro, o algoritmo calcula a direção do fluxo, com base nas oito direções possíveis para cada ponto:

É por isso que é chamado de D8, abreviatura de Deterministic 8 (Determinístico 8).

A direção escolhida é aquela com a menor elevação, exceto se a elevação da célula em si for menor do que a de todas as células adjacentes. Para nosso MDE proposto, podemos ver que isso ocorre neste pixel:

Quando isso acontece, chamamos de depressão ou sumidouro (ou, em inglês, de sink). Quando rodamos um algoritmo de remoção de sumidouros, o que pretendemos remover é esse tipo de ocorrência, pois ela cria descontinuidades no acúmulo de fluxo.

Por quê?

Porque, com base nessas 8 direções de fluxo possíveis, uma é escolhida, e é como se toda a precipitação acumulada do pixel original fluísse para este outro pixel com a elevação mais baixa. Para calcular a acumulação de fluxo (Flow accumulation), primeiro, o número de células contribuintes (ascendente) é contabilizado. Para o MDE exemplificado:

O que podemos notar aqui?

O fluxo não está saindo desta região por causa da depressão. Preenchendo essa depressão, obtemos um novo MDE:

Que gera esse acúmulo de fluxo com base no algoritmo D8:

Muito melhor! Agora o fluxo está saindo desta parte do MDE, como deveria. Alguns algoritmos multiplicam o número de células pela área de cada célula para calcular a área de drenagem à montante de cada célula.

3.2 Calculando o Flow Accumulation

O módulo Flow Accumulation (Flow Tracing) do SAGA pode ser usado no saga_cmd da seguinte forma:

saga_cmd ta_hydrology 2 -ELEVATION DEM_filledsinks.tif -CAREA flow_accumulation.tif -CAREA_UNIT 0 -METHOD 0

CAREA_UNIT é a unidade usada pela ferramenta para acumular o fluxo. Escolha 0 para o número de células e 1, para a área de drenagem à montante. Se você escolher a opção 0, é fácil obter a área de drenagem: se você estiver trabalhando com uma projeção cartográfica projetada como a UTM, basta multiplicar o acúmulo de fluxo nas células pela área de cada célula. Digamos que o tamanho da célula seja 30 m, então a área de cada célula é 30 * 30 = 900 m². Por outro lado, se você escolher a opção 1, poderá obter o número de células à montante dividindo os valores pela área de cada célula. Ambas as opções são adequadas para o cálculo do VDCN.

Dependendo da versão do SAGA GIS, CAREA pode ser FLOW e CAREA_UNIT, FLOW_UNIT.

METHOD indica o método de execução dos cálculos. A opção 0 é o D8, que é a minha sugestão neste post. Mas, se quiser, escolha outras opções, você que manda no seu VDCN!

No QGIS 3.16, a ferramenta mais semelhante é a Catchment Area. Ela não especifica a unidade em que o fluxo será acumulado, mas após um pequeno teste, acredito se tratar da área (equivalente à opção 1 do saga_cmd).

4. (opcional) Calcule o logaritmo natural do raster do acúmulo de fluxo

A saída do Flow Accumulation é bastante difícil de ver no QGIS ou mesmo em outros softwares de geoprocessamento, porque a escala em que é representado o raster, por padrão é linear, mas a distribuição dos valores nesse raster se aproxima mais de um exponencial. Para ser capaz de observar as nuances de um raster de acumulação de fluxo, recomendo calcular o logaritmo natural dele no Raster Calculator do QGIS. Basta usar:

ln (“Seu Raster”) 

e salvar com um nome diferente.

E, em seguida, abra no seu software de geoprocessamento e observe.

5. Escolha o seu limiar ou ponto de corte

Parte da vantagem de fazer esse processo passo a passo é controlar as pequenas partes dele, como o limite para o qual uma célula pode ser considerada parte da rede de drenagem. Para isso, sugiro abrir o raster do Flow Accumulation ou o seu logaritmo no QGis e ir para Properties, Symbology. Escolha o mínimo e o máximo para representação na tela. Por exemplo, se você escolher o limiar entre 7 e 8, coloque 7 como o valor mínimo e 8 como o valor máximo. Vá observando o mapa e escolha um número apropriado para o qual a rede de drenagem e apenas a rede de drenagem é destacada.

Pense assim, esta é a sua rede de drenagem? Se parecer muito abrangente, escolha um limiar mais alto; se não mostrar as partes mais altas da bacia hidrográfica, escolha um limiar mais baixo. Algumas pessoas preferem fazer esta etapa com auxílio de uma foto de VANT ou imagem de satélite sobreposta.

Após escolher seu limite, guarde esse número para usá-lo na próxima etapa.

6. Reclassificar por tabela

Usando QGIS, use a ferramenta do Processing Toolbox chamada Reclassify by table para reclassificar seu raster usando uma tabela.

Pesquise na Processing Toolbox por Reclassify by Table. Clique no botão “…” no lado direito de Reclassification table para gerar sua tabela.

Clique em Add Row (adicionar linha)

Coloque o limite escolhido na etapa anterior como mínimo e um número muito grande (fora dos limites do seu raster) como máximo. Esses pixels serão substituídos pelo valor 1.

Clique OK.

Marque a caixa de seleção “Use no data when no range matches value” traduzindo, “quando nenhum intervalo corresponder ao valor, substitua por No Data”. Isso fará com que todo o resto do seu raster se torne células sem dados (NaN, No Data).

Execute o algoritmo.

O resultado será um raster onde os únicos pontos que existem são aqueles da sua rede de drenagem.

7. Capture os valores da elevação sobre rede de drenagem

Para fazer isso, basta multiplicar o MDE pelo raster gerado na etapa 5 usando o Raster Calculator.

Você deve obter um raster com as elevações sobre a rede de drenagem.

8. Interpolar nível de base

Interpolando a elevação da rede de drenagem para o resto do raster, obtemos um raster com o nível do fluxo fluvial de base (subterrâneo) aproximado, com base no limiar escolhido por você (!).

Para fazer isso, volto para o Prompt e uso novamente o SAGA:

saga_cmd grid_spline 5 -GRID channels.tif -TARGET_OUT_GRID interp.tif

para esta interpolação, eu escolhi Multilevel B-Spline from Grid Points do SAGA (nome no QGIS: Multilevel B-Spline Interpolation from raster), mas você pode usar qualquer algoritmo que desejar para esta interpolação. Este que eu escolhi aqui gera um raster visualmente contínuo.

9. Calcule a distância vertical até a rede de drenagem (Sim, o VDCN!)

Este é o último passo (provavelmente)!

Calcule a distância até a rede de drenagem executando your_DEM.tif - interp.tif no Raster Calculator do QGIS. Salve isso em um arquivo. Aperte em Run. Pronto!

10. (opcional) Remova os valores negativos

Algumas pessoas argumentam que VDCN não deveria ter valores negativos, mas acredito que seu raster final possa ter alguns valores negativos em alguns pontos. Se você for uma dessas pessoas, reclassifique o raster da etapa 9 em Reclassify by table (o mesmo algoritmo de antes), mas coloque isso em sua tabela:

Mínimo: -9999999

Máximo: 0

Valor: 0

Portanto, ele substituirá todos os números negativos por zero. Além disso, lembre-se de desmarcar a caixa de seleção “Use no data when no range matches value”.