Determinando área superficial a partir de un DEM con Python mediante el procedimiento SHR (surface to horizontal area ratio)

En esta referencia existe una descripción más o menos pormenorizada de cómo aproximar el área superficial, por píxel, de un ráster DEM mediante el procedimiento SHR (surface to horizontal area ratio). Esta basado en la extracción de todos los valores de elevación para una ventana 3x3 para cada píxel y así determinar las áreas 3D de los 8 triángulos resultantes.

Esto puede observarse con más detalle (vista cenital) en la imagen siguiente:

Para facilitar la codificación, los triángulos están numerados entre paréntesis, los vértices con números enteros y las distancias en el plano del ráster precedidas con una d. Sin embargo, las distancias de los tríangulos que nos interesan son las que aparecen en la imagen siguiente y que pueden ser determinadas fácilmente por consideraciones eminentemente geométricas empleando las distancias referidas arriba y los valores de altitud de los 9 píxeles involucrados para cada ventana 3x3.

El área de los trángulos escalenos mostrados arriba puede ser calculada mediante la bien conocida fórmula de Herón. También puede observarse que el área estimada pertenece a 4 píxeles por lo que, a efectos de normalización, ésta debe finalmente dividirse entre cuatro para corresponder a un sólo píxel. Sin embargo, lo anterior es verdad cuando no existen problemas de borde. En estos últimos los valores de normalización serán menores (2 ó 1).

Para acceder a los valores de las ventanas 3x3 del ráster DEM se empleó el método 'generic_filter' del módulo python ndimage de scipy. En el código siguiente se puede observar que también se consideran los efectos de borde. Emplea un pequeño ráster de 9x9 para facilitar la depuración del código.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

from osgeo import gdal, osr import numpy as np from scipy import ndimage from os import path, remove from math import sqrt def func(x, subarrays): subarrays.append(x) return 1 dataset = gdal.Open("/home/zeito/pyqgis_data/dem_part_raster.tif") geotransform = dataset.GetGeoTransform() xsize = geotransform[1] ysize = -geotransform[5] data = dataset.ReadAsArray() subarrays = [] result = ndimage.generic_filter(data, func, size = 3, mode='constant', cval = np.NaN, extra_arguments = (subarrays,)) cols = dataset.RasterXSize rows = dataset.RasterYSize values = [] for subarray in subarrays: print subarray dataset = None

Cuando el código anterior se ejecuta en la Python Console de QGIS se obtiene el resultado siguiente:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82

