Open schoeller opened 2 years ago
Hi Sebastian,
Nice application! I think that should be possible but not with the pfafstetter method. I would suggest to try the following:
idxs_out
option. This returns a map with unique IDs for each basin.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!
@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.
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.
@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)
@schoeller Curious to hear how this is going.
@DirkEilander Efforts have been laying low since 1 month now. Have documented scripting efforts under https://gist.github.com/schoeller/ed0cfa4c163edaa4cde19d8cf8c891ef
@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
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.
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