Cálculo de la dispersión de mayor pendiente con área de drenaje con Python y Landlab - Tutorial

El análisis de superficies a nivel regional requiere un mayor grado de cálculo para evaluar varios temas de gestión del terreno. Por lo general, esos tipos de análisis sobre archivos ráster se limitaban al software de escritorio GIS con resultados para un solo ráster y con opciones limitadas para crear gráficos de valores de píxeles.

Python, Rasterio, Landlab y otros paquetes de alto rendimiento nos permiten tener otra vista de la elevación de nuestra superficie y evaluar diferentes variables que eran impensables con el software GIS tradicional. Este tutorial muestra el procedimiento completo para abrir un archivo raster con Rasterio, importar como model grid de Landlab, llenar sumideros y calcular el área de drenaje. Finalmente, a partir de los resultados de Landlab, se construyó un Pandas dataframe y las áreas de drenaje se clasificaron en intervalos para trazarse en un diagrama de caja con la pendiente más pronunciada.

Tutorial

Datos de entrada

Puedes descargar los datos de entrada desde este enlace.

Código

# import generic packages
import numpy as np
from matplotlib import pyplot as plt
import pandas as pd

# import geospatial packages
import rasterio
from rasterio.plot import show

# import landlab components
from landlab import RasterModelGrid, imshow_grid
from landlab.components import (
    DepressionFinderAndRouter,
    FlowAccumulator
)

#open raster image 
rasterObj = rasterio.open('../rst/clippedDem_12N.tif')
show(rasterObj)

<Axes: >

#extract array from raster
elevArray = rasterObj.read(1)
plt.imshow(elevArray)

<matplotlib.image.AxesImage at 0x1b7f41183a0>

#show array
elevArray

array([[1475, 1470, 1467, ..., 1809, 1809, 1808],
       [1469, 1464, 1460, ..., 1816, 1816, 1816],
       [1461, 1458, 1453, ..., 1815, 1819, 1820],
       ..., 
       [1618, 1616, 1616, ..., 1752, 1753, 1754],
       [1621, 1617, 1617, ..., 1743, 1745, 1749],
       [1621, 1617, 1617, ..., 1739, 1740, 1743]], dtype=int16)

#create grid from raster attributes
nrows = rasterObj.height  # number of raster cells on each side, initially 150
ncols = rasterObj.width
dxy = rasterObj.transform.a  # side length of a raster model cell, or resolution [m], initially 50

# define a landlab raster
mg = RasterModelGrid(shape=(nrows, ncols), 
                     xy_spacing=dxy,
                     xy_of_lower_left=(rasterObj.bounds[0],rasterObj.bounds[1]))

# create a dataset of zero values
zr = mg.add_zeros("topographic__elevation", at="node")
imshow_grid(mg, "topographic__elevation", shrink=0.5)

# apply cell elevation to defined arrray
zr += elevArray[::-1,:].ravel()
imshow_grid(mg, "topographic__elevation", shrink=0.5)

#run the component to accumulate flow and calculate area
fr = FlowAccumulator(mg, flow_director='D8')
fr.run_one_step()

#show the partial drainage area and check uncontinours streams
imshow_grid(mg, 'drainage_area', shrink=0.5)

#fill sinks
df = DepressionFinderAndRouter(mg)
df.map_depressions()

#plot corrected drainage area
fig = plt.figure(figsize=(12,6))
imshow_grid(mg, 'drainage_area', shrink=0.5)

#create a pandas dataframe for boxplot
dfMg = pd.DataFrame(columns=['slope','area(Ha)'])
dfMg['slope']=mg.at_node["topographic__steepest_slope"][mg.core_nodes]
dfMg['area(Ha)']=mg.at_node["drainage_area"][mg.core_nodes]/10000
dfMg.head()

	slope	area(Ha)
0	0.038986	0.065793
1	0.303241	0.065793
2	0.545807	0.065793
3	0.467834	0.065793
4	0.389862	0.065793

#create interval an apply into a new columns
areaInterval = pd.interval_range(0,dfMg['area(Ha)'].max()//100*100,25)
dfMg['areaInterval(Ha)'] = pd.cut(dfMg['area(Ha)'], bins=areaInterval)
dfMg.head()

	slope	area(Ha)	areaInterval(Ha)
0	0.038986	0.065793	(0.0, 580.0]
1	0.303241	0.065793	(0.0, 580.0]
2	0.545807	0.065793	(0.0, 580.0]
3	0.467834	0.065793	(0.0, 580.0]
4	0.389862	0.065793	(0.0, 580.0]

#plot the graph of maximum slope vs drainage area 
fig, ax = plt.subplots(figsize=(12, 6))
dfMg.boxplot('slope', by='areaInterval(Ha)', ax=ax, rot=90)
#ax.set_yscale('log')
ax.set_ylim(-0.01,0.8)