GRW algorithm ver1.5

3 April, 2014

1

The Global River Width Algorithm

2
3
4

Dai Yamazaki

5

School of Geographical Sciences, University of Bristol
Dai.Yamazaki@bristol.ac.uk

Note:

6

This document describes the Global River Width Algorithm (GRW Algorithm), which was7

used to develop the Global Width Database for Large Rivers (GWD-LR) [Yamazaki et al., in82014]. The Fortran90 code of the GRWA with sample dataset is available on request to the9developer (Dai Yamazaki).

10

1. Input and Output Datasets

11

1.1 Input Datasets

12

The GRW algorithm requires three input datasets: (1) a water body mask; (2) a flow

13

direction map; and (3) a drainage area map. These three input datasets should be prepared14

in a consistent grid coordination system and at the same spatial resolution. The SRTM

15

Water Body Data (SWBD) [NASA/NGA, 2003] and HydroSHEDS flow direction map [Lerner16

et al., 2013] are used as the input datasets to develop GWD-LR.

17

The water body mask describes whether each pixel is a land pixel (value: 0) or a water

18

body pixel (value: 1). Types of water body (i.e. river, lake, ocean) do not have to be19

distinguished, and all water bodies should be indicated by the value 1 (Figure 1.1).

20
21
22

Figure 1.1 Input water body mask. Blue: water body pixels, white: land pixels. A part of theCongo River (17.4-18.2E, 0.4-1.0S) is shown as an example.

23

The flow direction map describes the downstream direction of surface water flow at each24

pixel toward one of the eight neighboring pixels (1: north, 2: northeast, 3: east, 4: southeast,25

5: south, 6: southwest, 7: west, 8: northwest). The river mouth is indicated by the value 0,

26

while ocean pixels are represented by the value -9.

27
28

29
30
31
32

The drainage area map (or flow accumulation map) describes the accumulative drainage
area of each pixel. The drainage are map can be calculated from the flow direction map.

Figure 1.2: (a) Input flow direction map (HydroSHEDS). Colors represent flow directions(1: north, 2: northeast, 3: east, 4: southeast, 5: south, 6: southwest, 7: west, 8: northwest).
(b) Input drainage area map.

2

33

1.2 Output Data

34

The algorithm automatically calculates bank-to-bank river width (Figure 1.3a) and35

effective river width excluding islands (Figure 1.3b) for all water bodies in the calculation36

domain. The procedures of the algorithm are explained in Section 2. The modified flow

37

direction map is also outputted as a by-product (Figure 1.4b). The river width is calculated

38

along the modified flow direction map, thus it is straightforward to use the river width

39

database as a topographic parameter of large-scale hydrodynamic models.

40
41
42

Figure 1.3: (a) Bank-to-bank river width. (b) Effective river width excluding islands. Water
bodies are shown gray, while islands are represented by black.

43
44
45
46
47
48

Figure 1.4: (a) Original flow direction map. Major streams are shown by black lines, whileassociated flows are shown by blue lines. (b) Modified flow direction map. Centerline
pixels (black lines) and perpendicular flows to centerlines (blue lines) are illustrated.
Islands are represented by dark green. Note that flow directions of limited pixels areshown to represent the difference between two panels.

3

49

2 Calculation Steps50

2.1 Modification of Water Body Mask

51
52

Some modifications of the input water body mask are needed due to discrepancies
between the input water body mask and the input flow direction map.

53

[Step 1.1] Pixels whose drainage area is larger than a threshold value (default: 100 km2)

54

are changed to water body pixels. This modification is needed because some rivers in the55

input flow direction map may run outside of water body areas of the input water body mask

56

(see Figure 2.1).

57

The water body pixels added by this procedure are used for the calculation of58

bank-to-bank river width, but they are excluded in the calculation of effective river width.

59
60

Figure 2.1: Water mask modification using drainage area.

61[Step 1.2] Gaps between water masks whose area is smaller than the threshold value

62(default: 1000 km2) are filled as island (Figure 2.2). Water body pixels and island pixels are

63

termed as “in-bank pixels” which are used for the calculation of bank-to-bank river width.

64
65Bank-to-bank width is calculated for all in-bank area (water body and island pixels), whileisland pixels are excluded when effective river width is calculated (see Section 2.10).

66

The default threshold for island filling is set to a relatively large value (1000 km2). All

