Mapa de risco de erosão no QGIS com metodologia RUSLE

Da primeira vez que escrevi sobre esse tema aqui no blog, tratei o problema como um “índice de triagem”. Pesos arbitrários entre chuva e declividade, sem embasamento formal nenhum. Reli aquele texto depois de uns comentários no fórum apontando, com razão, que estava raso. Refiz do zero, com a metodologia que qualquer curso de conservação do solo ensina: a RUSLE, a versão revisada da USLE de Wischmeier e Smith.

A RUSLE não é perfeita. Foi desenvolvida para condições de clima temperado e depois adaptada para o Brasil por dezenas de estudos regionais, e isso importa para as escolhas metodológicas que vou justificar ao longo do texto. Mas é o padrão que qualquer engenheiro agrônomo, geógrafo ou avaliador ambiental vai reconhecer e conseguir questionar tecnicamente. É esse o ponto: quero que este mapa resista a uma banca, não só a um clique de curtir.

A equação é simples de escrever e complicada de implementar direito:

A = R × K × LS × C × P

A é a perda de solo estimada, em toneladas por hectare por ano. R é a erosividade da chuva. K é a erodibilidade do solo. LS combina comprimento de rampa e declividade. C é o uso e manejo do solo. P são as práticas conservacionistas. Vou construir cada um separadamente no QGIS e multiplicar os rasters no final.

Fator R: erosividade da chuva

O fator R representa a capacidade da chuva de causar erosão. O cálculo “correto”, segundo a literatura original, exige dados pluviográficos de alta resolução temporal (intensidade em intervalos de minutos) para calcular o índice EI30 chuva a chuva. Isso praticamente inviabiliza o trabalho para a maioria dos municípios brasileiros, porque essas séries pluviográficas são raras e concentradas em poucas estações históricas.

Na prática, o caminho viável, e amplamente usado na literatura nacional desde os anos 1990, é estimar o EI30 a partir de dados pluviométricos mensais e anuais, que são muito mais disponíveis. A equação empírica de Lombardi Neto e Moldenhauer (1992) resolve isso:

EI30(mês) = 67,355 × (p² / P)^0,85

R = Σ EI30(mês), para os 12 meses do ano

Aqui, p é a precipitação média do mês e P é a precipitação média anual, ambos em milímetros. Essa equação foi ajustada originalmente para Campinas, em São Paulo, e depois replicada e recalibrada em dezenas de outros estados. Se você está trabalhando fora do Sudeste, vale procurar se já existe uma equação de coeficiente de chuva calibrada para a sua região na literatura de ciência do solo local, antes de usar os coeficientes genéricos. É esse tipo de detalhe que separa um trabalho técnico defensável de um achismo com aparência de precisão.

Como implementar no QGIS:

1. Baixe as séries históricas de precipitação mensal da ANA (Agência Nacional de Águas e Saneamento) ou do INMET/BDMEP para todas as estações dentro e no entorno da sua área de estudo. Inclua estações fora do limite: a interpolação nas bordas fica mais confiável com pontos de apoio além do polígono.
2. Calcule p, P e o EI30 mensal para cada estação numa planilha, ou no Field Calculator do QGIS se você importar os dados já em formato tabular longo. Some os 12 valores mensais para obter o R anual de cada estação.
3. Importe as estações como camada de pontos com o valor de R como atributo. Reprojete para a UTM SIRGAS 2000 do fuso da área de estudo.
4. Interpole para toda a área usando Kriging. Aqui não use IDW: a escolha importa mais do que no post anterior, porque erosividade tem estrutura de autocorrelação espacial mais suave e contínua do que muitas variáveis pontuais, e o IDW tende a criar padrões em “olho de boi” ao redor de cada estação que não refletem a realidade física do fenômeno.

Fator K: erodibilidade do solo

K mede o quanto um tipo de solo específico é suscetível à erosão, independente da chuva ou da declividade. O cálculo de livro-texto usa o nomograma de Wischmeier et al. (1971), que combina porcentagem de silte mais areia muito fina, porcentagem de areia, teor de matéria orgânica, classe de estrutura e classe de permeabilidade do solo. Isso exige dados de laboratório que raramente estão disponíveis fora de um levantamento pedológico detalhado.

