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

原创已于 2025-01-15 17:42:38 修改 · 1.3k 阅读

12 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

收录于

python

Simulation

Flooding

于 2025-01-15 15:45:22 首次发布

Python3.8

Python 是一种高级、解释型、通用的编程语言，以其简洁易读的语法而闻名，适用于广泛的应用，包括Web开发、数据分析、人工智能和自动化脚本

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/#/

用于雨洪仿真的缺点

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

CPU计算，未使用GPU加速。如果你使用的DEM GSD比较高的情况下，比如1m甚至cm级别，会导致计算量非常大。例如1024*1024像素的数据，如果是城区，可能会需要几个小时甚至十几小时来处理，如果是海边可能需要2小时来处理。基于此，建议用户采用较低分辨率的DEM数据比如30m来看整体的情况。如果需要GPU加速的库来实现雨洪仿真，可以关注我的下一篇博客。
排水网和排水边界的设置不够便捷，可以查看官方手册和论文，对此有说明，但是使用时确实不够便利。下一篇博客中另一个可GPU加速的库对排水网的设置就比较便捷。
边界排水问题。这个问题并不是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), "input cannot be None"

        ext = os.path.splitext(watershed_dem)[1]

        dem_file_name = os.path.split(watershed_dem)[-1]

        output_dir = os.path.join(output_dir, os.path.splitext(dem_file_name)[0])

        if not os.path.exists(output_dir):
            os.mkdir(output_dir)

        output_flood_map = (output_dir + os.path.splitext(dem_file_name)[0] + '-result-' + str(task['starting_precip_mmhr'])
                            + "-" + str(task['storm_duration']) + "-" + str(task['elapsed_time']) + "-" + str(
                    task['model_run_time']) + ext)
        print('output_flood_map: ' + output_flood_map)

        assert ext in ext_list, "please check the input format, it can only be " + str(ext_list)
        if ext == ".asc":
            with open(watershed_dem) as fp:
                rmg = esri_ascii.load(fp, name="topographic__elevation", at="node")

            
            height = rmg.shape[0]
            width = rmg.shape[1]
            count = 1
            crs = 1
            transform = 2

        else:
            print('read image')
            time1 = time.time()
            rasterObj = rasterio.open(watershed_dem)
            elevArray = rasterObj.read(1)
            dxy = (rasterObj.transform.a,
                   -rasterObj.transform.e)  # side length of a raster model cell, or resolution [m], initially 50
            time2 = time.time()
            process_time['read_image'] = time2 - time1

            height = rasterObj.height
            width = rasterObj.width
            count = 1  # number of bands
            crs = rasterObj.crs
            transform = rasterObj.transform

            # define a landlab raster
            print('construct RasterModelGrid')
            time1 = time.time()
            rmg = RasterModelGrid(shape=(height, width),
                                  xy_spacing=dxy,
                                  xy_of_lower_left=(rasterObj.bounds[0], rasterObj.bounds[1]))

            # create a dataset of zero values
            zr = rmg.add_zeros("topographic__elevation", at="node")

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

            # clear empty values
            rmg.set_nodata_nodes_to_closed(zr, -9999)
            time2 = time.time()
            process_time['construct_RasterModelGrid'] = time2 - time1

        starting_precip_ms = starting_precip_mmhr * (2.77778 * 10 ** -7)

        ## The Flow Router calculates drainage area, which is helpful for
        ## calculating equilibrium discharge, which we illustrate later.
        print('start flowaccumulator calculation')
        time1 = time.time()
        fr = FlowAccumulator(rmg)  # Instantiate flow router
        fr.run_one_step()  # Drainage area calculated
        time2 = time.time()
        process_time['flowaccumulator'] = time2 - time1

        # ## Set boundary conditions on the grid
        # print('set boundary conditions on the grid')
        # time1 = time.time()
        # rmg.set_watershed_boundary_condition(rmg.at_node["topographic__elevation"])
        # time2 = time.time()
        # process_time['set_boundary_conditions'] = time2 - time1

        ## instantiate OverlandFlow object
        print('add zero surface_water__depth')
        rmg.add_zeros("surface_water__depth", at="node")
        print('start OverlandFlow calculation')
        time1 = time.time()
        of = OverlandFlow(rmg, alpha=0.45, steep_slopes=True)
        time2 = time.time()
        process_time['OverlandFlow'] = time2 - time1

        hydrograph_time = []

        ## Setting initial fields...
        rmg.at_node["surface_water__discharge"] = np.zeros(rmg.number_of_nodes)

        print('start the surface water calculation')
        time1 = time.time()
        with tqdm(total=model_run_time) as pbar:
            pbar.set_description(dem_file_name)
            while elapsed_time < model_run_time:
                # Setting the adaptive time step
                of.dt = of.calc_time_step()

                ## The storm starts when the model starts. While the elapsed time is less
                ## than the storm duration, we add water to the system as rainfall.
                if elapsed_time < (storm_duration):
                    of.rainfall_intensity = starting_precip_ms
                else:  # elapsed time exceeds the storm duration, rainfall ceases.
                    of.rainfall_intensity = 0.0

                of.run_one_step()  # Generating overland flow based on the deAlmeida solution.

                ## Append time and discharge to their lists to save data and for plotting.
                hydrograph_time.append(elapsed_time)
                q = rmg.at_link["surface_water__discharge"]

                elapsed_time += of.dt

                pbar.update(of.dt)
        time2 = time.time()
        process_time['surface_water'] = time2 - time1

        print('saving......')
        time1 = time.time()
        if ext == '.asc':
            write_esri_ascii(output_flood_map, rmg)
        else:
            # convert the resulting data to a numpy array
            zArray = rmg.at_node["surface_water__depth"].reshape((height, width))[::-1, :]

            # write the array to a new GeoTIFF file
            with rasterio.open(
                    output_flood_map,
                    'w',
                    driver='GTiff',
                    height=height,
                    width=width,
                    count=count,  # number of bands
                    dtype=zArray.dtype,
                    crs=crs,
                    transform=transform
            ) as dst:
                dst.write(zArray, 1)  # write the data to the first band
        time2 = time.time()
        process_time['saving'] = time2 - time1

        for key, value in process_time.items():
            print(key, "=", value)

        return True, output_flood_map
    except Exception as exc:
        return False, exc