67

island gaps except for very large ones (e.g. Ilha do Bananal in the Amazon River) are filled

68by this threshold. This large threshold is used to amalgamate bifurcated channels into one

69

merged channel, because the GWD-LR is mainly developed for application to large-scale

70

river models which cannot represent channel bifurcation.

4

71
72
73

Figure 2.2: Island gap filling. Gaps in water body mask whose area is smaller than 1000
2

km are filled as island pixels (gray).

74
75[Step 1.3] Water body pixels which represent very narrow channels (with 1-pixel width)76

are changed to land pixels (Figure 2.3). This modification is applied to improve the77

computational efficiency.

78

Most of these very nrrow channels are not represented in the original SWBD water mask,

79

given that the minimum channel width of the SWBD water mask is 183 m. Thus, most of80

them are introduced by the modification in Step 1.1.

81
82

Figure 2.3: Narrow channels with 1-pixel width are changed to land pixels.

5

83[Step 1.4] Pixels located within 1 km downstream from water bodies are changed to84

water body pixels. This modification is required because it is assumed in the proposed85

algorithm that river width is calculated for each “water body unit” which shares one outlet

86

pixel. If there are multiple outlet from one water body, the water body is treated as multiple

87

water body units, and river width is calculated separately for each water body unit. The

88

modification in downstream of a water body generates a shared outlet pixel for each water89

body unit (Figure 2.4).

90

The water body pixels added by this procedure are used for the calculation of91

bank-to-bank river width, but they are excluded when effective river width is calculated.

92
93

Figure 2.4: Modification on water body downstream.

94
95

2.2 Calculation of Water Body unit ID

96[Step 2] Water body and island pixels which share the same water body outlet (a water

97

body pixel with flow direction toward a land pixel) are treated as one water body unit.

98

Identical ID number is given to each water body unit.

6

99
100
101

Figure 2.5: Water body unit ID. Each color represents one water body unit. Island pixels

102

2.3 Calculation of Bank Distance

are represented by gray.

103[Step 3] Bank pixels are identified by searching eight neighboring pixels of each water

104

body pixel. If any of neighboring pixels is a land pixel, the water body pixels is identified as a

105bank pixel. Then, distance to the nearest bank pixel (hereafter “bank distance”) is calculated

106

for each in-bank pixel (Figure 2.6).

107
108
109

Figure 2.6: Bank distance. Distance to the nearest bank pixel is calculated for eachin-bank pixel.

7

110

Figure 2.7 is a schematic illustration of the definitions of bank distance, centerline

111

distance, and local distance between adjacent pixels. The bank distance is given as the112

distance to the nearest river bank pixel. The bank distance of the blue squared pixel is given

113

by the length of blue vector. The centerline distance is given as the normalized distance to114

the nearest centerline pixel (explained in Section 2.8). The centerline distance of the red

115

squared pixel is given by “the length of the red vector divided by the length of the black

116

vector” (Equation 2.2). The black vector represents the bank distance of the black squared

117

pixel, which is the nearest centerline pixel of the red squared pixel. The local distance118

between the two adjacent pixels (red and green squared pixels) is given by the length of the119green line.

120
121
122
123
124
125
126
127

Figure 2.7: Schematic illustration of the definitions of bank distancedistance

Db , centerline

Dc , and the local distance between adjacent pixels L . The bank distance of theblue squared pixel is given by the length of blue vector. The centerline distance of the redsquared pixel is given by “the length of the red vector (geometric centerline distance Dg )divided by the length of the black vector”. The local distance between the two adjacentpixels (red and green squared pixels) is given by the length of the green line.

2.4 Definition of Centerline Pixels

128[Step 4.1] Centerline pixels are defined by searching convex points in the bank distance

129

field. A pixel is judged to be a centerline pixel when the following two conditions are

130

satisfied: (1) the bank distance of the considered pixel is longer than the bank distance of six

131

or more neighboring pixels; and (2) the maximum gradient of bank distance between the132

considered pixel and its neighboring pixel is not larger than the threshold gradient (set to be

133

0.26, or ~tan(15deg) ). The gradient of bank distance is calculated by Equation (2.1):

134

Db x 

Dbj  Dbi
L(2.1),

8

135

where Db x is the gradient of bank distances, Dbi is the bank distance of the136