O caminho realista para o QGIS é outro: use o levantamento de solos disponível para a sua área (geralmente um shapefile de classes pedológicas da EMBRAPA, ou do zoneamento agroecológico estadual quando existir em escala mais detalhada) e atribua a cada classe de solo um valor de K obtido da literatura de calibração local. Valores de K para as principais classes de solo brasileiras já foram publicados em dezenas de estudos de bacias hidrográficas. Procure um trabalho de referência da sua região, ou de condições pedológicas semelhantes, antes de importar valores de tabelas americanas genéricas. Elas sistematicamente não refletem a mineralogia dos solos tropicais brasileiros, em especial os Latossolos, que se comportam de forma bem diferente dos solos temperados em termos de agregação e erodibilidade.

No QGIS, faça o join do valor de K à tabela de atributos do shapefile de solos (Propriedades, Junções, ou diretamente via Field Calculator criando um campo novo com expressão CASE WHEN por classe). Depois converta de vetor para raster (Conversão, Rasterizar) na mesma resolução espacial que você vai usar para os demais fatores. Esse alinhamento de resolução entre todos os rasters da equação é obrigatório antes da multiplicação final, e é o erro mais comum que vejo em projetos de RUSLE malfeitos.

Fator LS: comprimento de rampa e declividade

Esse é o fator mais sensível ao Modelo Digital de Elevação que você usar, e é onde a maioria dos tutoriais simplifica demais. A formulação clássica de Wischmeier e Smith assume rampas uniformes, o que raramente existe no terreno real. A versão que realmente vale a pena implementar é a de Desmet e Govers (1996), que calcula o LS a partir da área de contribuição por célula em vez do comprimento de rampa simples. Isso resolve boa parte da distorção em terrenos com convergência e divergência de fluxo.

No QGIS, dentro da Caixa de Ferramentas de Processamento:

1. Obtenha um Modelo Digital de Elevação de resolução adequada. Para trabalho municipal ou de bacia, o Topodata do INPE (refinamento do SRTM para 30 metros) é o mínimo aceitável. Se você tiver acesso a um MDE derivado de LiDAR ou de VANT para a área específica, use-o: a diferença na qualidade do fator LS é substancial.
2. Preencha depressões espúrias do MDE antes de qualquer cálculo hidrológico, usando “SAGA – Fill Sinks (Wang & Liu)” ou o algoritmo equivalente do GRASS (r.fill.dir). Pular esta etapa produz artefatos graves no cálculo de fluxo acumulado.
3. Rode “SAGA – LS-Factor” (o algoritmo já implementa a formulação de Desmet & Govers) apontando o MDE preenchido. Ou, se preferir, calcule fluxo acumulado com “SAGA – Flow Accumulation (Top-Down)” e declividade com “GRASS – r.slope.aspect” e combine manualmente pela fórmula de Moore e Burch (1986), no Raster Calculator:

LS = (fluxo_acumulado × tamanho_da_célula / 22,13)^0,4 × (sin(declividade_radianos) / 0,0896)^1,3

Se optar pela fórmula manual, cuidado com duas armadilhas. A declividade precisa estar em radianos, não em graus. O QGIS e o GRASS por padrão entregam em graus, então é preciso converter multiplicando por pi sobre 180 antes de aplicar o seno. E o fluxo acumulado do SAGA ou do GRASS geralmente vem em número de células, não em área. Multiplique pelo tamanho de célula ao quadrado se o algoritmo que você usou não fizer essa conversão internamente. Confira a documentação do algoritmo específico da sua versão antes de assumir a unidade.

Fator C: uso e manejo do solo

C representa o quanto a cobertura vegetal e o manejo agrícola protegem o solo contra a erosão, numa escala de 0 (proteção total) a 1 (solo completamente exposto, sem nenhuma proteção).

No QGIS, baixe o raster de uso e cobertura do MapBiomas para o ano de referência do seu estudo, recorte para a área de interesse (Raster, Extração, Cortar Raster pela Camada de Máscara) e reclassifique cada classe do MapBiomas para um valor de C usando a ferramenta Reclassificar por Tabela (Raster, Análise Raster, Reclassificar por Tabela). Os valores de referência mais usados na literatura brasileira giram em torno disto: floresta nativa e formações savânicas próximo de 0,001 a 0,005; pastagem em bom estado de manejo entre 0,01 e 0,1, com variação grande dependendo da cobertura; culturas anuais entre 0,1 e 0,5 conforme o ciclo e a cobertura do solo na época chuvosa; solo exposto ou área urbanizada próximo de 1 na porção não impermeabilizada, com classe separada (geralmente excluída do cálculo ou tratada à parte) para área efetivamente impermeabilizada.

