[GRASS-user] Classify basins as "narrow"

[Stefan Blumentrath]

Hei Ken,

What about a combination of r.grow.distance and r.neighbors?

1) Extract boundaries from raster basins
r.neighbors input=basins output=basin_diversity method=diversity
r.reclass input=basin_diversity output=basin_boundary rc=-
"1 = NULL
2 thru 9999999 = 1" (if you need the boundaries later on the vector solution is probably more efficient)

2) compute distance from boundaries
r.grow.distance input=basin_boundary output=basin_boundary_distance

3) get the maximum distance from the boundaries
r.neighbors method=maximum size="twice what you would consider narrow"

Depends if it is computationally more efficient... If narrow still means "a lot" of pixels in width, then r.neighbors might be relatively slow. But you might increase resolution in that case.

In order to get the skeleton of the basin you could use r.neighbors as well:
r.neighbors input=basin_boundary_distance method=maximum output=basin_boundary_distance_max
r.mapcalc expression="skeletons=if(basin_boundary_distance < basin_boundary_distance_max, null(), 1)"
This should approximate the center lines...
Note that you could skip r.neighbors and use the neighborhood modifier in r.mapcalc...

Just some ideas for the width issue...

Cheers,
Stefan

-----Original Message-----
From: grass-user <grass-user-bounces at lists.osgeo.org> On Behalf Of Ken Mankoff
Sent: onsdag 2. september 2020 20:03
To: GRASS user list <grass-user at lists.osgeo.org>
Subject: Re: [GRASS-user] Classify basins as "narrow"

Hi List,

On 2020-09-02 at 04:27 -07, Ken Mankoff <mankoff at gmail.com> wrote...
> I'd like to detect "narrow" features in GRASS. The attached screenshot 
> shows some basins (thick) and streams (thin) and some regions 
> (hatched). These regions are spurious because the basin is narrow 
> here. I'd like to estimate narrowness with an algorithm.
>
> I've looked into r.grow.distance r.distance and v.distance but haven't 
> been able to imagine a solution yet. Can anyone on this list suggest 
> something?

I've solved this, although it takes a lot of steps, and it only computes the width of the downstream region, not everywhere. The downstream/outflow region is what I'm interested in (but didn't specify this in my initial post) so this is OK with me. By this I mean that for a "Y" shaped catchment with three narrow regions and outflow at the bottom of the stem, the current algorithm will compute the width of one of the branches and the stem, but not the other branch.

I'm sharing this in case it helps someone else, or someone sees a way to do this more efficiently (perhaps all basins at once?). Here is the algorithm:

# elevation raster is input

# compute streams, outlets, and basins

r.stream.extract elevation=elevation threshold=11 memory=16384 direction=dir stream_raster=streams stream_vector=streams

r.mapcalc "outlets = if(dir < 0, 1, null())" # outlets

r.to.vect input=outlets output=outlets type=point

r.stream.basins -m direction=dir points=outlets basins=basins memory=16384 # basins

# compute distance from edge of each basin

r.to.vect -v input=basins output=basins type=area

v.to.lines input=basins output=bounds

v.to.rast -d input=bounds output=bounds use=val val=1

r.grow.distance input=bounds distance=edge_dist metric=euclidean

# Now the algorithm can only work on one basin at a time.
# Pick a basin with narrow outlet region manually, and mask to that basin.

basin_id=1174

r.mask raster=basins maskcats=${basin_id} --o

# Invert so center line is low, but make outlet lowest.
# Edge dist is 0 at edges, so this cost surface is a hole # Make the outlet the lowest point so it pours from there.
r.mapcalc "cost_surf = if(isnull(outlets), -edge_dist, -max(edge_dist)-100)"

# Route a stream down the center line

r.stream.extract elevation=cost_surf threshold=11 memory=16384 direction=cost_dir stream_raster=cost_streams stream_vector=cost_streams

# Find the main stream (hack = 1)

r.stream.order -m stream_rast=cost_streams direction=cost_dir elevation=cost_surf hack=hack memory=16384

r.mapcalc "hack_1 = if(hack == 1, 1, null())"

# Grow center stream to edges

r.grow.distance -m input=hack_1 distance=center_dist

# DONE

Now the value of center_dist at the basin boundary is the distance to the center line, or 1/2 the width.

_______________________________________________
grass-user mailing list
grass-user at lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user