considered pixel i , Dbj is the bank distance of the neighboring pixel j , and L is the137

local distance between the centers of the pixels i and j . Note that the distance between

138

two point is calculated as a function of longitude and latitude assuming the earth ellipsoid in

139

this study, so that the difference of actual distance in a degree coordination at different

140

latitude is considered.

141

The first condition was introduced because only the upstream and downstream pixels on

142

the centerline may have a bank distance longer than that of the considered centerline pixel.

143

The second condition was introduced to remove spurious centerlines detected by the first144

condition. These spurious centerlines are caused by curvatures of river banks (see blue

145

lines in Figure 2.8), and tend to extend from a river bank toward a true centerline. Thus, the146

gradient of bank distance tends to be larger on a spurious centerline than on a true147

centerline. The threshold gradient to distinguish true and spurious centerlines was set to be148

0.26 ( ~tan(15deg) ) by trial and error. We found that smaller threshold gradient produces less149

spurious centerlines, but some true centerlines are not detected when the threshold gradient150

is too small. The centerline pixels determined by these conditions are shown in gray lines in151

Figure 2.8, while spurious centerline pixels are shown in blue dots.

152
153
154
155

Figure 2.8: Defining centerline pixels by searching convex points in bank distance field.
Gray: defined centerline pixels, blue: spurious centerline pixels. The bank distance field
is shown by the background colors.

9

156

[Step 4.2] The centerline pixels detected by Step 4.1 have large gaps between them157

where river width is increasing because of the second condition. In order to improve the158

connectivity of centerlines, centerline pixels are extended by the following procedures: (1)159

for each centerline pixel, the pixel with maximum bank distance among eight neighboring160

pixels is selected; (2) if the selected neighboring pixel has a larger bank distance than the161

considered centerline pixel, the selected neighboring pixel is converted to a new centerline

162

pixel; (3) the extension procedure is repeated during the criteria (2) is true. The extended

163

centerline pixels are shown by red lines in Figure 2.9.

164
165
166

Figure 2.9: Centerline extension. Gray: original centerline pixels, red: extended centerline
pixels. The bank distance field is shown by the background colors.

167
168

2.5 Calculation of Outlet Distance169[Step 5] The riverline distance from the outlet pixel of each water body unit (hereafter

170“outlet distance”) is calculated. The outlet distance is calculated by accumulating the local

171

distance between pixels from the outlet pixel toward upstream (Figure 2.10).

172

In order to connect centerline pixels in a sequential downstream direction, the local

173

distance between adjacent pixels is weighted by 10 when the pixel in the upstream side is174

not a centerline pixel. (blue vectors in Figure 2.10). Sensitivity to this weighing parameter is

175

discussed in Section 3.

10

176

177
178
179
180
181
182
183
184

The calculated outlet distance is shown in Figure 2.11.

Figure 2.10: Schematic illustration of the outlet distance calculation. The pixel withnumber 0 denotes the outlet pixel of the water body unit. For easy explanation, the local
distance between adjacent pixels in an orthogonal position is set to 1, while the local
distance between adjacent pixels in a diagonal position is set to 1.4. The weighted local
distances are set to 10 for orthogonal direction and 14 for diagonal direction (bluevectors). The outlet distance is calculated cumulatively from the outlet pixel of a water
body unit toward upstream.

185
186
187

Figure 2.11: Outlet distance. The cumulative riverline distance with weight for
non-centerline pixels from the outlet of each water body unit is shown by periodic colors.

11

188

2.6 Setting Centerline Flow Directions

189[Step 6.1] The downstream direction of each centerline pixel is determined by choosing

190

the pixel with minimum outlet distance among its eight neighboring pixels (Figure 2.12). If

191

the selected downstream pixel is not a centerline pixel, the selected downstream pixel is192

changed to a new centerline pixel. The centerline extension is repeated until the193

downstream centerline pixel or the outlet pixel is connected (connected centerline pixels are

194

shown by orange boxes in Figure 2.12 and red lines in Figure 2.13). Thus, the connectivity

195

between all centerline pixels is ensured within a water body unit (Figure 2.13).

196
197

Figure 2.12: Determination of centerline flow directions.

198
199
200

Figure 2.13: Centerline connection. Gray: original centerline pixels, red: connected
centerline pixels. The outlet distance field is shown by the background colors.

12

201

