Flow Direction Raster format

[FarisSquared]

The documentation in the input files mentions that we can generate flow direction from ArcMap and QGIS.

I am using QGIS and GRASS's r.fill.dir, the issue I am facing is that the coding of (agnps) mode is different to the one required in GIStoSWMM.

• 1 in agnps is north
• 1 required in GIStoSWMM is northeast

so,
Could you please share with us how you generated the raster?

as a temporary solution, I am writing a quick script to fix my coding. but curious about how to do it because it will simplify the workflow.

Thanks,
Faris

[FarisSquared]

it also seems like ArcGIS pro uses a different encoding for direction.

[FarisSquared]

Leaving this here in case someone needs it,

import rasterio 
import numpy as np 

with rasterio.open("flow_direction_raster.tiff") as flow:
    flow_arr = flow.read(1).copy() # one here means band number 1 (the only band)
    #now we define the encoding needed by GIStoSWMMM
    GIStoSWMM_mapping = {
        1:2,
        2:3,
        3:4,
        4:5,
        5:6,
        6:7,
        7:8,
        8:1,
        255:255 #nodata
    }
    # now iterate over matrix and replace values with GIStoSWMM encoding  
    for cell in np.nditer(flow_arr, op_flags=['readwrite']):
        cell[...] = GIStoSWMM_mapping[cell.item()] 

    # saving the changed array into a raster 
    with rasterio.open("flow_direction_raster_GIStoSWMM.tif", 'w+', **flow.profile) as out:
        out.write_band(1,flow_arr)

[schoeller]

@fariscs
@tjniemi
Good day. I have reformatted the flowdirection raster originating from pyflwdir module as below based on definitions as per

https://deltares.github.io/pyflwdir/latest/flwdir.html?highlight=d8
https://github.com/AaltoUrbanWater/GisToSWMM5/blob/master/input_files.md

with rasterio.open("raster_dem.tif", "r") as src:
    elevtn = src.read(1)
    nodata = src.nodata
    transform = src.transform
    crs = src.crs
    extent = np.array(src.bounds)[[0, 2, 1, 3]]
    prof = src.profile
    
    flw = pyflwdir.from_dem(
        data=elevtn,
        nodata=src.nodata,
        transform=transform,
        outlets="min",
    )

d8_data = flw.to_array(ftype="d8")
prof.update(dtype=d8_data.dtype, nodata=255)
with rasterio.open("raster_flowdir.tif", "w", **prof) as src:
    src.write(d8_data, 1)

# Mapping flow direction values to GIStoSWMM    
with rasterio.open("raster_flowdir.tif") as flow:
    flow_arr = flow.read(1).copy() # one here means band number 1 (the only band)
    #now we define the encoding needed by GIStoSWMMM
    GIStoSWMM_mapping = {
        16:5,
        8:6,
        4:7,
        32:4,
        0:9,
        2:8,
        64:3,
        128:2,
        1:1,
        247:255 #nodata
    }
    # now iterate over matrix and replace values with GIStoSWMM encoding  
    for cell in np.nditer(flow_arr, op_flags=['readwrite']):
        cell[...] = GIStoSWMM_mapping[cell.item()] 

    # saving the changed array into a raster 
    with rasterio.open("raster_flowdir_gistoswmm.tif", 'w+', **flow.profile) as out:
        out.write_band(1,flow_arr)

I have tried to map D8/pyflwdir

32	64	128		↖️	⬆️	↗️
16	X	1		⬅️	X	➡️
8	4	2		↙️	⬇️	↘️

to GisToSwmm D8 definition

3	2	1		↖️	⬆️	↗️
4	X	8		⬅️	X	➡️
5	6	7		↙️	⬇️	↘️

The flowdir raster looks correct, but still throws errors like

-> Error, no data to calculate slope between cells 168214 and 167743.

The file header looks like

ncols        470
nrows        358
xllcorner    3476805.690999999642
yllcorner    5532655.959999999031
cellsize     1.000000000000
NODATA_value 255
 255 255 255 255 255 255 255

If you have any idea how this behaviour might arrive, please kindly share: hints are very welcome.

[FarisSquared]

Hi @schoeller,

First, I could not run this project using my data as I kept hitting a "segmentation error" that I could not figure out.

Second, it is okay to get the no slope error if you have no_data values in your raster.

Third, double-check the mapping. I think you should have 128:1 rather than 128:2 and so on.

Finally, Best of luck!
Please let me know if you succeed with running this (If you do, I will give it another try with my dataset)

I hope @tjniemi has time to review all this, as I like the idea of this project.

[schoeller]

@fariscs
I altered

    GIStoSWMM_mapping = {
        16:5,
        8:6,
        4:7,
        32:4,
        0:9,
        2:8,
        64:3,
        128:2,
        1:1,
        247:255 #nodata

to

        32:3,
        16:4,
        8:5,
        4:6,
        2:7,
        1:8,
        128:1,
        64:2,
        0:9,
        247:255 #nodata

Then I had to set NODATA_valueto 0 within landuse. I did likely to dem and flowdir. I then ran a modified batch and received

Report:
-> Running time: 0.0876675 minutes
-> Total number of cells: 146
-> Number of active cells: 146
-> Catchment area: 5.6865 ha
-> Catchment average elevation: 0 m
-> Catchment average slope: 0
-> Landuse information (code, number of cells, area in ha, % of catchment):
        10      54      2.3597  41.4965
        12      51      1.9333  33.9981
        14      41      1.3935  24.5054
        20      0       0       0
        30      0       0       0
        32      0       0       0
        34      0       0       0
        60      0       0       0
        62      0       0       0
        64      0       0       0

I expect non-zero average catchment elevation and suspect an error in my input ASC files. If @tjniemi could throw an eye I be more than glad.

[FarisSquared]

Excellent work!

I expect non-zero average catchment elevation and suspect an error in my input ASC files. If @tjniemi could throw an eye I be more than glad.

As for the zero here, it makes sense. since in the line linked. it checks before committing the addition for two things

1. that the cell's landuse value is larger than BUILT_AREA. judging from the report you posted it is not the case for all cells you have. unless the report is also incorrect?
2. that it is larger than no_data which assigned from the DEM,
   relevant lines:

• reading the noData from the raster as string
• converting noValueData from Raster to double then assigning it to each cell

I did not catch that during my time with this project so good job.
I think the second check should be equality because zero or minus values are valid as an elevation value.

[siamak6566]

@schoeller

Have a good day
I am working with GisToSWMM5 and somehow new in. During my SWMM set up I got a report which is very similar to yours in terms of non-zero average catchment elevation and recognizing the landuse classes. I also guess that there is something wrong with my ASC files. I received the following report:

Report:
-> Running time: 21.4897 minutes
-> Total number of cells: 2494
-> Number of active cells: 2494
-> Catchment area: 32.9963 ha
-> Catchment average elevation: 0 m
-> Catchment average slope: 0
-> Landuse information (code, number of cells, area in ha, % of catchment):
31 0 0 0
34 0 0 0
30 0 0 0
10 2494 32.9963 100
33 0 0 0
32 0 0 0
60 0 0 0
61 0 0 0

As can be seen, only building with code 10 is recognized by the tool as sub-catchments. Additionally, catchment average elevation and slope are 0 which made me suspicious about the ASC files.
Would you let me know how you did overcome with these issues in your project?

Thank you in advance