Subcatchment routing hierarchy by landuse about pyflwdir HOT 8 OPEN

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!

from pyflwdir.

@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.

from pyflwdir.

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.

from pyflwdir.

@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)

from pyflwdir.

@schoeller Curious to hear how this is going.

from pyflwdir.

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

from pyflwdir.

@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

from pyflwdir.

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.

from pyflwdir.