>>>execfile(u'/home/zeito/pyqgis_scripts/calculating_surface_area3.py'.encode('UTF-8')) [ nan nan nan nan 2703. 2781. nan 2682. 2767.] [ nan nan nan 2703. 2781. 2836. 2682. 2767. 2819.] [ nan nan nan 2781. 2836. 2894. 2767. 2819. 2863.] [ nan nan nan 2836. 2894. 2928. 2819. 2863. 2887.] [ nan nan nan 2894. 2928. 2966. 2863. 2887. 2916.] [ nan nan nan 2928. 2966. 3002. 2887. 2916. 2950.] [ nan nan nan 2966. 3002. 3037. 2916. 2950. 3004.] [ nan nan nan 3002. 3037. 3078. 2950. 3004. 3033.] [ nan nan nan 3037. 3078. nan 3004. 3033. nan] [ nan 2703. 2781. nan 2682. 2767. nan 2652. 2729.] [ 2703. 2781. 2836. 2682. 2767. 2819. 2652. 2729. 2776.] [ 2781. 2836. 2894. 2767. 2819. 2863. 2729. 2776. 2809.] [ 2836. 2894. 2928. 2819. 2863. 2887. 2776. 2809. 2840.] [ 2894. 2928. 2966. 2863. 2887. 2916. 2809. 2840. 2880.] [ 2928. 2966. 3002. 2887. 2916. 2950. 2840. 2880. 2932.] [ 2966. 3002. 3037. 2916. 2950. 3004. 2880. 2932. 2963.] [ 3002. 3037. 3078. 2950. 3004. 3033. 2932. 2963. 2983.] [ 3037. 3078. nan 3004. 3033. nan 2963. 2983. nan] [ nan 2682. 2767. nan 2652. 2729. nan 2671. 2703.] [ 2682. 2767. 2819. 2652. 2729. 2776. 2671. 2703. 2738.] [ 2767. 2819. 2863. 2729. 2776. 2809. 2703. 2738. 2775.] [ 2819. 2863. 2887. 2776. 2809. 2840. 2738. 2775. 2820.] [ 2863. 2887. 2916. 2809. 2840. 2880. 2775. 2820. 2858.] [ 2887. 2916. 2950. 2840. 2880. 2932. 2820. 2858. 2878.] [ 2916. 2950. 3004. 2880. 2932. 2963. 2858. 2878. 2919.] [ 2950. 3004. 3033. 2932. 2963. 2983. 2878. 2919. 2952.] [ 3004. 3033. nan 2963. 2983. nan 2919. 2952. nan] [ nan 2652. 2729. nan 2671. 2703. nan 2671. 2703.] [ 2652. 2729. 2776. 2671. 2703. 2738. 2671. 2703. 2738.] [ 2729. 2776. 2809. 2703. 2738. 2775. 2703. 2738. 2775.] [ 2776. 2809. 2840. 2738. 2775. 2820. 2738. 2775. 2820.] [ 2809. 2840. 2880. 2775. 2820. 2858. 2775. 2820. 2858.] [ 2840. 2880. 2932. 2820. 2858. 2878. 2820. 2858. 2878.] [ 2880. 2932. 2963. 2858. 2878. 2919. 2858. 2878. 2919.] [ 2932. 2963. 2983. 2878. 2919. 2952. 2878. 2919. 2952.] [ 2963. 2983. nan 2919. 2952. nan 2919. 2952. nan] [ nan 2671. 2703. nan 2671. 2703. nan 2644. 2670.] [ 2671. 2703. 2738. 2671. 2703. 2738. 2644. 2670. 2704.] [ 2703. 2738. 2775. 2703. 2738. 2775. 2670. 2704. 2747.] [ 2738. 2775. 2820. 2738. 2775. 2820. 2704. 2747. 2784.] [ 2775. 2820. 2858. 2775. 2820. 2858. 2747. 2784. 2812.] [ 2820. 2858. 2878. 2820. 2858. 2878. 2784. 2812. 2857.] [ 2858. 2878. 2919. 2858. 2878. 2919. 2812. 2857. 2901.] [ 2878. 2919. 2952. 2878. 2919. 2952. 2857. 2901. 2935.] [ 2919. 2952. nan 2919. 2952. nan 2901. 2935. nan] [ nan 2671. 2703. nan 2644. 2670. nan 2616. 2654.] [ 2671. 2703. 2738. 2644. 2670. 2704. 2616. 2654. 2687.] [ 2703. 2738. 2775. 2670. 2704. 2747. 2654. 2687. 2717.] [ 2738. 2775. 2820. 2704. 2747. 2784. 2687. 2717. 2751.] [ 2775. 2820. 2858. 2747. 2784. 2812. 2717. 2751. 2788.] [ 2820. 2858. 2878. 2784. 2812. 2857. 2751. 2788. 2829.] [ 2858. 2878. 2919. 2812. 2857. 2901. 2788. 2829. 2873.] [ 2878. 2919. 2952. 2857. 2901. 2935. 2829. 2873. 2907.] [ 2919. 2952. nan 2901. 2935. nan 2873. 2907. nan] [ nan 2644. 2670. nan 2616. 2654. nan 2633. 2672.] [ 2644. 2670. 2704. 2616. 2654. 2687. 2633. 2672. 2696.] [ 2670. 2704. 2747. 2654. 2687. 2717. 2672. 2696. 2725.] [ 2704. 2747. 2784. 2687. 2717. 2751. 2696. 2725. 2751.] [ 2747. 2784. 2812. 2717. 2751. 2788. 2725. 2751. 2780.] [ 2784. 2812. 2857. 2751. 2788. 2829. 2751. 2780. 2811.] [ 2812. 2857. 2901. 2788. 2829. 2873. 2780. 2811. 2855.] [ 2857. 2901. 2935. 2829. 2873. 2907. 2811. 2855. 2896.] [ 2901. 2935. nan 2873. 2907. nan 2855. 2896. nan] [ nan 2616. 2654. nan 2633. 2672. nan 2636. 2672.] [ 2616. 2654. 2687. 2633. 2672. 2696. 2636. 2672. 2697.] [ 2654. 2687. 2717. 2672. 2696. 2725. 2672. 2697. 2723.] [ 2687. 2717. 2751. 2696. 2725. 2751. 2697. 2723. 2749.] [ 2717. 2751. 2788. 2725. 2751. 2780. 2723. 2749. 2782.] [ 2751. 2788. 2829. 2751. 2780. 2811. 2749. 2782. 2812.] [ 2788. 2829. 2873. 2780. 2811. 2855. 2782. 2812. 2843.] [ 2829. 2873. 2907. 2811. 2855. 2896. 2812. 2843. 2866.] [ 2873. 2907. nan 2855. 2896. nan 2843. 2866. nan] [ nan 2633. 2672. nan 2636. 2672. nan nan nan] [ 2633. 2672. 2696. 2636. 2672. 2697. nan nan nan] [ 2672. 2696. 2725. 2672. 2697. 2723. nan nan nan] [ 2696. 2725. 2751. 2697. 2723. 2749. nan nan nan] [ 2725. 2751. 2780. 2723. 2749. 2782. nan nan nan] [ 2751. 2780. 2811. 2749. 2782. 2812. nan nan nan] [ 2780. 2811. 2855. 2782. 2812. 2843. nan nan nan] [ 2811. 2855. 2896. 2812. 2843. 2866. nan nan nan] [ 2855. 2896. nan 2843. 2866. nan nan nan nan]

Corresponde a los 81 sub arreglos correspondientes. Los problemas de borde se manifiestan con la impresión de 'nan'.

El código anterior sirvió de base para la determinación de las áreas superficiales por píxel. Éstas se almacenaron como ráster empleando los recursos de GDAL. El resultado puede observarse en la imagen siguiente: