Preparación de Datos Geoespaciales para el Modelado Ambiental con InVEST

Procesamiento de Datos para Módulos Hídricos y de Erosión en InVEST

Este documento detalla los procedimientos experimentales para la preparación y el procesamiento de datos requeridos por los módulos del modelo InVEST, centrándose en el rendimiento hídrico, la retención de agua y la erosión del suelo. Se comparten los pasos de experimentación y se agradece cualquier sugerencia para su mejora.

  1. Módulo de Rendimiento Hídrico

1.1. Configuración del Entorno

Se recomienda establecer un espacio de trabajo claro y organizar los archivos por sufijos para facilitar el procesamiento por lotes.

1.2. Datos de Precipitación

Para la precipitación, se pueden utilizar conjuntos de datos mensuales con una resolución de 1 km para China (1901-2023), disponibles en el Centro Nacional de Datos Científicos del Sistema Terrestre (www.geodata.cn). Los datos originales provienen de CRU, y la unidad de precipitación es 0.1 mm. Complementariamente, se pueden emplear datos mensuales de temperatura media cercana a la superficie de 1 km para China (2000-2020), derivados de ERA5 y disponibles en el mismo centro (www.geodata.cn).

1.3. Datos de Evapotranspiración Potencial (ETP)

Se pueden obtener conjuntos de datos mensuales de ETP para China (1901-2023) con una resolución espacial de aproximadamente 1 km (www.geodata.cn). Al igual que la precipitación, la unidad es 0.1 mm.

1.4. Procesamiento de Datos Climáticos (Precipitación y ETP)

Para la ingestión de datos en InVEST, es común que los archivos NetCDF deban convertirse a un formato de imagen ráster, como GeoTIFF. El siguiente script de MATLAB facilita esta conversión, extrayendo datos mensuales y aplicando factores de escala.

Asegúrese de ajustar la ruta del directorio de trabajo, la ruta de un archivo GeoTIFF de ejemplo (para obtener la referencia espacial) y el nombre de la variable dentro del archivo NetCDF (por ejemplo, 'pre' para precipitación, 'tmp' para temperatura, 'etp' para evapotranspiración).


clc;
clear;

% --- Configuración de Rutas y Archivos ---
rutaTrabajo = 'E:\PET\AOrigion'; % Directorio que contiene los archivos .nc
muestraReferencia = 'E:\PET\Sample2.tif'; % Ruta a un GeoTIFF de muestra para la referencia espacial

% --- Lectura de Información de Referencia Espacial ---
% Utilizar el archivo de muestra para obtener la información de georreferenciación
[~, R] = geotiffread(muestraReferencia);
infoGeotiff = geotiffinfo(muestraReferencia);
dimensiones = infoGeotiff.RasterSize; % [filas, columnas]

% --- Preparación de la Carpeta de Salida ---
cd(rutaTrabajo);
archivosNc = dir('pet_*.nc'); % Buscar archivos NetCDF que comiencen con 'pet_'

carpetaSalida = fullfile(rutaTrabajo, 'output_etp_mensual');
if ~exist(carpetaSalida, 'dir')
    mkdir(carpetaSalida);
end

% --- Procesamiento de Archivos NetCDF ---
for i = 1:length(archivosNc)
    [~, nombreBase, ~] = fileparts(archivosNc(i).name);
    anio = str2double(nombreBase(5:8)); % Extraer el año del nombre del archivo (ej. "pet_2014.nc")
  
    % Construir la ruta completa al archivo .nc
    rutaArchivoNc = fullfile(rutaTrabajo, archivosNc(i).name);
    
    % Leer variables del archivo NetCDF
    % Ajuste el nombre de la variable ('etp', 'pre', 'tmp') según sea necesario
    evapotranspiracion = ncread(rutaArchivoNc, 'etp'); 
    
    % Manejar valores nodata (si aplica)
    evapotranspiracion(evapotranspiracion == -32768) = NaN;
    
    % Aplicar factor de escala (ej. 0.1 para convertir de 0.1 mm a mm)
    % Si no hay factor de escala, comente o elimine esta línea
    evapotranspiracion = evapotranspiracion * 0.1; 

    % Extraer y guardar datos mensuales
    numMeses = size(evapotranspiracion, 3); % Número de capas (meses)
    
    for idxMes = 1:numMeses
        % Extraer la capa de datos para el mes actual
        datosMensuales = evapotranspiracion(:, :, idxMes)'; % Transponer si es necesario
                                                            % para que coincida con el formato geotiffread
        
        % Definir el nombre del archivo de salida
        nombreArchivoSalida = sprintf('etp_%d_mes_%02d.tif', anio, idxMes);
        rutaArchivoSalida = fullfile(carpetaSalida, nombreArchivoSalida);

        % Escribir el GeoTIFF
        geotiffwrite(rutaArchivoSalida, datosMensuales, R, ...
            'TiffType', 'bigtiff', 'GeoKeyDirectoryTag', infoGeotiff.GeoTIFFTags.GeoKeyDirectoryTag);      
    end
end
disp('Procesamiento de datos climáticos completado.');

1.5. Profundidad de la Capa Limitante Radicular

Este dato, necesario en milímetros (mm), se puede obtener de diversas fuentes:

  • SoilGrids250m 2017-03 - Depth to bedrock (R horizon): Disponible en ISRIC.
  • Mapa de profundidad del lecho rocoso de China con una resolución espacial de 100 metros: Datos disponibles en figshare, con su publicación asociada en Scientific Data. Es crucial recordar que estos datos están en metros (m) y deben multiplicarse por 1000 para convertirlos a milímetros.

1.6. Contenido de Agua Disponible para las Plantas (CADP)

1.6.1. Cálculo Basado en Documentación Oficial

Una aproximación es utilizar los datos de SoilGrids250m 2017-03 - Derived available soil water capacity (volumetric fraction) until wilting point.

  1. Descarga: La descarga de estos archivos puede ser lenta.

  2. Recorte: Renombre los archivos según la nomenclatura requerida y recorte las áreas de interés, preferiblemente con un búfer. Los archivos resultantes deben seguir un patrón como: AWC_sl1_clip.tif, AWC_sl2_clip.tif, ..., AWC_sl7_clip.tif.

  3. Cálculo: Cargue todos los archivos recortados en su software SIG (ej. ArcGIS) y utilice la Calculadora Ráster con la siguiente fórmula: ```arcgis

    0.005 * 0.5 * (5 * ("AWC_sl1_clip.tif" + "AWC_sl2_clip.tif") + ... 10 * ("AWC_sl2_clip.tif" + "AWC_sl3_clip.tif") + ... 15 * ("AWC_sl3_clip.tif" + "AWC_sl4_clip.tif") + ... 30 * ("AWC_sl4_clip.tif" + "AWC_sl5_clip.tif") + ... 40 * ("AWC_sl5_clip.tif" + "AWC_sl6_clip.tif") + ... 100 * ("AWC_sl6_clip.tif" + "AWC_sl7_clip.tif"))

  4. Adaptación a InVEST: El resultado debe estar en el rango de 0 a 100. Conviértalo a porcentaje (si es necesario) y realice los siguientes pasos para su compatibilidad con InVEST:

    • Reproyección: A un sistema de coordenadas adecuado (ej. UTM WGS84 48N).
    • Remuestreo: A la resolución deseada (ej. 30 m) utilizando un método como el bilineal.
    • Extracción por Máscara: Asegurando que la extensión y el número de filas/columnas coincidan con los datos de uso del suelo.

1.6.2. Estimación Mediante Fórmula Empírica

Las fórmulas empíricas ofrecen una alternativa para estimar el CADP, como las propuestas por Zhou Wenzo et al. (2003, 2005) para China. Estas fórmulas requieren datos de porcentaje de arena (SAND), limo (SILT), arcilla (CLAY) y carbono orgánico (OC) del suelo. Si se dispone de múltiples capas de suelo (ej. HWSD 2.0 con 7 capas), se recomienda utilizar un promedio ponderado.


clc;
clear;

% --- Configuración de Rutas y Archivos de Entrada/Salida ---
rutaBase = 'F:\Gannan\Gannan_auxiliary\HWSD1_Gannan\ccopy\';
nombreArchivoSalida = [rutaBase, 'CADP_T.TIF'];

% Rutas a los archivos GeoTIFF de propiedades del suelo
rutaArena = [rutaBase, 'HWSD1_T_SAND.tif'];
rutaLimo = [rutaBase, 'HWSD1_T_SILT.tif'];
rutaArcilla = [rutaBase, 'HWSD1_T_CLAY.tif'];
rutaMateriaOrganica = [rutaBase, 'HWSD1_T_OM_oc1.724.tif']; % O rutaOC si tiene carbono orgánico directo

% --- Lectura y Preparación de Datos ---
valorNoData = -128; % Valor que representa "no data" en los rásteres

[arena, R] = geotiffread(rutaArena);
arena(arena == valorNoData) = 0; % Asignar 0 a valores nodata

[limo, R] = geotiffread(rutaLimo);
limo(limo == valorNoData) = 0;

[arcilla, R] = geotiffread(rutaArcilla); 
arcilla(arcilla == valorNoData) = 0;

[materiaOrganica, R] = geotiffread(rutaMateriaOrganica);
materiaOrganica(materiaOrganica < 0) = 0; % Asegurar que la MO no sea negativa

% --- Cálculo del CADP (PAWC) con Fórmula Empírica ---
% Fórmula basada en estudios como Zhou Wenzuo et al.
% CADP = 54.509 - 0.132*Arena - 0.003*Arena^2 - 0.055*Limo - 0.006*Limo^2 - 0.738*Arcilla + 0.007*Arcilla^2 - 2.688*MO + 0.501*MO^2
cadpCalculado = 54.509 ...
             - 0.132 .* arena - 0.003 .* arena.^2 ...
             - 0.055 .* limo - 0.006 .* limo.^2 ...
             - 0.738 .* arcilla + 0.007 .* arcilla.^2 ...
             - 2.688 .* materiaOrganica + 0.501 .* materiaOrganica.^2;

% --- Manejo de Valores Anómalos ---
% Los resultados del CADP deben estar dentro de un rango esperado (ej. 0-100).
% Si hay valores excepcionalmente altos debido a altos niveles de materia orgánica,
% se pueden limitar utilizando funciones como 'Con' en la Calculadora Ráster (ej. Con(PAWC_CAL > 30, 30, PAWC_CAL)).
% Opcionalmente, se puede aplicar una escala si los valores son muy grandes.

% --- Escritura del Archivo GeoTIFF de Salida ---
geotiffwrite(nombreArchivoSalida, cadpCalculado, R);
disp('Cálculo del Contenido de Agua Disponible para las Plantas completado.');

1.6.3. Manejo de Valores Vacíos

Si el CADP calculado presenta áreas sin datos, se puede crear una capa de ceros (ej. convirtiendo un shapefile del área de estudio a ráster y reclasificándolo a cero). Luego, utilice la herramienta "Mosaico a Nuevo Ráster" (Mosaic To New Raster) en su SIG, configurando la propiedad "First" para la superposición, asegurándose de que la capa de CADP esté por encima de la capa de ceros para rellenar los vacíos.

1.7. Datos de Uso del Suelo

La preparación de los datos de uso del suelo es un paso estándar y no se detallará aquí.

1.8. Tabla de Parámetros Biofísicos

La tabla biofísica es fundamental para el módulo de rendimiento hídrico. Consulte el manual de InVEST para conocer los campos requeridos y proporcione los datos correspondientes para cada tipo de uso del suelo.

1.9. Parámetro Z

El parámetro Z es un factor empírico que requiere calibración. Su valor debe ajustarse iterativamente basándose en la literatura local o meidante un proceso de validación.

1.10. Datos de Cuencas Hidrográficas

Utilice los límites de su área de estudio (ej. límites de municipios o condados fusionados). Es esencial añadir un campo de tipo entero llamado "ws_id" a la tabla de atributos de este shapefile.

1.11. Datos de Subcuencas

Para la delimitación de subcuencas, se pueden emplear modelos hidrológicos especializados como el modelo SWAT.

1.12. Ejecución (RUN)

Antes de ejecutar el modelo, realice una verificación final de todos los datos de entrada:

  • Precipitación anual total
  • Evapotranspiración potencial anual total
  • Profundidad de la capa limitante radicular
  • Contenido de agua disponible para las plantas
  • Datos de uso del suelo
  • Tabla de parámetros biofísicos
  • Parámetro Z
  • Datos de cuencas hidrográficas
  • Datos de subcuencas (si aplica)

1.13. Validación

La validación se realiza comparando los resultados del modelo con datos observacionales, como los informes oficiales de recursos hídricos. Se puede construir una tabla en Excel para probar diferentes valores del parámetro Z y evaluar la precisión del modelo.

  1. Módulo de Retención de Agua

2.1. Conductividad Hidráulica Saturada del Suelo (Ks)

La Ks es un insumo clave para este módulo. Se puede calcular utilizando software como SPAW (Soil-Plant-Atmosphere-Water), disponible en NRCS USDA.

El proceso incluye:

  1. Recortar los datos de propiedades del suelo (ej. BIL) a su área de estudio para generar archivos TIFF.
  2. Crear un mapa de valores únicos para las combinaciones de propiedades del suelo.
  3. Utilizar una tabla de búsqueda (ej. con VLOOKUP o MATLAB) para vincular los porcentajes de arena y arcilla (SAND y CLAY) con los valores de Ks calculados previamente en SPAW para cada combinación única.

2.2. Índice de Humedad Topográfica (TWI)

El TWI (Topographic Wetness Index) se calcula a partir de un Modelo Digital de Elevación (DEM):

  1. Rellenar Sumideros (Fill Sinks): Aplique esta herramienta al DEM para eliminar pequeños errores topográficos.

  2. Dirección de Flujo (Flow Direction): Calcule la dirección del flujo para cada celda a partir del DEM sin sumideros.

  3. Acumulación de Flujo (Flow Accumulation): Calcule el número de celdas que fluyen hacia cada celda, utilizando el resultado de la dirección de flujo.

  4. Calculadora Ráster: Utilice la siguiente fórmula para calcular el TWI: ```arcgis

    TWI = ln ( [Acumulación de Flujo] / tan ( [Pendiente en radianes] ) )

    
 Donde SCA (Specific Catchment Area) es la acumulación de flujo, y la pendiente debe estar en radianes.

  5. Módulo de Erosión del Suelo (RUSLE)

