How does this package implement grid runoff calculation in a distributed hydrological model?

[chooron]

I will describe this issue in more detail. Firstly, we have a geographic elevation dataset, which we will convert into flow direction (D8 flow direction). This means that a standard grid cell has eight neighbors. Based on each cell's unique D8, we can determine that the current state of the cell will transfer to one of its eight neighbors. I checked the documentation, and Neighbors seems to be one of the keys to solving this problem. However, I am unsure how to assign a state to a specific neighbor based on the flow direction. I suspect that rules should be constructed for the eight flow directions, and when a flow direction matches a rule, that rule is executed. Below is a simple implementation of my requirement without using this package, but this computation does not seem to be differentiable:

using PaddedViews

function d8_pads(d8_shape)
    x_shape, y_shape = d8_shape[1] + 2, d8_shape[2] + 2
    d8_pads_dims = [
        (1:x_shape-2, 1:y_shape-2),
        (1:x_shape-2, 2:y_shape-1),
        (1:x_shape-2, 3:y_shape),
        (2:x_shape-1, 1:y_shape-2),
        (2:x_shape-1, 3:y_shape),
        (3:x_shape, 1:y_shape-2),
        (3:x_shape, 2:y_shape-1),
        (3:x_shape, 3:y_shape),
    ]
    d8_pads_dims
end

function grid_routing(flow::AbstractMatrix, flwdir::AbstractMatrix)
    d8_codes = [0, 2, 4, 8, 16, 32, 64, 128]
    pad_dims = d8_pads(size(flow))
    pad_shape = (1:size(flow)[1]+2, 1:size(flow)[2]+2)
    flow_routed = sum([collect(PaddedView(eltype(flow)(0.0), flow .* (flwdir .== code), pad_shape, pad_dim))
                       for (code, pad_dim) in zip(d8_codes, pad_dims)])
    flow_routed[2:size(flow)[1]+1, 2:size(flow)[2]+1]
end

In addition, at each time point, each cell will generate a new state (rainfall generating flow) and combine the input values from other cells in the previous period to construct a new state. I am also unsure how to implement this functionality.

Finally, I have another question about this package: how can we set the computation order of each cell (in some cases, it might need to be determined based on the cumulative flow value of each cell generated by the flow direction)?

[rafaqz]

The fast way to do this is to assign the value of the current cell based on the neighbours, like a convolution. For that you can just use a simple Neighbours rule and return a new cell value based on its neighbours. So water flows inwards or outwards but you only assign the change to the current cell.

If you need to do it the other way around and set neighbour values for some reason, you can use SetNeighbors and update the neighbouring cell values with add!. You would iterate over them with for n in indices(neighbourhood, I). And finally update the centre cell.

[rafaqz]

Finally, I have another question about this package: how can we set the computation order of each cell (in some cases, it might need to be determined based on the cumulative flow value of each cell generated by the flow direction)?

This is not something you can (or should) do in DynamicGrids.jl. The concept is the same as cellular automata: everything happens at the same time! then updates are done for the next time-step.

If you need an iterative model, just write it using Stencils.jl directly (it does all the windowing here) then you can just extract stencils wherever you like in a loop.

[chooron]

Thank you very much for your suggestion. I will try to improve my code based on your recommendations. Once completed, I can upload the implemented code, hoping it will be useful for others.