Determinando área superficial a partir de un DEM con Python mediante el gradiente de la pendiente (Slope Gradient)

En esta referencia también existe una descripción para aproximar el área superficial, por píxel, de un ráster DEM mediante el gradiente de la pendiente (Slope Gradient). También está basada en la extracción de todos los valores de elevación para una ventana 3x3 para cada píxel para así determinar las áreas 3D de los 8 rectángulos resultantes.

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

Las distancias d1, d3, d5 y d7 no corresponden a la hipotenusa del triángulo rectángulo respectivo sino al tamaño del píxel ya que están referidas al círculo inscrito en el cuadrado. Lo que interesa en este caso es determinar el área de los ocho rectángulos ABEF, tal como se representan en la imagen siguiente, que debido a la pendiente han incrementado su magnitud con relación al cuadrado base ABCD.

El código completo propuesto es el 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

from osgeo import gdal, osr import numpy as np from scipy import ndimage 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: #distances d1 = d2 = d3 = d4 = d5 = d6 = d7 = d8 = xsize a1 = sqrt((subarray[0] - subarray[4])**2 + d1**2)*xsize a2 = sqrt((subarray[1] - subarray[4])**2 + d2**2)*xsize a3 = sqrt((subarray[2] - subarray[4])**2 + d3**2)*xsize a4 = sqrt((subarray[5] - subarray[4])**2 + d4**2)*xsize a5 = sqrt((subarray[8] - subarray[4])**2 + d5**2)*xsize a6 = sqrt((subarray[7] - subarray[4])**2 + d6**2)*xsize a7 = sqrt((subarray[6] - subarray[4])**2 + d7**2)*xsize a8 = sqrt((subarray[3] - subarray[4])**2 + d8**2)*xsize areas = [a1, a2, a3, a4, a5, a6, a7, a8] new_areas = [ item for item in areas if str(item) != 'nan' ] SAp = np.sum(new_areas)/len(new_areas) values.append(SAp) values = np.reshape(np.array(values), (rows,cols)) # Create gtif file driver = gdal.GetDriverByName("GTiff") output_file = "/home/zeito/pyqgis_data/surface3D_raster_slope.tif" dst_ds = driver.Create(output_file, cols, rows, 1, gdal.GDT_Float32) #writting output raster dst_ds.GetRasterBand(1).WriteArray( values ) #setting extension of output raster # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution dst_ds.SetGeoTransform(geotransform) wkt = dataset.GetProjection() # setting spatial reference of output raster srs = osr.SpatialReference() srs.ImportFromWkt(wkt) dst_ds.SetProjection( srs.ExportToWkt() ) dataset = None dst_ds = None

Al ejecutarse en la Python Console de QGIS produce el resultado siguiente; el cual es comparable al de la determinación del área superficial mediante el procedimiento SHR (surface to horizontal area ratio).