3.1. Configuración del Entorno

Al igual que en otros módulos, establezca un espacio de trabajo organizado y considere el uso de sufijos para el procesamiento por lotes de archivos.

3.2. Modelo Digital de Elevación (DEM)

El DEM es un insumo fundamental para el cálculo de los factores topográficos.

3.3. Erosividad de la Lluvia (Factor R)

El índice de erosividad de la lluvia (R) cuantifica la capacidad de la lluvia para causar erosión. Consulte la literatura regional, como los trabajos de Zhang Wenbo y Fu Jinsheng (2003), para métodos de estimación basados en datos pluviométricos.

3.4. Erodibilidad del Suelo (Factor K)

El factor K representa la susceptibilidad del suelo a la erosión. Se puede estimar utilizando el modelo EPIC (Williams y Arnold, 1997).

3.4.1. Cálculo con la Calculadora Ráster de ArcGIS

Primero, calcule el factor K intermedio (K_EPIC):


(0.2 + 0.3 * Exp(-0.0256 * "HWSD1_T_SAND.tif" * (1 - "HWSD1_T_SILT.tif"/100))) * ...
Power(("HWSD1_T_SILT.tif" / ("HWSD1_T_CLAY.tif" + "HWSD1_T_SILT.tif")), 0.3) * ...
(1.0 - ((0.25 * "HWSD1_T_OC.tif") / ("HWSD1_T_OC.tif" + Exp(3.72 - 2.95 * "HWSD1_T_OC.tif")))) * ...
(1.0 - ((0.7 * (1 - "HWSD1_T_SAND.tif"/100)) / ((1 - "HWSD1_T_SAND.tif"/100) + Exp(22.9 * (1 - "HWSD1_T_SAND.tif"/100) - 5.51))))