[Step 6.2] The topmost centerline pixels are modified to non-centerline pixels. Most of the202

topmost centerline pixels which do not have an upstream centerline pixel are caused by the203

determination of flow directions (Figure 2.14) and they don’t have to be treated as centerline

204

pixels any more.

205

206
207
208

Figure 2.14: Modification on topmost centerline pixels. Topmost centerline pixels (orange)
are converted to non-centerline.

209
210

2.7 Calculation of bank-to-bank river width

211[Step 7.1] Centerline pixels are classified into actual centerlines and virtual connecting

212

centerlines which connect a small tributary to a wide main channel. Centerline pixels are

213

considered to be actual centerline when the gradient of bank distances toward downstream

214

centerline pixel (calculated by Equation 2.1) is not larger than 0.57 ( ~tan(30deg) ), while it’s215

considered to be virtual connecting centerlines when the gradient is larger than 0.57 (Figure

2162.15). Virtual connecting centerlines have rapidly increasing bank distance because they217

run from bank toward centerline in a wide main channel (blue lines in Figure 2.15). The

218

threshold value of 0. 57 [ ~tan(30deg) ] is decided by trial and error.

13

219
220
221

Figure 2.15: Centerline classification. Gray: actual centerline. Blue: virtual connecting

222[Step 7.2] The bank-to-bank river width of actual centerline pixels is set to twice of its

223

bank distance. The bank-to-bank distance of the virtual connecting centerlines (blue lines in224

Figure 2.15) is set to the same value as the bank-to-bank river width of its nearest upstream

225

actual centerline pixel. The calculated bank-to-bank river width is shown in Figure 2.16.

226
227

Figure 2.16 Bank-to-bank river width. Gray: non-centerline water body. Black: island.

centerline. Background color represents bank distance.

14

228

2.8 Calculation of Centerline Distance229[Step 8] For every in-bank pixels, the distance to its nearest centerline pixel (hereafter

230“centerline distance”) is calculated. Centerline distance is normalized by the bank distance231of the nearest centerline pixel in order to avoid unrealistic accumulation of flow from an area

232

outside of a tributary’s width, within the zone where the tributary merges into its main

233

channel (Explained in Section 3.2).

234

The scaled centerline distance is given by Equation (2.2):

235

Dc 

236

where Dc is the normalized centerline distance, D g is the geometric distance to the237

centerline pixel, and Dbc is the bank distance of the nearest centerline pixel. The

238

normalized centerline distance of a considered pixel becomes smaller when its nearest239

centerline pixel has larger bank distance and vice versa. The weighted centerline distance is240

shown by colors in Figure 2.17.

Dg
Dbc(2.2),

241
242
243

Figure 2.17: Centerline distance. The background colors represent centerline distance (i.e.
normalized distance to its nearest centerline). Centerlines are illustrated by black lines.

15

244
245
246

2.9 Determination of flow direction[Step 9] The flow direction of each non-centerline pixel is decided based on the gradientof centerline distance given by Equation 2.3.

Dcj  Dci

247

Dc x 

248

where Dc x is the gradient of centerline distance, Dci is the centerline distance of249

the considered pixel i , Dcj is the centerline distance of the neighboring pixel j , and L

250

is the local distance between the pixels i and j .

L(2.3),

251

Then, the flow direction is determined for each pixel by choosing the maximum gradient252

among the eight neighboring pixels. The modified flow directions are shown in Figure 2.18.

253
254
255
256

Figure 2.18 Calculation of Flow Directions. Flow directions are represented by blue lines
while centerlines are indicated by black lines. Background colors represent normalized
centerline distance. Note that flow directions of limited pixels are shown.

16

257

2.10 Calculation of effective river width

258[Step 10.1] The effective river segment for each centerline pixel is defined for the259

calculation of effective river width excluding islands. The effective river segment for a given

260

centerline pixel is defined as; the longitudinal reach within the bank distance length either

261

side of the pixel (the red dashed line in Figure 2.19), and includes the none-centerline pixels

262

draining to that segment of centerline. (the area shaded with green in Figure 2.19).

263
264
265
266
267

Figure 2.19: Calculation of effective channel segment. The green shaded area represents
the effective river segment of the centerline pixel marked by the yellow square. The
dashed red vector represents the effective centerline of the considered centerline pixel,
which is determined by the bank distance of the considered pixels (orange vectors).

268
269
270

