Delimitación de Límite de Cuenca y Red Hidrica con Python y Pysheds - Tutorial

Nunca antes el proceso de delimitación de cuencas fue tan fácil. Ahora puedes delimitar cuencas y redes hídricas online con Hatari Utils:

utils.hatarilabs.com

Conoce más en este enlace:

gidahatari.com/ih-es/delimitacin-online-de-cuencas-y-redes-hdricas-con-hatari-utils-tutorial

Que pasaría si todos los procesos que realizamos en software de SIG los hicieramos en Python? Que pasaría si tratáramos a los datos espaciales como objetos y variables en un script…. enconces nos preguntaríamos si en realidad es necesario inventar de nuevo la rueda. Para qué hacer cosas de manera diferente, si lo que existe ya funciona?

La respuesta a esta interrogante es muy simple: Más control. Es tan simple y tan sencillo sobre eso.

Trabajar en Python nos da más control sobre el geoprocesamiento al dejar el entorno visual de clicks sobre iconos. Con Python en un entorno de Jupyter Notebook podemos referenciar los archivos a importar, definir los geoprocesos y sus opciones, hacer representaciones de los datos espaciales intermedios y finales, y exportar los resultados en formatos compatibles con cualquier entorno de SIG. Existen otras ventajas del análisis espacial en Python como son la reproducibilidad y la velocidad de cálculo.

Este enfoque de cálculo es nuevo entre los profesionales y solo algunos hydro-geeks se sentirán motivados en introducirse en este tema. Lo que si creemos es que varios procesos y modelos hidrológicos se desarrollarán en Python a partir de las herramientas disponibles ahora.

El paquete Pysheds es la librería que te permite la delimitación de cuencas y la extracción de la red hídrica en Python 3. Esta librería requiere el uso de otras librerías avanzadas de análisis espacial como Numpy, Pandas, Scipy, Scikit-image, Rasterio entre otras. Este tutorial muestra el procedimiento completo dentro de un Jupyter Notebook para:

la importación de un modelo digital de elevación corregido (hidrológicamente correcto)
el cálculo de raster de dirección,
la delimitación de cuencas,
el análisis de acumulación de flujo
la determinación de la red hídrica
la generación de archivos espaciales raster y vectoriales

Dado que el tutorial se realiza en entorno Windows, también se ha realizado un tutorial para la instalación de este paquete y los paquetes complementarios.

Tutorial

Tutorial para la instalación de Rasterio, GDAL, Geopandas y Pysheds para Anaconda en Windows

Tutorial para la delimitación de límite de cuenca y red hidrica con Python y Pysheds

Código en Python

Este es el código completo en Python 3 para la realización de este análisis espacial

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import geopandas as gpd
from pysheds.grid import Grid
import mplleaflet
%matplotlib inline

grid = Grid.from_raster('../Rst/LocalDem.tif', data_name='dem')

def plotFigure(data, label, cmap='Blues'):
    plt.figure(figsize=(12,10))
    plt.imshow(data, extent=grid.extent, cmap=cmap)
    plt.colorbar(label=label)
    plt.grid()

plotFigure(grid.dem, 'Elevation (m)')

#N    NE    E    SE    S    SW    W    NW
dirmap = (64,  128,  1,   2,    4,   8,    16,  32)

grid.flowdir(data='dem', out_name='dir', dirmap=dirmap)

plotFigure(grid.dir,'Flow Directiom','viridis')

# Specify pour point
x, y = 285612.017,2936416.682

# Delineate the catchment
grid.catchment(data='dir', x=x, y=y, dirmap=dirmap, out_name='catch',
               recursionlimit=15000, xytype='label', nodata_out=0)

# Clip the bounding box to the catchment
grid.clip_to('catch')

# Get a view of the catchment
demView = grid.view('dem', nodata=np.nan)

plotFigure(demView,'Elevation')

#export selected raster
grid.to_raster(demView, '../Output/clippedElevations.tif')

grid.accumulation(data='catch', dirmap=dirmap, pad_inplace=False, out_name='acc')

accView = grid.view('acc', nodata=np.nan)
plotFigure(accView,"Cell Number",'PuRd')

streams = grid.extract_river_network('catch', 'acc', threshold=200, dirmap=dirmap)
streams["features"][:2]

def saveDict(dic,file):
    f = open(file,'w')
    f.write(str(dic))
    f.close()

#save geojson as separate file
saveDict(streams,'../Output/streams.geojson')

streamNet = gpd.read_file('../Output/streams.geojson')
streamNet.crs = {'init' :'epsg:32613'}

# The polygonize argument defaults to the grid mask when no arguments are supplied
shapes = grid.polygonize()

# Plot catchment boundaries
fig, ax = plt.subplots(figsize=(6.5, 6.5))

for shape in shapes:
    coords = np.asarray(shape[0]['coordinates'][0])
    ax.plot(coords[:,0], coords[:,1], color='cyan')

ax.set_xlim(grid.bbox[0], grid.bbox[2])
ax.set_ylim(grid.bbox[1], grid.bbox[3])
ax.set_title('Catchment boundary (vector)')
gpd.plotting.plot_dataframe(streamNet, None, cmap='Blues', ax=ax)

#ax = streamNet.plot()
mplleaflet.display(fig=ax.figure, crs=streamNet.crs, tiles='esri_aerial')

# The polygonize argument defaults to the grid mask when no arguments are supplied
shapes = grid.polygonize()

# Plot catchment boundaries
fig, ax = plt.subplots(figsize=(6.5, 6.5))

for shape in shapes:
    coords = np.asarray(shape[0]['coordinates'][0])
    ax.plot(coords[:,0], coords[:,1], color='cyan')

ax.set_xlim(grid.bbox[0], grid.bbox[2])
ax.set_ylim(grid.bbox[1], grid.bbox[3])
ax.set_title('Catchment boundary (vector)')
gpd.plotting.plot_dataframe(streamNet, None, cmap='Blues', ax=ax)

#ax = streamNet.plot()
mplleaflet.display(fig=ax.figure, crs=streamNet.crs, tiles='esri_aerial')