Luego, calcule el factor K final:


(-0.01383 + 0.51575 * "Soil_Erodibility_K_epic.tif") * 0.1317

3.4.2. Cálculo con MATLAB

El siguiente script de MATLAB realiza el cálculo del factor K a partir de las propiedades del suelo (arena, limo, arcilla, carbono orgánico).


clc;
clear;

% --- Configuración de Rutas y Archivos de Entrada/Salida ---
rutaBase = 'F:\Gannan\Gannan_auxiliary\HWSD1_Gannan\TRS\';
nombreArchivoSalida = [rutaBase, 'Factor_K_Erodibilidad.TIF'];

% Rutas a los archivos GeoTIFF de propiedades del suelo
rutaArena = [rutaBase, 'HWSD1_T_SAND.tif'];
rutaLimo = [rutaBase, 'HWSD1_T_SILT.tif'];
rutaArcilla = [rutaBase, 'HWSD1_T_CLAY.tif'];
rutaCarbonoOrganico = [rutaBase, 'HWSD1_T_OC.tif'];

% --- Lectura y Preparación de Datos ---
% SAN, SIL, CLA, C representan el porcentaje de arena, limo, arcilla y
% carbono orgánico (0-30cm) respectivamente.
valorNoData = -128; 

[arena, R] = geotiffread(rutaArena);
arena(arena == valorNoData) = 0;

