En un tutorial anterior aprendimos a calcular el Índice de Erosividad de la lluvia, uno de los factores que intervienen en la Ecuación Universal de la Pérdida de suelo (USLE). En esta ocasión, vamos a realizar un análisis de erosión con ArcGIS, para lo cual calcularemos todos los parámetros que intervienen en dicha ecuación.
Análisis de la Erosión con ArcGIS en una cuenca hidrológica
La erosión puede desencadenar la disminución de las zonas cultivables, la fertilidad de los suelos, la capacidad de retención de zonas inundables o incluso puede la colmatación y eutrofización de masas de agua. Por tanto, el estudio de la erosión del suelo resulta especialmente importante en ámbitos como la agricultura, la explotación forestal, la hidrología y demás ciencias relacionadas con el medio ambiente.
Uno de los métodos más utilizados para su cálculo y que ha tenido mayor aceptación y difusión ha sido la Ecuación Universal de la Pérdida de Suelo (USLE), que establece seis factores diferentes como los responsables de las pérdidas anuales de suelo:
A=R∙K∙L∙S∙C∙P
Las pérdidas anuales de suelo A [t/ha.año] se calculan en base a un índice de erosividad de la lluvia R [MJ.mm/ha.h], la erosionabilidad del suelo K [t/ha por unidad de R] , un factor topográfico producto de los factores longitud de la ladera L [adimensional] y su pendiente S [adimensional], un factor de cubierta vegetal C [adimensional] y un factor de prácticas de conservación de suelos o prácticas de manejo P [adimensional].
Esquema del proceso
Nosotros vamos a obtener y analizar el riesgo de erosión de una cuenca hidrológica y el proceso que vamos a seguir para obtener cada uno de los factores que intervienen en la USLE es el siguiente:
- Factor R: contamos con un raster de la erosividad de la lluvia en nuestro ámbito de estudio. Este raster se ha obtenido mediante el cálculo del Índice Modificado de Fournier (IMF) (calculado a partir de las precipitaciones mensuales) al que se le ha aplicado la ecuación genérica del I.C.O.N.A. para estimar el factor de erosividad de la lluvia.
- Factor K: disponemos de una capa vectorial que contiene los valores de erosionabilidad del suelo en nuestro ámbito de estudio de manera que habrá que obtener el raster correspondiente.
- Factor LS: a partir del MDT calcularemos de forma independiente la pendiente (factor S) y la longitud de ladera (factor L) y una vez obtenidos aplicaremos la fórmula de Moore & Burch (1986) para generar el raster de factores topográficos.
- Factor C: disponemos de una capa vectorial que contiene el factor de cubierta vegetal de nuestra cuenca. Para obtener esta capa se ha partido de la información de usos de suelo en la región de estudio y se han utilizado los valores del factor C en función del tipo de cubierta establecidos por el I.C.O.N.A. en 1983 . Por tanto, lo que tendremos que hacer es obtener el raster correspondiente.
- Factor P: en nuestra zona no existen prácticas de manejo, de manera que asumiremos que el factor P es 1.
El objetivo final es obtener el raster de la pérdida de suelo que tiene lugar en nuestra cuenca hidrológica por lo que debemos disponer de la información de cada uno de los factores en formato raster.
Índice de Erosionabilidad (K)
El factor K (erosionabilidad/erodabilidad del suelo) refleja su vulnerabilidad a la erosión y se calcula en base a sus propiedades físicas. Nosotros contamos con una capa vectorial que ya contiene la estimación de estos valores así que lo único que debemos hacer es pasar esta información a formato raster.
Para ello, utilizamos la herramienta de conversión para trasformar esta capa de polígonos (formato vectorial ) a formato raster ArcToolbox < Conversion Tools < To Raster < Polygon to Raster.
Seleccionamos la capa que queremos convertir y el campo que queremos representar (‘Value Field’), que en este caso, será el campo denominado K. Por último indicamos un nombre para el raster de salida.
Tras aceptar, obtenemos un raster que recoge los valores de K para el ámbito de estudio; dichos valores se mueven en el rango de 0,15 a 0,4.
Factores Topográficos (LS)
Los factores topográficos se refieren a la longitud de la ladera (L) y su pendiente (S) y ambos forman el factor conocido como LS. Nosotros vamos a calcular el factor LS mediante la fórmula de Moore & Burch (1986):
LS = (Flow accumulation*(cell size)⁄22,13)^(0,4 ) ∙ 〖(sinslope⁄0,0896)〗^1,3
donde Flow Accumulation es el número de celdas que contribuyen al flujo en una celda dada, cell size es la longitud del tamaño de un lado de las celdas y sin slope es el seno de la pendiente en radianes.
Pendiente (factor S)
En primer lugar calcularemos la pendiente en grados para después transformarla a radianes. Utilizamos ahora la herramienta específica para el cálculo de pendientes y para ello es necesario que dispongamos del MDT sobre el que se van a calcular; en nuestro caso contamos con el MDT de la zona de estudio; lo cargamos en el proyecto y seleccionamos la herramienta «Slope» ArcToolbox < Spatial Analyst Tools < Surface < Slope.
Seleccionamos el MDT y asignamos un nombre a la capa de salida. La herramienta nos permitirá calcular la pendiente en grados (‘Degree’) o en % (‘Percent Rise’); seleccionamos grados.
Este raster recoge la pendiente (en grados) pero la fórmula que vamos a emplear en el cálculo del factor LS requiere que la pendiente esté calculada en radianes así que tenemos que multiplicar la pendiente obtenida por π/180. Abrimos la calculadora raster para escribir la expresión correspondiente ArcToolbox < Spatial Analyst Tools < Map Algebra < Raster calculator.
Longitud de ladera (factor L)
El cálculo de la longitud de ladera se realiza en base al raster de acumulación de flujo. Este mapa de acumulación de flujo representa las celdas en las que se acumula el agua al fluir desde las celdas con mayor valor de altitud.
Para obtenerlo habrá que seguir los siguientes pasos:
a) Relleno de sumideros:
Un sumidero es una celda o un conjunto de celdas conectadas espacialmente a la que no se le puede asignar una dirección de flujo, por lo tanto para garantizar la representación correcta de cuencas y arroyos hay que realizar el relleno de sumideros. Abrimos la herramienta correspondiente ArcToolbox < Spatial Analyst Tools < Hydrology < Fill y lo único que debemos hacer es seleccionar el raster a rellenar (‘Input surface raster’), es decir, nuestro MDT, y asignar un nombre al nuevo raster que se va a crear (‘Output surface raster’). Al hacerlo obtenemos un nuevo raster aparentemente idéntico al MDT que habíamos generado pero aunque no lo apreciaremos, este nuevo raster será el MDT corregido, donde se habrán rellenado los sumideros en caso de que existieran.
b) Mapa de Direcciones de Flujo
El siguiente paso consiste en obtener el Mapa de Direcciones de Flujo. A partir de un MDT que contiene los datos de altitud de cada celda, se obtiene la dirección que seguirá el flujo y se le asigna un valor concreto en función de la dirección que presenta.
Obtenemos el mapa de direcciones de flujo de nuestra cuenca ArcToolbox < Spatial Analyst Tools < Hydrology < Flow Direction. Indicamos el MDT sobre el que se quiere obtener el mapa de direcciones de flujo (‘Input surface raster’); en este caso debemos especificar que es el denominado «Fill» ya que este raster es el que tiene los sumideros corregidos y como siempre asignamos un nombre a la capa de salida (‘Output flow direction raster’).
c) Mapa de acumulación de Flujo
El mapa de acumulación de flujo permite conocer cuáles son las celdas en las que se acumula el agua al fluir desde las celdas con mayor valor de altitud. Así, las celdas con acumulación de flujo altas serán áreas donde el flujo se concentra, y las celdas con valores de acumulación de flujo igual a 0 serán alturas topográficas locales. Procedemos a realizar el análisis ArcToolbox < Spatial Analyst Tools < Hydrology < Flow Accumulation
Indicamos el mapa de direcciones que vamos a utilizar para su cálculo (‘Input flow direction raster’) y asignamos un nombre a la capa de salida (‘Output accumulation raster’). Dejamos por defecto «Float» como tipo de dato de salida (‘Output data type’). Obtenemos el mapa donde observamos en negro las zonas con una acumulación de flujo 0 (es decir, zonas donde no existe acumulación de flujo) y en blanco aquellas zonas donde sí se produce la acumulación.
Factor LS
Volvemos a la fórmula de cálculo del factor LS:
LS = (Flow accumulation ∙(cell size) ⁄ 22,13)^0,4 ∙ 〖(sinslope ⁄ 0,0896)〗^1,3
La fórmula implica que hay que multiplicar la acumulación de flujo por el tamaño de la celda lo que representa la longitud de escurrimiento. Nuestro máximo valor de acumulación de flujo es de 235426 celdas, lo que multiplicado por el tamaño de la celda (50 m) daría como resultado una longitud máxima de escurrimiento de 11.771.300 metros. Utilizar estos datos sería erróneo pues estaríamos sobreestimando el valor de la longitud de pendiente y podríamos obtener valores exagerados de erosión. Para evitar esta sobreestimación debemos restringir la longitud estableciendo por ejemplo una longitud máxima de 250 metros, lo que equivale a 5 celdas. Así que debemos obtener un raster de acumulación de flujo cuyo valor máximo sea 5. Esto significa que las celdas que en el raster de acumulación tienen un valor menor de 5 deben conservar su valor original pero aquellas que presentan un valor mayor, deberán tomar todas el valor de 5 (valor máximo que hemos decidido fijar).
Lo que tendremos que hacer será crear varias capas a partir de la original:
- Capa A: Si Fac ≤5, tomará el valor de 1 y si Fac > 5, tomará el valor 0.
- Capa B: Si Fac ≤ 5, tomará el valor de 0 y si Fac > 5 tomará el valor 5.
Después lo que haremos será multiplicar la Fac original por la Capa A con lo que conseguiremos obtener un raster que conservará los valores originales a excepción de aquellas celdas donde el valor original era mayor de 5, que ahora será 0. Una vez hecho esto, a esta nueva capa le sumaremos la Capa B de manera que esas celdas que tienen un valor de 0 tomarán el valor de 5 que es el máximo que queremos fijar.
Para obtener las que hemos denominado «Capa A» y «Capa B» tendremos que «modificar» los valores del raster de Fac original para lo que vamos a realizar una reclasificación ArcToolbox < Spatial Analyst Tools < Reclass < Reclassify.
Seleccionamos el raster de acumulación de flujo (‘Fac’) y como campo a reclasificar (‘Reclass Field’) tenemos que elegir el valor (‘Value’). El programa agrupará por defecto los valores originales en varios rangos (‘Old Values’) y nosotros debemos asignar los nuevos (‘New values’). Lo que haremos será agrupar los valores en dos categorías:
- Fac < 5, es decir de 0 a 5
- Fac > 5, es decir de 5 al valor máximo que toma el raster original.
Para hacer esta nueva agrupación pulsamos sobre ‘Classify…’.
Seleccionamos dos clases e indicaremos que la primera de ellas será el valor 5 dejando el máximo que aparece. Pulsamos en ‘OK’ y el programa nos devuelve a la pantalla anterior. Ahora los ‘Old values’ se agrupan en dos rangos diferentes: 0 – 5 y 5 – 235426.
Tal y como hemos explicado lo que queremos es que la Fac menor o igual a 5 (‘Old values’ 0 – 5) tomen el valor de 1 y la Fac mayor de 5 (‘Old values’ 5 – 235426) tomen el valor de 0 . Para conseguirlo cambiamos los New values y asignamos el nombre a la capa de salida.
Se genera un raster en el que cada celda únicamente presenta valores 0 o 1. El valor 1 son aquellas celdas cuyo Fac era ≤ 5 y el valor 0 son aquellas que presentaban Fac > 5.
Para generar la “Capa B” hacemos de exactamente lo mismo que anteriormente, establecer dos clasificaciones para el Fac, 5 y el valor máximo. Lo que será diferente es que en esta ocasión daremos un nuevo valor de «0» a aquellos que estaban por debajo de 5 y un valor de 5 a todos aquellos que estaban por encima.
En este caso el raster de 0 y 1 representa con un valor de 0 todas las celdas que en la capa original tenían un Fac menor o igual a 5 y un valor de 5 para todas aquellas que en el original tenían Fac por encima de 5.
Ahora que tenemos la capa original, la capa A y la capa B, volvemos a la calculadora raster y realizamos la operación que explicamos anteriormente:
Fac original × Capa A + Capa B
Obtenemos la capa final que es la que utilizaremos para calcular el factor LS.
Ahora únicamente tendremos que utilizar la calculadora raster de nuevo para calcular el factor LS mediante la fórmula correspondiente:
LS = (Flow accumulation ∙(cell size) ⁄ 22,13)^0,4 ∙ 〖(sinslope ⁄ 0,0896)〗^1,3.
Debemos tener cuidado para utilizar las capas correctas: como «Flow accumulation» tendremos que usar la última capa que hemos obtenido (la que contiene la acumulación de flujo limitada a un máximo de 5 celdas) y como «slope» tendremos que utilizar la capa que contiene la pendiente en radianes. Por último asignaremos un nombre al raster de salida.
Así hemos obtenido el raster que refleja el factor LS en nuestra cuenca y que presenta valores de 0 a 48.
Índice de Cobertura Vegetal (C)
El factor C representa la efectividad de las plantas como cubierta protectora del suelo frente a la energía de impacto de las gotas de lluvia y a la fuerza del flujo superficial.
Para nuestro estudio nos han facilitado una capa vectorial que contiene varios polígonos y para cada uno de ellos tenemos información del uso de suelo que presenta y el dato del factor C. Por lo tanto lo único que debemos hacer es obtener el raster correspondiente ArcToolbox < Conversion Tools < To Raster < Polygon to Raster.
Seleccionamos la capa e indicamos el campo por el que queremos convertirla, que en este caso es el denominado “Factor_C”.
El raster obtenido refleja el factor C en nuestra cuenca cuyos valores se encuentran en el rango de 0 a 0,4.
Prácticas de manejo (P)
Según nos han indicado, no se han detectado prácticas de manejo en la región por lo que asumiremos que el factor P es 1 para todo el ámbito, así que no necesitamos crear ningún raster que represente este factor.
APLICACIÓN DE LA ECUACIÓN USLE
Como ya sabemos, las pérdidas anuales de suelo A [t/ha.año] se calculan en base al índice de erosividad de la lluvia R [MJ.mm/ha.h], la erosionabilidad del suelo K [t/ha por unidad de R], el factor topográfico producto de los factores longitud de la ladera L [adimensional] y su pendiente S [adimensional], el factor de cubierta vegetal C [adimensional] y el factor de prácticas de conservación de suelos P [adimensional].
A=R∙K∙L∙S∙C∙P
Ya disponemos del raster de cada uno de los factores así que vamos a aplicar la ecuación para obtener la pérdida de suelo que se produce en la cuenca.
Utilizamos una vez más la calculadora raster ArcToolbox < Spatial Analyst Tools < Map Algebra < Raster Calculator, para escribir la fórmula teniendo en cuenta que el factor P es igual a la unidad por lo que no será necesario incluirlo en el cálculo.
Ahora, para cuantificar el rango de pérdida de suelo vamos a utilizar alguna de las clasificaciones existentes a este respecto. Por ejemplo, vamos a editar la leyenda para establecer la clasificación por niveles erosivos establecida en el «Mapa del grado de erosión Hídrica de los suelos de la CAPV» (Gobierno Vasco, 2005), que establece lo siguiente:
A la vista de los resultados podemos concluir que se trata de una región en la que los procesos erosivos tienen una gran importancia. En general, toda la zona presenta una tasa anual de pérdida de suelo considerable, especialmente marcada en la zona norte, coincidiendo con las zonas donde la erosividad de la lluvia es mayor, las pendientes son más escarpadas y existe una menor cobertura vegetal.
< ¿QUIERES APLICAR TUS CONOCIMIENTOS DE GIS EN EL SECTOR AMBIENTAL?
Échale un vistazo a nuestro Curso ArcGIS aplicado a la Gestión Ambiental