17

271

[Step 10.2] The total in-water area (water body and island) and the total water body area272

within the effective river segment are calculated. Then, effective river width is calculated by273

Equation (2.4).

274

A 
We  Wb   w 
 At 

275

where We is the effective river width, Wb is the bank-to-bank river width, At is the276(2.4),

total in-water area (i.e. water body and islands), Aw is the total water body area.

277

The effective river width of the virtual connecting centerlines (blue lines in Figure 2.15) is278

modified to the same value as the effective river width of its nearest upstream actual279

centerline pixel.

280

281
282
283

The calculated effective river width is shown in Figure 2.20.

Figure 2.20. Effective river width. Non-centerline water bodies are represented by gray,while islands are represented by black.

284

18

285
286

3. Sensitivity to Parameters
Sensitivity of the river width calculation to parameters used in the global river width287

algorithm is discussed in this section.

288

3.1 Weight for outlet distance289

The weight on the outlet distance is introduced to achieve sequential downstream

290

connection of centerline pixels (Section 2.6). If the weight is too small, centerline pixels are291

not appropriately connected (the weight is set to 1 and 3 in Figures 3.1a and 3.1b,292

respectively). Intermittent centerlines are preferred to be connected to their nearest

293

downstream centerline, but it’s difficult to find the nearest downstream centerline when the294weight too small. When the weight is large enough, the result is not sensitive to the weight295

value (the weight is set to 10 and 30 in Figures 3.1c and 3.1d, respectively). Calculation time296

becomes longer when a larger weight is used, thus we decided to use the weight of 10.

297
298
299
300
301

Figure 3.1: Sensitivity of centerline connection to the weight on outlet distancecalculation. The original centerlines calculated from bank distance are shown by graylines, while connected centerlines calculated by outlet distance are shown by red lines.
Background colors represent outlet distance.

19

302

3.2 Normalized centerline distance

303

Centerline distance (i.e. distance to the nearest centerline, see Section 2.8) is normalized304

by the bank distance of the nearest centerline pixel (see, Equation 2.2). The normalization is305

introduced in order to avoid unrealistic accumulation of flow from an area outside of a306

tributary’s width, within the zone where the tributary merges into its main channel. If

307

centerline distance is not normalized, small tributaries unrealistically gather flows from its

308

main channel (Figure 3.2b).

309

310
311
312
313
314

Figure 3.2 Modified flow direction derived from (a) normalized centerline distance, and (b)geometric (non-normalized) centerline distance. Flow directions are shown by blue lines,
while centerlines are illustrated by black lines. Background colors represent centerline
distance.

315
316
317
318

3.3 Threshold for island gap filling

319

Relatively large threshold (1000 km2) for island gap filling was used in order to320

amalgamate bifurcated channels into one merged channel (Figure 3.3a), because the main321

target of GWD-LR is application to large-scale hydrodynamic models which cannot

322

represent channel bifurcations.

20

323

In case width values are needed for each bifurcated channels, smaller threshold value for324

island gap filling may be preferred. By taking the smaller threshold, effective river width can325

be calculated separately for bifurcated channels (Figure 3.3b). Note that accuracy of326

effective river width calculation for bifurcated channels depend on the quality of the input

327

water mask and input flow direction map. Manual correction may be needed because small

328

channels are generally not well represented both in the water mask and the flow direction

329

map.

330
331
332
333

Figure 3.3: Effective river width calculated with a different threshold value for island gap
2

2

filling. (a) 1000 km threshold, (b) 10 km threshold. Water masks are shown by gray, while
filled islands are shown by black.

334
335

References

336
337
338

Lehner, B, and G. Grill (2013), Global river hydrography and network routing: baseline dataand new approaches to study the world's large river systems, Hydrol. Proc., 27,
2171-2186, doi:10.1002/hyp.9740.

339
340

NASA/NGA (2003), SRTM Water Body Data Product Specific Guidance, Version 2.0,
available online: http://dds.cr.usgs.gov/srtm/version2_1/SWBD/SWBD_Documentation/

341
342
343

Yamazaki D., F. O’Loughlin, M. A. Trigg, Z. F. Miller, T. M. Pavelsky, and P. D. Bates (2014),
Development of the Global Width Database for Large Rivers, Water Resources
Research, vol.50, in print, DOI: 10.1002/2013WR014664

21