[limo, R] = geotiffread(rutaLimo);
limo(limo == valorNoData) = 0;

[arcilla, R] = geotiffread(rutaArcilla); 
arcilla(arcilla == valorNoData) = 0;

[carbonoOrganico, R] = geotiffread(rutaCarbonoOrganico);
carbonoOrganico(carbonoOrganico < 0) = 0;

% Convertir a tipo double para asegurar cálculos de punto flotante
arena = double(arena);
limo = double(limo);
arcilla = double(arcilla);
carbonoOrganico = double(carbonoOrganico);

% Calcular SNI (Soil Non-Sand Index, 1 - Arena/100)
indiceArenaNo = 1 - arena / 100;

% --- Cálculo del Factor K_EPIC ---
% Inicializar kFactorEpic con NaN, del mismo tamaño que las capas de suelo
kFactorEpic = NaN(size(arena)); 

% Máscara para celdas válidas (ej. donde Arena > 0 para evitar divisiones por cero implícitas)
mascaraValida = arena > 0; 

% Aplicar la fórmula del modelo EPIC para las celdas válidas
kFactorEpic(mascaraValida) = (0.2 + 0.3 .* exp(-0.0256 .* arena(mascaraValida) .* (1 - limo(mascaraValida) / 100))) ...
    .* (limo(mascaraValida) ./ (arcilla(mascaraValida) + limo(mascaraValida))).^0.3 ...
    .* (1.0 - (0.25 .* carbonoOrganico(mascaraValida)) ./ (carbonoOrganico(mascaraValida) + exp(3.72 - 2.95 .* carbonoOrganico(mascaraValida)))) ...
    .* (1.0 - (0.7 .* indiceArenaNo(mascaraValida)) ./ (indiceArenaNo(mascaraValida) + exp(22.9 .* indiceArenaNo(mascaraValida) - 5.51)));

% --- Cálculo del Factor K final ---
% Conversión de unidades o ajuste final
kFactorFinal = (-0.01383 + 0.51575 .* kFactorEpic) .* 0.1317;

% --- Escritura del Archivo GeoTIFF de Salida ---
% Especificar el código de referencia espacial si es conocido
geotiffwrite(nombreArchivoSalida, kFactorFinal, R, 'CoordRefSysCode', 32648); 
disp('Cálculo del Factor K de erodibilidad completado.');

Ambos métodos, ArcGIS y MATLAB, ofrecen resultados similares. Sin embargo, el script de MATLAB proporciona mayor flexibilidad para la automatización y el manejo de grandes volúmenes de datos.

Etiquetas: InVEST Modelado Ambiental SIG procesamiento de datos Teledetección

Publicado el 7-15 18:16

Etiquetas populares

©2026 Friki Work