Subcatchment routing hierarchy by landuse

[schoeller]

Dear all,

I am trying to resemble a functionality similar to GisToSWMM using pyflwdir. In my understanding of this algorithm a subcatchment hierarchy is determined considering landuse. Every (sub-)subcatchment contains only one landuse and routes either to another (sub-)subcatchment or an outlet. Outlets are given in the form of enlarged nodes, thus rasterized sink information derived from outlet nodes and DEM sinks.

Reading through pyflwdir documentation I figure that the Pfafstetter method allows for deriving hierarchical subcatchment information. If understood correctly it permits to figure out a form of subcatchment routing, thus which subcatchment drains to another subcatchment or outlet.

I am seeking ways to mask this routing method using landuse, so that every subcatchment contains only one landuse. Any ideas if and/or how to achieve this goal are highly appreciated.

Kind regards

Sebastian

[DirkEilander]

Hi Sebastian,

Nice application! I think that should be possible but not with the pfafstetter method. I would suggest to try the following:

1. Create regions from the landuse map, using e.g. scipy.ndimage.label. This returns a map with unique IDs for connected regions with the same landuse.
2. Find the most downstream "outflow pixel" of each region using pyflwdir.FlwdirRaster.inflow_idxs Note: should still be added to the API Reference in the docs. This returns a list with linear indices of outflow pixels.
3. Derive subbasins using pyflwdir.FlwdirRaster.basins based on the indices of the outflow pixels from step 2, supplied using the idxs_out option. This returns a map with unique IDs for each basin.
4. Find the relations between the basins using pyflwdir.FlwdirRaster.streams, again based on the indices of the outflow pixels from step 2, supplied using the idxs_out option. This returns geofeatures with river segments between outflow pixels and per segment the index of the next downstream segment (and thus basin).

Hope this is helpful!

[schoeller]

@DirkEilander many thanks for the detailed response.

I am stumbling in step 1. I have read a geodataframe for landuse from a gpkg layer and converted to xarray + calling of label as follows:

lu_grid = make_geocube(
    vector_data=gdf_lu,
    output_crs=gdf_lu.crs.to_epsg(),
    resolution=(-1, 1),
)
s = generate_binary_structure(2,2)
lu_labeled_array, lu_num_features = label(lu_grid.type, s)

lu_num_features contains one solutions only which is not desired taking below variable content:

print(lu_grid.type)
<xarray.DataArray 'type' (y: 231, x: 204)>
array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       ...,
       [nan, nan, nan, ..., 30., 61., nan],
       [nan, nan, nan, ..., 61., 30., nan],
       [nan, nan, nan, ..., 61., 61., 30.]])
Coordinates:
  * y            (y) float64 230.5 229.5 228.5 227.5 226.5 ... 3.5 2.5 1.5 0.5
  * x            (x) float64 17.5 18.5 19.5 20.5 ... 217.5 218.5 219.5 220.5
    spatial_ref  int32 0
Attributes:
    name:        type
    long_name:   type
    _FillValue:  nan

Hints on where to dig further are highly welcome.

[DirkEilander]

Note that this goes beyond pyflwdir, so I'm just guessing here. You might have to pass a numpy.ndarray instead of xarray.DataArray to the label method, so lu_labeled_array, lu_num_features = label(lu_grid.type.values, s). I'm also not sure how the label method deals with nan values.

[schoeller]

@DirkEilander looks as if this works to deliver a dataarray with reasonable results. Crawling along...

from xrspatial.zonal import regions
lu_regions = regions(raster=lu_grid.type,neighborhood=8)

[DirkEilander]

@schoeller Curious to hear how this is going.

[schoeller]

@DirkEilander Efforts have been laying low since 1 month now. Have documented scripting efforts under https://gist.github.com/schoeller/ed0cfa4c163edaa4cde19d8cf8c891ef

[schoeller]

@DirkEilander: I have started over my efforts. Being nooby I have stumbled over your first point. I have used the original raster files from the reference project. The file looks as below:

I have tried the following code:

import rasterio
import scipy
with rasterio.open('GIS/raster_landuse.tif', nodata=0) as src:
    labeled_array, num_features = scipy.ndimage.label(src.read(1).astype(int))
    print(num_features)

num_features throws 1 as a result, which I do not consider reasonable. Maybe you can spare time to pinpoint me into the right direction.

Best

Sebastian

[DirkEilander]

Hi Sebastian,

From the docs of scipy.ndimage.label method: "Any non-zero values in input are counted as features and zero values are considered the background." Seeing the picture if looks like most non-zero (greyscale) pixels are connected right? If you use want just want to label the black areas you could try adjusting the input the label method using e.g. "src.read(1).astype(int) > x".
Note that this is not pyflwdir functionality so for help on this method it's best to use Stack Overflow or the scipy issue trackers.