Modelamiento de Evolución de Superficie a Escala de Cuenca con Python y Landlab - Tutorial

El clima cambia, la gente cambia y la tierra también cambia con el tiempo. No podemos creer que las redes fluviales seguirán siendo las mismas durante los próximos 1000 años o que las montañas y las depresiones tendrán la misma elevación en los próximos 10000 años. Pero los cambios no están relacionados con grandes periodos de tiempo, también pueden ocurrir en décadas o años a tasas más bajas. Para evaluar esos cambios necesitamos modelar los componentes clave de la evolución del suelo: fluvial, pendiente de ladera y levantamiento. Hemos desarrollado un tutorial con Python y la librería Landlab para simular la evolución del suelo a escala de cuenca durante 100 mil años; las entradas provienen de rásters geoespaciales y los datos de salida se exportan como archivos ráster ASCII.

Necesitará un entorno conda con herramientas geoespaciales para el tutorial. Crea el entorno siguiendo este enlace:

gidahatari.com/ih-es/como-instalar-python-geopandas-en-windows-bajo-un-entorno-en-conda-tutorial

Tutorial

Código

Modeling Land Evolution at Basin Scale with Python and Landlab

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

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

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

# Open raster image 
rasterObj = rasterio.open('../data/DEM_18S_clip.tif')
show(rasterObj)

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

# define parameters for fluvial, hillslope and uplift components

uplift_rate = 0.001  # [m/yr], initially set at 0.001
K_sp = 1.0e-5  # units vary depending on m_sp and n_sp, initially set at 1e-5
m_sp = 0.5  # exponent on drainage area in stream power equation, initially 0.5
n_sp = 1.0  # exponent on slope in stream power equation, initially 1.
K_hs = 0.05  # [m^2/yr], initially 0.05

#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]))

# show number of rows, cols and resolution

print(nrows, ncols, dxy)

167 179 90.3109901477272

# define temporal distribution 

dt = 1000  # time step [yr]
total_time = 0  # amount of time the landscape has evolved [yr]
tmax = 100000  # time for the model loop to run [yr]

t = np.arange(0, tmax, dt)

# create a dataset of zero values

zr = mg.add_zeros("topographic__elevation", at="node")
imshow_grid(mg, "topographic__elevation")

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

# initialize the process components of the equation

frr = FlowAccumulator(mg)  # intializing flow routing
spr = StreamPowerEroder(
    mg, K_sp=K_sp, m_sp=m_sp, n_sp=n_sp, threshold_sp=0.0
)  # initializing stream power incision
dfn = LinearDiffuser(
    mg, linear_diffusivity=K_hs, deposit=False
)  # initializing linear diffusion
df = DepressionFinderAndRouter(mg) # Initializing the pit finder

# run the process for every time step 

for ti in t:
    zr[mg.core_nodes] += uplift_rate * dt  # uplift the landscape
    dfn.run_one_step(dt)  # diffuse the landscape
    frr.run_one_step()  # route flow
    # df.map_depressions()
    spr.run_one_step(dt)  # fluvial incision
    total_time += dt  # update time keeper
    print("\r"+str(total_time), end=' ')

# plot final surface elevation

fig = plt.figure(figsize=(18,12))

imshow_grid(
    mg, "topographic__elevation", grid_units=("m", "m"), var_name="Elevation (m)", shrink=0.5
)

title_text = f"$K_{{sp}}$={K_sp}; $K_{{hs}}$={K_hs}; $time$={total_time}; $dx$={dxy}"

plt.title(title_text)

plt.show()

# export data as geospatial rasters

from landlab.io.esri_ascii import write_esri_ascii

files = write_esri_ascii('../Out/basinLandEvolution.asc', mg)