Notebook for A Vector Approach to Drainage Network Analysis Based on LiDAR Data

$Fangzheng$ $Lyu^{1,2}$, $Zewei$ $Xu^{1,2}$, $Xinlin$ $Ma^{3}$, $Shaohua$ $Wang^{1,2}$, $Zhiyu$ $Li^{1,2}$, $Shaowen$ $Wang^{1,2}$

$^{1}$$Department$ $of$ $Geography$ $and$ $Geographic$ $Information$ $Science$, $University$ $of$ $Illinois$ $at$ $Urbana-Champaign$, $Urbana$, $IL$, $USA$
$^{2}$$CyberGIS$ $Center$ $for$ $Advanced$ $Digital$ $and$ $Spatial$ $Studies$, $University$ $of$ $Illinois$ $at$ $Urbana-Champaign$, $Urbana$, $IL$, $USA$
$^{3}$$Department$ $of$ $Management$ $and$ $Urban$ $Planning$, $Tsinghua$ $University$, $Beijing$, $China$
$Corresponding$ $Author:$ $shaowen@illinois.edu$

Drainage network analysis is fundamental to understanding the characteristics of surface hydrology. Based on elevation data, drainage network analysis is often used to extract key hydrological features like drainage networks and streamlines. Limited by raster-based data models, conventional drainage network algorithms typically allow water to flow in 4 or 8 directions (surrounding grids) from a raster grid. To resolve this limitation, this paper describes a new vector-based method for drainage network analysis that allows water to flow in any direction around each location. The method is enabled by rapid advances in Light Detection and Ranging (LiDAR) remote sensing and high-performance computing. The drainage network analysis is conducted using a high-density point cloud instead of Digital Elevation Models (DEMs) at coarse resolutions. Our computational experiments show that the vector-based method can better capture water flows without limiting the number of directions due to imprecise DEMs. Our case study applies the method to Rowan County watershed, North Carolina in the US. After comparing the drainage networks and streamlines detected with corresponding reference data from US Geological Survey generated from the Geonet software, we find that the new method performs well in capturing the characteristics of water flows on landscape surfaces in order to form an accurate drainage network.

This notebook is a sample notebook for running a small size dataset (from the LiDAR dataset of Rowan Watershed) with our new method for drainage system analysis. The algorithm is implemented with an execuation sample dataset.

Paper:

Fangzheng Lyu, Zewei Xu, Xinlin Ma, Shaohua Wang, Zhiyu Li, Shaowen Wang, A vector-based method for drainage network analysis based on LiDAR data, Computers & Geosciences, Volume 156, 2021, 104892, ISSN 0098-3004

https://www.sciencedirect.com/science/article/pii/S0098300421001849

Github repo:

https://github.com/cybergis/Drainage-System-Analysis-Research

Psudocode

Set up

The Library used for in the algorithm is set up here

Set Up Environment and Import Library

In [1]:
%load_ext pycodestyle_magic
#%pycodestyle_on
%flake8_on --ignore E501,E251,E231,E225

Import library

In [2]:
import os
import datetime
import laspy
import numpy as np
import pandas as pd
from numpy.linalg import inv
import math
import sys

Set up input variables

In [3]:
# var1 is the starting location in the x-axis
var1 = 5
# var2 is the starting location in the y-axis
var2 = 5
# var3 is the angle different we want
var3 = 30

Prepare sample LiDAR data

In [4]:
import os
from shutil import copyfile
import gdown

# download sample data if it is not local
local_sample_path = "./3_2.las"
if not os.path.isfile(local_sample_path):
    print("Missing sample data!")
    sample_shared = "/home/jovyan/shared_data/data/drainage_system_analysis_research/3_2.las"
    if os.path.join(sample_shared):
        print("Copying sample data from shared folder...")
        copyfile(sample_shared, local_sample_path)
    else:
        print("Copying sample data from shared google drive...")
        url_gdrive = 'https://drive.google.com/uc?id=1JOl1IylIZg72QdMM89xs10NLFk4rObml'
        gdown.download(url_gdrive, local_sample_path, quiet=False)
else:
    print("Sample data already exists")
if not os.path.isfile(local_sample_path):
    print("Can not retrieve sample data!")
