基于Landlab的雨洪分析，python语言

Landlab简介

来自Landlab官网的介绍

Landlab is a Python-based modeling environment that allows scientists and students to build numerical landscape models. Designed for disciplines that quantify earth surface dynamics such as geomorphology, hydrology, vegetation ecology, glaciology, and stratigraphy, it can also be used for any application that needs 2D grid-based numerical models.

Landlab provides components to compute flows (such as water, sediment, glacial ice, volcanic material, or landslide debris) across a gridded terrain. With its robust, reusable components, Landlab allows scientists to quickly build landscape model experiments and compute mass balance across scales.

Landlab库代码地址

https://github.com/landlab/landlab

安装和用户手册

https://landlab.github.io/#/

用于雨洪仿真的缺点

基于博主的使用效果，其主要缺点在于：

1. CPU计算，未使用GPU加速。如果你使用的DEM GSD比较高的情况下，比如1m甚至cm级别，会导致计算量非常大。例如1024*1024像素的数据，如果是城区，可能会需要几个小时甚至十几小时来处理，如果是海边可能需要2小时来处理。基于此，建议用户采用较低分辨率的DEM数据比如30m来看整体的情况。如果需要GPU加速的库来实现雨洪仿真，可以关注我的下一篇博客。
2. 排水网和排水边界的设置不够便捷，可以查看官方手册和论文，对此有说明，但是使用时确实不够便利。下一篇博客中另一个可GPU加速的库对排水网的设置就比较便捷。
3. 边界排水问题。这个问题并不是Landlab本身问题，而是我在使用浅水方程处理洪水时都面临的问题，如果有更好解决方案的兄弟欢迎给我留言。这个问题是说DEM是被分成一块一块来进行计算的，这导致你对图片边界排水能力的定义是有问题的，例如一个凹地被左右两景DEM数据切分，那左图的右边界和右图的坐边界排水能力如何计算呢？如果定义每小时排XX的水，那这样子并不符合凹的的定义，但是简单将图片边界定义为不排水，那仿真结果一定是边界区域存在一定程度的积水，当把多张DEM分析结果和在一起放在地图上时，会看到明显的每块DEM的边界存在积水，这并不符合真实的世界情况。

本文目的

Landlab是非常强大的地标动力学仿真库。本文只是使用其中的雨洪仿真(Rainfall Flood Simulation)相关函数（基于浅水方程），实现的效果就是用户输入降雨强度，降雨时长和排水网信息，基于已有的DEM给出不同区域积水情况。
Landlab可以做的事情远不止Rainfall Flood Simulation这点小事，地表径流仿真，水库仿真，沉积仿真，等等也不局限于水。

本文代码及解释

代码中task是一个json方便API调用的，其定义为：
{‘watershed_dem’: ‘dem.tif’, ‘starting_precip_mmhr’: 142.0, ‘storm_duration’: 60.0, ‘elapsed_time’: 0.0, ‘model_run_time’: 60.0, ‘GPU’: 0}
这里各自的格式要求和单位为：

• watershed_dem: DEM的地址，tif、tiff、asc都可以，具体你可以在下文代码中的ext_list来做格式限制
• starting_precip_mmhr: 降雨强度，单位是 毫米/小时
• storm_duration: 降雨时长，单位是 秒
• elapsed_time: 可被认为是开始的时刻 单位是秒
• model_run_time: 整个仿真运行的时间，不小于storm_duration，单位是秒。这里其实是整个事件的时间，比如说24小时，其中只有1个小时下雨了，因此 model_run_time >= storm_duration
• GPU: 这个参数只有0和1，只是为了配合上一篇博客使用的，在下面代码中无用

下述代码如果和我的上一个博客一起使用的话，那么不要考虑try catch的问题，因为上一个博客已经考虑了错误处理，这里处理错误的话反而没有很好将报错信息传递给celery

下面代码中

import os
import time
import rasterio
import numpy as np


from tqdm import tqdm
from landlab import RasterModelGrid
from landlab.io import esri_ascii, write_esri_ascii
from landlab.components import FlowAccumulator, OverlandFlow


def landlab_rainfall_simulation_cpu(task, output_dir='/home/xxx/data/rainfall_simulation/'):


    process_time = {
   }

    ext_list = ['.tiff', ".tif", ".asc"]


    try:
        watershed_dem = task['watershed_dem']

        # the intensity of starting_precip_mmhr mm/hr, with a duration of storm_duration seconds.
        starting_precip_mmhr = task['starting_precip_mmhr']
        storm_duration = task['storm_duration']

        elapsed_time = task['elapsed_time']  # s
        model_run_time = task['model_run_time']  # s

        assert os.path.exists(watershed_dem