Um ponto que costuma passar batido: C não é estático ao longo do ano em áreas agrícolas, porque a cobertura do solo muda com o calendário de plantio. Se o seu objetivo é um diagnóstico anual médio, o MapBiomas é suficiente, já que é uma classificação anual consolidada. Se o objetivo é identificar picos de vulnerabilidade, por exemplo solo exposto entre a colheita e o plantio da safra seguinte coincidindo com o início do período chuvoso, você precisa de uma abordagem mensal ou sazonal com séries de imagens de satélite de maior frequência temporal. Isso foge do escopo deste tutorial, mas vale deixar registrado como limitação.

Fator P: práticas conservacionistas

P representa o efeito de práticas de manejo que reduzem a erosão sem alterar a cobertura vegetal em si: terraceamento, plantio em contorno, cultivo em faixas. Na ausência de um cadastro confiável dessas práticas para a sua área, o que é a realidade da maioria dos municípios brasileiros fora de grandes projetos agrícolas mapeados, o procedimento padrão e aceito na literatura é assumir P igual a 1 para todo o raster. Isso resulta numa estimativa conservadora, no sentido de “pior cenário”, superestimando levemente a perda de solo em áreas que de fato têm alguma prática de manejo não mapeada. Documente essa premissa explicitamente no seu relatório ou no rodapé do mapa publicado. É o tipo de simplificação que precisa estar visível, não escondida.

Combinando os fatores no Raster Calculator

Com R, K, LS, C e P como rasters alinhados na mesma resolução, mesmo CRS e mesma extensão espacial, a etapa final é simplesmente:

A = “R@1” * “K@1” * “LS@1” * “C@1” * “P@1”

no Raster Calculator do QGIS. O resultado sai em toneladas por hectare por ano. Antes de rodar, confirme visualmente, no Alinhador de Rasters ou pelas propriedades de cada camada, que todos têm exatamente a mesma resolução de célula e a mesma origem de grade. Rasters “quase alinhados”, com uma pequena diferença de origem, produzem multiplicação incorreta célula a célula sem gerar nenhum erro visível no QGIS. Esse é provavelmente o defeito mais comum e mais silencioso em mapas de RUSLE malfeitos que já revisei.

Classificação final e leitura do resultado

Classifique o raster de perda de solo em faixas de severidade. A literatura brasileira de conservação do solo costuma usar algo próximo disto: menor que 10 toneladas por hectare por ano como tolerável, entre 10 e 20 como moderada, entre 20 e 50 como alta, acima de 50 como muito alta a crítica. Ajuste essas faixas para a tolerância de perda de solo específica do tipo de solo predominante na sua área, se esse dado estiver disponível na literatura pedológica regional, em vez de aplicar um corte genérico sem verificar se faz sentido para o contexto local.

Validação, a parte que a maioria dos tutoriais pula

Um mapa de RUSLE sem nenhuma forma de validação é uma estimativa de crença, não um produto técnico. Duas validações mínimas dá para fazer sem sair a campo.

Validação cruzada da interpolação do fator R. Ao rodar o Kriging no QGIS via SAGA, ative a opção de validação cruzada leave-one-out quando disponível na versão do algoritmo, ou rode manualmente: remova uma estação por vez, interpole sem ela, compare o valor previsto no ponto removido com o valor observado, e calcule o erro médio quadrático do conjunto. Isso dá uma medida objetiva de quão confiável é a superfície de erosividade que você gerou. É exatamente o tipo de número que qualquer TCC ou parecer técnico sério sobre o tema costuma pedir.

Comparação com pontos de erosão conhecidos. Se você tem acesso a registros de ocorrências de erosão já documentadas na área (voçorocas mapeadas, laudos de Defesa Civil, cadastros municipais de áreas de risco), sobreponha esses pontos ao mapa classificado e verifique se eles caem predominantemente nas classes de maior severidade. Não é uma validação estatística formal, mas é um teste de sanidade rápido que frequentemente revela problemas grosseiros de metodologia antes de qualquer coisa ser publicada.

O que eu deixaria escrito, em destaque, junto ao mapa publicado

• Ano-base de cada camada de entrada (série pluviométrica, cobertura MapBiomas, levantamento de solos), porque isso muda a interpretação do resultado.
• A premissa de P igual a 1, se for o caso.
• A resolução espacial final do raster e a fonte do MDE usado.
• Uma frase explícita: isso é uma estimativa de potencial de perda de solo por processo laminar e em sulcos. Não inclui erosão por voçorocamento em estágio avançado nem movimentos de massa. Esses processos exigem modelagem geotécnica separada.