GRASS GIS: Computing center line of raster feature representing a river

Question

My goal:

Convert a raster layer representing a river feature (see picture) to the center line of the river stored as a polyline feature.

Limitation:

• The whole process has to be automated.
• I can only use GRASS GIS for this application.

The initial input is a raster resulting from a scanner. From this scan I have reclassified the river.

I have tried the following so far:

• Convert to polygon --> convert to line (out-lines e.g. river banks) --> convert river banks to point --> apply voronoi diagrams and clip this output to the extent of your river --> see output. I tried using v.clean to clean the output e.g. remove dangles but this has not worked yet.
• r.grow to thin my raster or r.thin to thin it --> convert thinned raster to a line --> v.clean to clean the result. This result (black lines) contains "loops" which I did not manage to delete.

Especially here "v.generalize" with a high threshold yielded some better results but still not all loops are deleted. So, you can consider this tool as well in a possible solution!

Any recommendations on a different approach or to alter my previous attempts? In an other post here (Finding the Center Line from a set of 3D Points) they mention a possible interesting approach suggested by Joseph.

EDIT: 28.05.2019 Will work on this tomorrow again and test again my approach which is also mentioned by Micha. Jons' solution seems to be the most robust approach however with my skill and time budget I am not sure if I can reproduce a algorithm which will work in GRASS GIS.

Answer 1

I have solved this problem numerous times. I have a Matlab solution posted here, and a Python solution is included in my (soon to be released) RivGraph package here. The function you'd be interested in is called "mask_to_centerline." It takes the approach of skeletonizing (or thinning) and finding the longest, shortest path between end-locations as the centerline. There are a few other tricks in there as well.

You might want to look at polygonizing your raster, then using Voroni polygons to approximate a centerline. This approach seems to work well (I have never implemented it) and is the basis of Wes Lauer's ArcGIS Planform Statistics Toolbox's centerline function, although that requires banklines as inputs.

As you are married to GRASS, I cannot provide you some cut-and-pasteable code, but you can follow the logic of the scripts (the RivGraph Python script in particular) to maybe craft your own. One tricky challenge of this problem is to handle the edges, where the thinning tends to connect to the corner of the river object. I use mirroring to overcome this, but it doesn't always work.

Here is the output of RivGraph's mask_to_centerline.py on your image:

from rivgraph.rivers import river_utils as ru
from scipy import misc
import numpy as np
from matplotlib import pyplot as plt

# Load image
I = misc.imread(r"C:\Users\Jon\Desktop\test\yO67j.png")

# Make mask
Imask = np.zeros((I.shape[0], I.shape[1]), dtype=np.bool)
Imask[I[:,:,0] != 255] = True

# Get centerline
coords, widths = ru.mask_to_centerline(Imask, 'ew')

# Plot
plt.imshow(Imask)
plt.plot(coords[:,0], coords[:,1])

I have used the mask_to_centerline.py for dozens of rivers, and just yesterday on 80 different tidal meanders that were given as polygons in Supplementary Info of this paper to extract tidal meander centerlines. So I can say that it's fairly robust, and if you follow its logic you should end up with a similar result.

If you use Anaconda for package management, you may also try to install the RivGraph package from conda using conda install -c jschwenk rivgraph BUT I must warn you that I am still trying to debug that installation for non-Windows OSes (I also recommend installing in its own environment).

Answer 2

I think that the GRASS commands r.thin, and then r.to.vect ... type=line should do what you want. You have to prepare the scanned raster first by setting all the pixels that cover the river to a single value (it seems to be so from your image), and all the surrounding pixels must be NULL (r.null setnull=...)

(Edit:) Here's the procedure I tried and the result. I started with a "fat" river.

# Make an integer copy of the original (for r.thin)
r.mapcalc "str_tmp2 = int(str_tmp)" --o
# Run r.thin
r.thin str_tmp2 out=str_tmp_thin --o
# Convert to vector
r.to.vect str_tmp_thin out=str_tmp_thin type=line --o

The red line is the final vector. Dark colors are the original, and the light green is the thinned raster.

HTH