Sample data already exists
In [5]:
!ls {local_sample_path} -alh
-rw-r--r-- 1 jovyan users 159M Sep  9 15:25 ./3_2.las

Read LiDAR data using laspy

In [6]:
# Store the LiDAR data as infile
infile = laspy.file.File(local_sample_path, mode="r")
# Get the value for x axis, y axis, and elevation for LiDAR point cloud
ground_x = infile.x
ground_y = infile.y
ground_z = infile.z
# Normalize the output for the LiDAR file
ground_x_2 = ground_x-ground_x.min()
ground_y_2 = ground_y-ground_y.min()

The total amount of LiDAR point in the LiDAR file

In [7]:
len(ground_z)
Out[7]:
5546582

The x, y, z value for points in the LiDAR file are stored

In [8]:
x = ground_x_2
y = ground_y_2
z = ground_z

Use a dictionary to store the datset

In [9]:
threedarray = np.vstack((x,y,z)).T
dictionary = pd.Series(threedarray.tolist(), index=map(lambda a: round(a,2), x.tolist()))

An example of the format of the dataset

In [10]:
dictionary[1.11]
Out[10]:
1.11     [1.1100000001024455, 238.38000000000466, 690.64]
1.11     [1.1100000001024455, 357.63000000000466, 691.02]
1.11     [1.1100000001024455, 493.71999999997206, 695.84]
1.11     [1.1100000001024455, 510.03000000002794, 695.94]
1.11    [1.1100000001024455, 543.7900000000373, 694.55...
1.11      [1.1100000001024455, 608.5899999999674, 694.42]
1.11      [1.1100000001024455, 668.1500000000233, 698.37]
1.11      [1.1100000001024455, 804.5600000000559, 698.38]
1.11      [1.1100000001024455, 803.4599999999627, 698.38]
1.11      [1.1100000001024455, 795.9000000000233, 698.37]
1.11      [1.1100000001024455, 829.4000000000233, 698.32]
1.11      [1.1100000001024455, 843.5400000000373, 697.71]
1.11       [1.1100000001024455, 850.109999999986, 697.59]
1.11                 [1.1100000001024455, 899.75, 699.48]
1.11      [1.1100000001024455, 927.9899999999907, 699.64]
1.11      [1.1100000001024455, 1039.5100000000093, 699.7]
1.11     [1.1100000001024455, 1207.2900000000373, 703.09]
1.11                [1.1100000001024455, 1221.25, 704.29]
1.11     [1.1100000001024455, 1261.1600000000326, 707.38]
1.11      [1.1100000001024455, 1298.140000000014, 707.98]
1.11      [1.1100000001024455, 1284.0899999999674, 707.6]
1.11      [1.1100000001024455, 1332.920000000042, 708.77]
1.11      [1.1100000001024455, 1338.140000000014, 708.79]
1.11      [1.1100000001024455, 1457.170000000042, 711.67]
1.11      [1.1100000001024455, 1442.109999999986, 711.32]
1.11      [1.1100000001024455, 1463.6600000000326, 711.9]
1.11    [1.1100000001024455, 1756.0500000000466, 712.6...
dtype: object

Create a hash table data structure to store the LiDAR data

Create the index for each LiDAR data point

In [11]:
list_a = map(lambda a: 10000*int(a), x.tolist())
list_b = map(lambda a: int(a), y.tolist())
index_list = [sum(x) for x in zip(list_a, list_b)]

Generate the hash table and use dictionary data structure to store the dataset

In [12]:
grid_dictionary = pd.Series(threedarray.tolist(), index=index_list)

Get a overview of the dataset

In [13]:
grid_dictionary
Out[13]:
16        [0.02000000001862645, 16.46999999997206, 688.92]
6              [0.0, 6.800000000046566, 688.9300000000001]
28                       [0.0, 28.400000000023283, 688.96]
280008     [28.010000000009313, 8.900000000023283, 689.47]
60002       [6.370000000111759, 2.2800000000279397, 689.0]
                                ...                       
352542      [35.6500000001397, 2542.0200000000186, 702.96]
322542     [32.199999999953434, 2542.030000000028, 702.85]
372542       [37.40999999991618, 2542.930000000051, 703.2]
372547                  [37.75, 2547.359999999986, 703.51]
482542      [48.47999999998137, 2542.0200000000186, 704.2]
Length: 5546582, dtype: object

Function to find the elevation

This is the function used in the algorithm to find the elevation of any given point (x, y) using bilinear interpolation.

In [14]:
def find_elevation_new(x, y):
    # bilinear interpolation
    storage = []
    # curr_index = 10000*int(x)+int(y)
    diff = 0
    while(len(storage)<4):
        # Find all the data within the nearby grid
        temp = []
        # For loop here to find candidate points for bilinear interpolation
        for i in range(int(x-diff), int(x+diff+1)):
            for j in range(int(y-diff), int(y+diff+1)):
                if (i==int(x-diff) or i==int(x+diff) or j==int(y-diff) or j==int(y+diff)):
                    try:
                        rt = grid_dictionary[10000*i+j]
                        if (type(rt)==list):
                            temp.append(rt)
                        else:
                            for it in range(0,len(rt)):
                                temp.append(rt[it])
                    except Exception:
                        pass
        temp.sort(key = lambda e: (e[0]-x)*(e[0]-x)+(e[1]-y)*(e[1]-y))
        if (len(storage)+len(temp) <= 4):
            # add them
            for i in range(0,len(temp)):
                storage.append(temp[i])
        else:
            k = 0
            while(len(storage) != 4):
                storage.append(temp[k])
                k = k+1
        diff = diff+1
    # find the points used for data interpolation
    new_storage = storage[:4]
    a = [[1,new_storage[0][0], new_storage[0][1], new_storage[0][0]*new_storage[0][1]],
         [1, new_storage[1][0], new_storage[1][1], new_storage[1][0]*new_storage[1][1]],
         [1, new_storage[2][0], new_storage[2][1], new_storage[2][0]*new_storage[2][1]],
         [1, new_storage[3][0], new_storage[3][1], new_storage[3][0]*new_storage[3][1]]]
    b = [new_storage[0][2], new_storage[1][2], new_storage[2][2], new_storage[3][2]]
    try:
        # Conduct bilinear interpolation
        coef_matrix = np.matmul(inv(a), b)
        rt = coef_matrix[0]+coef_matrix[1]*x+coef_matrix[2]*y+coef_matrix[3]*x*y
        return rt
    except Exception:
        return 10000

Simulate the water flow

This is the major step in the algorithm that is used for simulation of the water flow.

This cell may run up to ~15 mins.</style>

In [15]:
# The area_length variable here represents how large the area the users want to calculate.
area_length = 10

# Print the current time before the execution of the algorithm
print(datetime.datetime.now())
min_x = var1*100
min_y = var2*100
max_x = var1*100+area_length
max_y = var2*100+area_length
increment = 1
angle = var3
x_coord = min_x
y_coord = min_y
rt = []
while x_coord!=(max_x+1):
    while y_coord!=(max_y+1):
        # find the elevation of the current coordinate
        curr_elevation = find_elevation_new(x_coord, y_coord)
        curr_x = x_coord
        curr_y = y_coord
        # print("Starting Point:")
        print((curr_x,curr_y))
        curr_array = []
        while ((curr_x>=min_x and curr_x<=max_x) and (curr_y>=min_y and curr_y<=max_y)):
            # print((curr_x, curr_y, curr_elevation))
            curr_array.append((curr_x, curr_y))
            rt_x = curr_x
            rt_y = curr_y
            rt_elevation = curr_elevation
            angel_diff = 0
            # find the elevation of all candidate points
            while(angel_diff<360):
                new_x = curr_x+increment*math.sin(math.pi/180*angel_diff)
                new_y = curr_y+increment*math.cos(math.pi/180*angel_diff)
                new_elevation = find_elevation_new(new_x, new_y)
                # print((angel_diff,new_x,new_y,new_elevation,rt_x,rt_y,rt_elevation))
                if (new_elevation<rt_elevation):
                    rt_x = new_x
                    rt_y = new_y
                    rt_elevation = new_elevation
                angel_diff = angel_diff+angle
            # Check if the data is a pit or not
            if (rt_elevation<curr_elevation):
                curr_x = rt_x
                curr_y = rt_y
                curr_elevation = rt_elevation
            else:
                break
        rt.append(curr_array)
        # print("Result For:")
        # print((x_coord,y_coord))
        # print(curr_array)
        y_coord=y_coord+1
    y_coord = min_y
    x_coord=x_coord+1
# Print the end time of the execution
print(datetime.datetime.now())
2021-09-16 22:16:01.143523
(500, 500)
(500, 501)
(500, 502)
(500, 503)
(500, 504)
(500, 505)
(500, 506)
(500, 507)
(500, 508)
(500, 509)
(500, 510)
(501, 500)
(501, 501)
(501, 502)
(501, 503)
(501, 504)
(501, 505)
(501, 506)
(501, 507)
(501, 508)
(501, 509)
(501, 510)
(502, 500)
(502, 501)
(502, 502)
(502, 503)
(502, 504)
(502, 505)
(502, 506)
(502, 507)
(502, 508)
(502, 509)
(502, 510)
(503, 500)
(503, 501)
(503, 502)
(503, 503)
(503, 504)
(503, 505)
(503, 506)
(503, 507)
(503, 508)
(503, 509)
(503, 510)
(504, 500)
(504, 501)
(504, 502)
(504, 503)
(504, 504)
(504, 505)
(504, 506)
(504, 507)
(504, 508)
(504, 509)
(504, 510)
(505, 500)
(505, 501)
(505, 502)
(505, 503)
(505, 504)
(505, 505)
(505, 506)
(505, 507)
(505, 508)
(505, 509)
(505, 510)
(506, 500)
(506, 501)
(506, 502)
(506, 503)
(506, 504)
(506, 505)
(506, 506)
(506, 507)
(506, 508)
(506, 509)
(506, 510)
(507, 500)
(507, 501)
(507, 502)
(507, 503)
(507, 504)
(507, 505)
(507, 506)
(507, 507)
(507, 508)
(507, 509)
(507, 510)
(508, 500)
(508, 501)
(508, 502)
(508, 503)
(508, 504)
(508, 505)
(508, 506)
(508, 507)
(508, 508)
(508, 509)
(508, 510)
(509, 500)
(509, 501)
(509, 502)
(509, 503)
(509, 504)
(509, 505)
(509, 506)
(509, 507)
(509, 508)
(509, 509)
(509, 510)
(510, 500)
(510, 501)
(510, 502)
(510, 503)
(510, 504)
(510, 505)
(510, 506)
(510, 507)
(510, 508)
(510, 509)
(510, 510)
2021-09-16 22:30:15.184070

Print the simulated water flow

In [16]:
print(rt)
    
data = []
# Print the simulated waterflow
for i in range(0,area_length+1):
    data.append([0]*(area_length+1))
[[(500, 500)], [(500, 501)], [(500, 502)], [(500, 503)], [(500, 504)], [(500, 505)], [(500, 506)], [(500, 507)], [(500, 508)], [(500, 509)], [(500, 510)], [(501, 500), (500.0, 500.0)], [(501, 501), (500.5, 501.8660254037844)], [(501, 502), (500.1339745962156, 501.5)], [(501, 503), (500.1339745962156, 502.5)], [(501, 504), (500.0, 504.0)], [(501, 505), (500.5, 504.1339745962156)], [(501, 506), (500.0, 506.0)], [(501, 507), (500.1339745962156, 507.5)], [(501, 508), (500.1339745962156, 508.5)], [(501, 509), (500.1339745962156, 509.5)], [(501, 510)], [(502, 500), (501.0, 500.0), (500.0, 500.0)], [(502, 501), (501.0, 501.0), (500.5, 501.8660254037844)], [(502, 502), (502.5, 501.1339745962156)], [(502, 503), (501.0, 503.0), (500.1339745962156, 502.5)], [(502, 504), (501.1339745962156, 503.5), (500.1339745962156, 503.5)], [(502, 505), (501.5, 504.1339745962156), (500.5, 504.1339745962156)], [(502, 506), (502.5, 506.8660254037844)], [(502, 507), (501.1339745962156, 507.5), (500.26794919243116, 508.0)], [(502, 508), (501.5, 508.8660254037844), (500.6339745962156, 509.3660254037844)], [(502, 509), (501.1339745962156, 509.5), (500.26794919243116, 510.0)], [(502, 510)], [(503, 500), (502.0, 500.0), (501.0, 500.0), (500.0, 500.0)], [(503, 501), (502.0, 501.0), (501.0, 501.0), (500.5, 501.8660254037844)], [(503, 502), (502.5, 501.1339745962156)], [(503, 503), (502.0, 503.0), (501.0, 503.0), (500.1339745962156, 502.5)], [(503, 504), (502.0, 504.0), (501.1339745962156, 503.5), (500.1339745962156, 503.5)], [(503, 505), (502.5, 504.1339745962156), (501.5, 504.1339745962156), (500.5, 504.1339745962156)], [(503, 506), (502.5, 506.8660254037844)], [(503, 507), (502.5, 507.8660254037844), (502.5, 506.8660254037844)], [(503, 508), (502.5, 508.8660254037844), (501.5, 508.8660254037844), (500.6339745962156, 509.3660254037844)], [(503, 509), (502.1339745962156, 509.5), (501.26794919243116, 510.0)], [(503, 510)], [(504, 500)], [(504, 501), (503.0, 501.0), (502.0, 501.0), (501.0, 501.0), (500.5, 501.8660254037844)], [(504, 502), (503.1339745962156, 502.5), (502.26794919243116, 503.0), (501.40192378864674, 502.5), (500.5358983848623, 502.0)], [(504, 503), (503.0, 503.0), (502.0, 503.0), (501.0, 503.0), (500.1339745962156, 502.5)], [(504, 504), (503.0, 504.0), (502.0, 504.0), (501.1339745962156, 503.5), (500.1339745962156, 503.5)], [(504, 505), (503.1339745962156, 504.5), (502.26794919243116, 505.0), (502.76794919243116, 504.1339745962156)], [(504, 506), (503.0, 506.0), (502.5, 506.8660254037844)], [(504, 507), (503.1339745962156, 507.5), (502.6339745962156, 508.3660254037844), (502.1339745962156, 509.23205080756884), (501.26794919243116, 509.73205080756884)], [(504, 508), (503.1339745962156, 508.5), (502.26794919243116, 509.0), (501.40192378864674, 509.5)], [(504, 509), (503.1339745962156, 509.5), (502.26794919243116, 510.0)], [(504, 510)], [(505, 500), (506.0, 500.0)], [(505, 501), (504.1339745962156, 501.5), (503.26794919243116, 502.0), (502.40192378864674, 502.5), (501.40192378864674, 502.5), (500.5358983848623, 502.0)], [(505, 502), (504.0, 502.0), (503.1339745962156, 502.5), (502.26794919243116, 503.0), (501.40192378864674, 502.5), (500.5358983848623, 502.0)], [(505, 503), (504.1339745962156, 503.5), (503.26794919243116, 504.0), (502.26794919243116, 504.0), (501.26794919243116, 504.0), (500.26794919243116, 504.0)], [(505, 504), (504.0, 504.0), (503.0, 504.0), (502.0, 504.0), (501.1339745962156, 503.5), (500.1339745962156, 503.5)], [(505, 505), (504.0, 505.0), (503.1339745962156, 504.5), (502.26794919243116, 505.0), (502.76794919243116, 504.1339745962156)], [(505, 506), (504.1339745962156, 506.5), (503.6339745962156, 507.3660254037844), (502.76794919243116, 507.8660254037844), (502.26794919243116, 508.73205080756884), (501.76794919243116, 509.59807621135326)], [(505, 507), (504.1339745962156, 506.5), (503.6339745962156, 507.3660254037844), (502.76794919243116, 507.8660254037844), (502.26794919243116, 508.73205080756884), (501.76794919243116, 509.59807621135326)], [(505, 508), (504.1339745962156, 507.5), (503.26794919243116, 508.0), (502.40192378864674, 508.5), (501.90192378864674, 509.3660254037844)], [(505, 509), (504.1339745962156, 509.5), (503.26794919243116, 510.0)], [(505, 510)], [(506, 500)], [(506, 501), (506.0, 500.0)], [(506, 502), (506.8660254037844, 502.5)], [(506, 503), (506.8660254037844, 502.5)], [(506, 504), (505.0, 504.0), (504.0, 504.0), (503.0, 504.0), (502.0, 504.0), (501.1339745962156, 503.5), (500.1339745962156, 503.5)], [(506, 505), (505.0, 505.0), (504.0, 505.0), (503.1339745962156, 504.5), (502.26794919243116, 505.0), (502.76794919243116, 504.1339745962156)], [(506, 506), (505.1339745962156, 505.5), (504.1339745962156, 505.5), (504.1339745962156, 506.5), (503.6339745962156, 507.3660254037844), (502.76794919243116, 507.8660254037844), (502.26794919243116, 508.73205080756884), (501.76794919243116, 509.59807621135326)], [(506, 507), (505.1339745962156, 507.5), (504.1339745962156, 507.5), (503.26794919243116, 508.0), (502.40192378864674, 508.5), (501.90192378864674, 509.3660254037844)], [(506, 508), (505.1339745962156, 507.5), (504.1339745962156, 507.5), (503.26794919243116, 508.0), (502.40192378864674, 508.5), (501.90192378864674, 509.3660254037844)], [(506, 509), (505.5, 509.8660254037844)], [(506, 510)], [(507, 500), (506.0, 500.0)], [(507, 501), (506.1339745962156, 500.5)], [(507, 502), (506.1339745962156, 502.5), (505.6339745962156, 501.6339745962156), (504.76794919243116, 502.1339745962156), (503.76794919243116, 502.1339745962156), (502.90192378864674, 502.6339745962156), (502.0358983848623, 503.1339745962156), (501.1698729810779, 503.6339745962156), (500.3038475772935, 504.1339745962156)], [(507, 503), (506.0, 503.0), (506.8660254037844, 502.5)], [(507, 504), (506.0, 504.0), (505.0, 504.0), (504.0, 504.0), (503.0, 504.0), (502.0, 504.0), (501.1339745962156, 503.5), (500.1339745962156, 503.5)], [(507, 505), (507.8660254037844, 505.5)], [(507, 506), (507.8660254037844, 505.5)], [(507, 507), (507.5, 507.8660254037844)], [(507, 508), (507.5, 508.8660254037844), (507.5, 507.8660254037844)], [(507, 509), (506.1339745962156, 509.5)], [(507, 510), (506.0, 510.0)], [(508, 500)], [(508, 501), (507.5, 500.1339745962156)], [(508, 502), (507.0, 502.0), (506.1339745962156, 502.5), (505.6339745962156, 501.6339745962156), (504.76794919243116, 502.1339745962156), (503.76794919243116, 502.1339745962156), (502.90192378864674, 502.6339745962156), (502.0358983848623, 503.1339745962156), (501.1698729810779, 503.6339745962156), (500.3038475772935, 504.1339745962156)], [(508, 503), (507.0, 503.0), (506.0, 503.0), (506.8660254037844, 502.5)], [(508, 504), (507.1339745962156, 503.5), (506.26794919243116, 503.0), (505.26794919243116, 503.0), (504.40192378864674, 503.5), (503.5358983848623, 504.0), (502.5358983848623, 504.0), (501.5358983848623, 504.0), (500.6698729810779, 503.5), (500.1698729810779, 504.3660254037844)], [(508, 505), (507.5, 505.8660254037844), (506.6339745962156, 506.3660254037844), (506.1339745962156, 507.23205080756884), (505.1339745962156, 507.23205080756884)], [(508, 506), (507.1339745962156, 506.5), (507.6339745962156, 505.6339745962156)], [(508, 507), (507.5, 507.8660254037844)], [(508, 508), (507.5, 508.8660254037844), (507.5, 507.8660254037844)], [(508, 509), (509.0, 509.0)], [(508, 510), (507.5, 509.1339745962156), (507.0, 508.26794919243116)], [(509, 500), (508.0, 500.0)], [(509, 501), (508.1339745962156, 500.5)], [(509, 502), (508.0, 502.0), (507.0, 502.0), (506.1339745962156, 502.5), (505.6339745962156, 501.6339745962156), (504.76794919243116, 502.1339745962156), (503.76794919243116, 502.1339745962156), (502.90192378864674, 502.6339745962156), (502.0358983848623, 503.1339745962156), (501.1698729810779, 503.6339745962156), (500.3038475772935, 504.1339745962156)], [(509, 503), (508.0, 503.0), (507.0, 503.0), (506.0, 503.0), (506.8660254037844, 502.5)], [(509, 504), (508.1339745962156, 504.5), (507.26794919243116, 504.0), (506.26794919243116, 504.0), (505.26794919243116, 504.0), (504.26794919243116, 504.0), (503.26794919243116, 504.0), (502.26794919243116, 504.0), (501.26794919243116, 504.0), (500.26794919243116, 504.0)], [(509, 505), (508.1339745962156, 505.5), (507.26794919243116, 506.0), (507.76794919243116, 505.1339745962156)], [(509, 506), (508.1339745962156, 506.5), (507.6339745962156, 505.6339745962156)], [(509, 507), (508.1339745962156, 507.5), (507.1339745962156, 507.5), (507.1339745962156, 508.5)], [(509, 508), (509.0, 509.0)], [(509, 509)], [(509, 510), (509.0, 509.0)], [(510, 500)], [(510, 501), (509.0, 501.0), (508.1339745962156, 500.5)], [(510, 502), (509.0, 502.0), (508.0, 502.0), (507.0, 502.0), (506.1339745962156, 502.5), (505.6339745962156, 501.6339745962156), (504.76794919243116, 502.1339745962156), (503.76794919243116, 502.1339745962156), (502.90192378864674, 502.6339745962156), (502.0358983848623, 503.1339745962156), (501.1698729810779, 503.6339745962156), (500.3038475772935, 504.1339745962156)], [(510, 503), (509.1339745962156, 502.5), (508.1339745962156, 502.5), (507.1339745962156, 502.5), (506.26794919243116, 503.0), (505.26794919243116, 503.0), (504.40192378864674, 503.5), (503.5358983848623, 504.0), (502.5358983848623, 504.0), (501.5358983848623, 504.0), (500.6698729810779, 503.5), (500.1698729810779, 504.3660254037844)], [(510, 504), (509.1339745962156, 504.5), (508.6339745962156, 505.3660254037844), (507.76794919243116, 505.8660254037844)], [(510, 505), (509.1339745962156, 505.5), (508.1339745962156, 505.5), (507.26794919243116, 506.0), (507.76794919243116, 505.1339745962156)], [(510, 506), (510.0, 507.0), (510.0, 508.0)], [(510, 507), (510.0, 508.0)], [(510, 508)], [(510, 509), (509.0, 509.0)], [(510, 510), (509.1339745962156, 509.5)]]
2:1: W293 blank line contains whitespace

Track the flow direction

In [17]:
for i in range(0,len(rt)):
    for j in range(0,len(rt[i])):
        data[int(math.floor(rt[i][j][0]))-var1*100][int(math.floor(rt[i][j][1]))-var2*100]=data[int(math.floor(rt[i][j][0]))-var1*100][int(math.floor(rt[i][j][1]))-var2*100]+1

Visulization

The result after tracking the water direction

In [18]:
print(data)
[[4, 6, 8, 9, 13, 1, 2, 2, 3, 4, 2], [3, 4, 4, 14, 7, 1, 1, 2, 3, 10, 2], [2, 5, 6, 9, 14, 4, 5, 5, 9, 4, 2], [1, 2, 8, 2, 12, 1, 2, 5, 5, 2, 2], [1, 2, 6, 4, 5, 4, 4, 4, 1, 2, 1], [1, 5, 1, 3, 4, 3, 1, 4, 1, 2, 1], [5, 1, 10, 6, 3, 1, 2, 2, 1, 2, 2], [2, 1, 5, 4, 2, 9, 4, 6, 5, 2, 1], [4, 1, 4, 2, 2, 4, 2, 2, 1, 1, 1], [1, 2, 3, 1, 2, 2, 1, 1, 1, 6, 1], [1, 1, 1, 1, 1, 1, 1, 2, 3, 1, 1]]

Simple visulization of the result, with brighter area having more water flowing through

In [19]:
%matplotlib inline
import matplotlib.pyplot as plt

plt.imshow(data, interpolation='none')
plt.show()