| 1 | #!/usr/bin/env python
|
|---|
| 2 | # -*- coding: utf-8
|
|---|
| 3 |
|
|---|
| 4 | """
|
|---|
| 5 | MODULE: r.valley.bottom
|
|---|
| 6 |
|
|---|
| 7 | AUTHOR(S): Helmut Kudrnovsky <alectoria AT gmx at>
|
|---|
| 8 | Steven Pawley <dr.stevenpawley@gmail.com>
|
|---|
| 9 |
|
|---|
| 10 | PURPOSE: Calculates Multi-resolution Valley Bottom Flatness (MrVBF) index
|
|---|
| 11 | Index is inspired by
|
|---|
| 12 | John C. Gallant and Trevor I. Dowling 2003.
|
|---|
| 13 | A multiresolution index of valley bottom flatness for mapping depositional areas.
|
|---|
| 14 | WATER RESOURCES RESEARCH, VOL. 39, NO. 12, 1347, doi:10.1029/2002WR001426, 2003
|
|---|
| 15 |
|
|---|
| 16 | COPYRIGHT: (C) 2018 by the GRASS Development Team
|
|---|
| 17 |
|
|---|
| 18 | This program is free software under the GNU General Public
|
|---|
| 19 | License (>=v2). Read the file COPYING that comes with GRASS
|
|---|
| 20 | for details.
|
|---|
| 21 | """
|
|---|
| 22 |
|
|---|
| 23 | #%module
|
|---|
| 24 | #% description: Calculation of Multi-resolution Valley Bottom Flatness (MrVBF) index
|
|---|
| 25 | #% keyword: raster
|
|---|
| 26 | #% keyword: terrain
|
|---|
| 27 | #%end
|
|---|
| 28 |
|
|---|
| 29 | #%option G_OPT_R_ELEV
|
|---|
| 30 | #% key: elevation
|
|---|
| 31 | #% description: Name of elevation raster map
|
|---|
| 32 | #% required: yes
|
|---|
| 33 | #%end
|
|---|
| 34 |
|
|---|
| 35 | #%option G_OPT_R_OUTPUT
|
|---|
| 36 | #% key: mrvbf
|
|---|
| 37 | #% description: Name of output MRVBF raster map
|
|---|
| 38 | #% required: yes
|
|---|
| 39 | #%end
|
|---|
| 40 |
|
|---|
| 41 | #%option G_OPT_R_OUTPUT
|
|---|
| 42 | #% key: mrrtf
|
|---|
| 43 | #% description: Name of output MRRTF raster map
|
|---|
| 44 | #% required: no
|
|---|
| 45 | #%end
|
|---|
| 46 |
|
|---|
| 47 | #%option
|
|---|
| 48 | #% key: t_slope
|
|---|
| 49 | #% description: Initial Threshold for Slope
|
|---|
| 50 | #% required: no
|
|---|
| 51 | #% answer: 16
|
|---|
| 52 | #% end
|
|---|
| 53 |
|
|---|
| 54 | #%option
|
|---|
| 55 | #% key: t_pctl_v
|
|---|
| 56 | #% description: Threshold (t) for transformation of Elevation Percentile (Lowness)
|
|---|
| 57 | #% required: no
|
|---|
| 58 | #% answer: 0.4
|
|---|
| 59 | #% end
|
|---|
| 60 |
|
|---|
| 61 | #%option
|
|---|
| 62 | #% key: t_pctl_r
|
|---|
| 63 | #% description: Threshold (t) for transformation of Elevation Percentile (Upness)
|
|---|
| 64 | #% required: no
|
|---|
| 65 | #% answer: 0.3
|
|---|
| 66 | #% end
|
|---|
| 67 |
|
|---|
| 68 | #%option
|
|---|
| 69 | #% key: t_vf
|
|---|
| 70 | #% description: Threshold (t) for transformation of Valley Bottom Flatness
|
|---|
| 71 | #% required: no
|
|---|
| 72 | #% answer: 0.3
|
|---|
| 73 | #% end
|
|---|
| 74 |
|
|---|
| 75 | #%option
|
|---|
| 76 | #% key: t_rf
|
|---|
| 77 | #% description: Threshold (t) for transformation of Ridge Top Flatness
|
|---|
| 78 | #% required: no
|
|---|
| 79 | #% answer: 0.35
|
|---|
| 80 | #% end
|
|---|
| 81 |
|
|---|
| 82 | #%option
|
|---|
| 83 | #% key: p_slope
|
|---|
| 84 | #% description: Shape Parameter (p) for Slope
|
|---|
| 85 | #% required: no
|
|---|
| 86 | #% answer: 4
|
|---|
| 87 | #% end
|
|---|
| 88 |
|
|---|
| 89 | #%option
|
|---|
| 90 | #% key: p_pctl
|
|---|
| 91 | #% description: Shape Parameter (p) for Elevation Percentile
|
|---|
| 92 | #% required: no
|
|---|
| 93 | #% answer: 3
|
|---|
| 94 | #% end
|
|---|
| 95 |
|
|---|
| 96 | #%option
|
|---|
| 97 | #% key: min_cells
|
|---|
| 98 | #% description: Minimum number of cells in generalized DEM
|
|---|
| 99 | #% required: no
|
|---|
| 100 | #% answer: 1
|
|---|
| 101 | #% end
|
|---|
| 102 |
|
|---|
| 103 | #%flag
|
|---|
| 104 | #% key: s
|
|---|
| 105 | #% description: Use square moving window instead of circular moving window
|
|---|
| 106 | #%end
|
|---|
| 107 |
|
|---|
| 108 |
|
|---|
| 109 | from grass.pygrass.gis.region import Region
|
|---|
| 110 | from grass.pygrass.modules.shortcuts import general as g
|
|---|
| 111 | from grass.pygrass.modules.shortcuts import raster as r
|
|---|
| 112 | import sys
|
|---|
| 113 | import os
|
|---|
| 114 | import grass.script as grass
|
|---|
| 115 | import math
|
|---|
| 116 | import atexit
|
|---|
| 117 | import random
|
|---|
| 118 | import string
|
|---|
| 119 |
|
|---|
| 120 | if "GISBASE" not in os.environ:
|
|---|
| 121 | print("You must be in GRASS GIS to run this program.")
|
|---|
| 122 | sys.exit(1)
|
|---|
| 123 |
|
|---|
| 124 |
|
|---|
| 125 | def cleanup():
|
|---|
| 126 | grass.message("Deleting intermediate files...")
|
|---|
| 127 |
|
|---|
| 128 | for k, v in TMP_RAST.iteritems():
|
|---|
| 129 | for f in v:
|
|---|
| 130 | if len(grass.find_file(f)['fullname']) > 0:
|
|---|
| 131 | grass.run_command(
|
|---|
| 132 | "g.remove", type="raster", name=f, flags="f", quiet=True)
|
|---|
| 133 |
|
|---|
| 134 | Region.write(current_region)
|
|---|
| 135 |
|
|---|
| 136 |
|
|---|
| 137 | def cell_padding(input, output, radius=3):
|
|---|
| 138 | """Mitigates edge effect by growing an input raster map by radius cells
|
|---|
| 139 |
|
|---|
| 140 | Args
|
|---|
| 141 | ----
|
|---|
| 142 | input, output : str
|
|---|
| 143 | Names of GRASS raster map for input, and padded output
|
|---|
| 144 | radius : int
|
|---|
| 145 | Radius in which to expand region and grow raster
|
|---|
| 146 |
|
|---|
| 147 | Returns
|
|---|
| 148 | -------
|
|---|
| 149 | input_grown : str
|
|---|
| 150 | GRASS raster map which has been expanded by radius cells"""
|
|---|
| 151 |
|
|---|
| 152 | region = Region()
|
|---|
| 153 |
|
|---|
| 154 | g.region(n=region.north + (region.nsres * radius),
|
|---|
| 155 | s=region.south - (region.nsres * radius),
|
|---|
| 156 | w=region.west - (region.ewres * radius),
|
|---|
| 157 | e=region.east + (region.ewres * radius))
|
|---|
| 158 |
|
|---|
| 159 | r.grow(input=input,
|
|---|
| 160 | output=output,
|
|---|
| 161 | radius=radius+1,
|
|---|
| 162 | quiet=True)
|
|---|
| 163 |
|
|---|
| 164 | return (region)
|
|---|
| 165 |
|
|---|
| 166 |
|
|---|
| 167 | def focal_expr(radius, window_square=False):
|
|---|
| 168 | """Returns array offsets relative to centre pixel (0,0) for a matrix of
|
|---|
| 169 | size radius
|
|---|
| 170 |
|
|---|
| 171 | Args
|
|---|
| 172 | ----
|
|---|
| 173 | radius : int
|
|---|
| 174 | Radius of the focal function
|
|---|
| 175 | window_square : bool. Optional (default is False)
|
|---|
| 176 | Boolean to use a circular or square window
|
|---|
| 177 |
|
|---|
| 178 | Returns
|
|---|
| 179 | -------
|
|---|
| 180 | offsets : list
|
|---|
| 181 | List of pixel positions (row, col) relative to the center pixel
|
|---|
| 182 | ( 1, -1) ( 1, 0) ( 1, 1)
|
|---|
| 183 | ( 0, -1) ( 0, 0) ( 0, 1)
|
|---|
| 184 | (-1, -1) (-1, 0) (-1, 1)"""
|
|---|
| 185 |
|
|---|
| 186 | offsets = []
|
|---|
| 187 |
|
|---|
| 188 | # generate a list of spatial neighbourhood offsets for the chosen radius
|
|---|
| 189 | # ignoring the centre cell
|
|---|
| 190 | if window_square:
|
|---|
| 191 |
|
|---|
| 192 | for i in range(-radius, radius+1):
|
|---|
| 193 | for j in range(-radius, radius+1):
|
|---|
| 194 | if (i,j) != (0,0):
|
|---|
| 195 | offsets.append((i, j))
|
|---|
| 196 |
|
|---|
| 197 | else:
|
|---|
| 198 |
|
|---|
| 199 | for i in range(-radius, radius+1):
|
|---|
| 200 | for j in range(-radius, radius+1):
|
|---|
| 201 | row = i + radius
|
|---|
| 202 | col = j + radius
|
|---|
| 203 |
|
|---|
| 204 | if pow(row - radius, 2) + pow(col - radius, 2) <= \
|
|---|
| 205 | pow(radius, 2) and (i, j) != (0,0):
|
|---|
| 206 | offsets.append((j, i))
|
|---|
| 207 |
|
|---|
| 208 | return offsets
|
|---|
| 209 |
|
|---|
| 210 |
|
|---|
| 211 | def get_percentile(L, input, radius=3, window_square=False):
|
|---|
| 212 | """Calculates the percentile whichj is the ratio of the number of points of
|
|---|
| 213 | lower elevation to the total number of points in the surrounding region
|
|---|
| 214 |
|
|---|
| 215 | Args
|
|---|
| 216 | ----
|
|---|
| 217 | L : int
|
|---|
| 218 | Processing step (level)
|
|---|
| 219 | input : str
|
|---|
| 220 | GRASS raster map (elevation) to perform calculation on
|
|---|
| 221 | radius : int
|
|---|
| 222 | Neighborhood radius (in pixels)
|
|---|
| 223 | window_square : bool. Optional (default is False)
|
|---|
| 224 | Boolean to use square or circular neighborhood
|
|---|
| 225 |
|
|---|
| 226 | Returns
|
|---|
| 227 | -------
|
|---|
| 228 | PCTL : str
|
|---|
| 229 | Name of GRASS cell-padded raster map with elevation percentile for
|
|---|
| 230 | processing step L"""
|
|---|
| 231 |
|
|---|
| 232 | PCTL = "PCTL{L}".format(L=L+1)
|
|---|
| 233 | TMP_RAST[L].append(PCTL)
|
|---|
| 234 | input_grown=input
|
|---|
| 235 |
|
|---|
| 236 | # get offsets for given neighborhood radius
|
|---|
| 237 | offsets = focal_expr(radius=radius, window_square=window_square)
|
|---|
| 238 |
|
|---|
| 239 | # generate grass mapcalc terms and execute
|
|---|
| 240 | n_pixels = float(len(offsets))
|
|---|
| 241 |
|
|---|
| 242 | # create mapcalc expr
|
|---|
| 243 | # if pixel in neighborhood contains nodata, attempt to use opposite neighbor
|
|---|
| 244 | # if opposite neighbor is also nodata, then use center pixel
|
|---|
| 245 | terms = []
|
|---|
| 246 | for d in offsets:
|
|---|
| 247 | valid = ','.join(map(str, d))
|
|---|
| 248 | invalid = ','.join([str(-d[0]), str(-d[1])])
|
|---|
| 249 | terms.append("if(isnull({input}[{d}]), if(isnull({input}[{e}]), 1, {input}[{e}]<={input}), {input}[{d}]<={input})".format(
|
|---|
| 250 | input=input_grown, d=valid, e=invalid))
|
|---|
| 251 |
|
|---|
| 252 | expr = "{x} = ({s}) / {n}".format(x=PCTL, s=" + ".join(terms), n=n_pixels)
|
|---|
| 253 | grass.mapcalc(expr)
|
|---|
| 254 |
|
|---|
| 255 | return(PCTL)
|
|---|
| 256 |
|
|---|
| 257 |
|
|---|
| 258 | def get_slope(L, elevation):
|
|---|
| 259 |
|
|---|
| 260 | region = Region()
|
|---|
| 261 |
|
|---|
| 262 | slope = 'tmp_slope_step{L}'.format(L=L+1) + \
|
|---|
| 263 | ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
|
|---|
| 264 | TMP_RAST[L].append(slope)
|
|---|
| 265 |
|
|---|
| 266 | slope_expr = """{s} = 100 * (sqrt( \
|
|---|
| 267 | ((if(isnull({a}[0,-1]), {a}[0,0], {a}[0,-1]) - if(isnull({a}[0, 1]), {a}[0,0], {a}[0, 1])) / (2*{ewres}))^2 + \
|
|---|
| 268 | ((if(isnull({a}[-1,0]), {a}[0,0], {a}[-1,0]) - if( isnull({a}[1, 0]), {a}[0,0], {a}[1, 0])) / (2*{nsres}))^2)) \
|
|---|
| 269 | """.format(
|
|---|
| 270 | s=slope, a=elevation, nsres=region.nsres, ewres=region.ewres)
|
|---|
| 271 |
|
|---|
| 272 | grass.mapcalc(slope_expr)
|
|---|
| 273 |
|
|---|
| 274 | return (slope)
|
|---|
| 275 |
|
|---|
| 276 |
|
|---|
| 277 | def get_flatness(L, slope, t, p):
|
|---|
| 278 | """Calculates the flatness index
|
|---|
| 279 | Flatness F1 = 1 / (1 + pow ((slope / t), p)
|
|---|
| 280 |
|
|---|
| 281 | Equation 2 (Gallant and Dowling, 2003)"""
|
|---|
| 282 |
|
|---|
| 283 | F = "tmp_F{L}".format(L=L+1) + \
|
|---|
| 284 | ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
|
|---|
| 285 | TMP_RAST[L].append(F)
|
|---|
| 286 | grass.mapcalc("$g = 1.0 / (1.0 + pow(($x / $t), $p))",
|
|---|
| 287 | g=F, x=slope, t=t, p=p)
|
|---|
| 288 |
|
|---|
| 289 | return F
|
|---|
| 290 |
|
|---|
| 291 |
|
|---|
| 292 | def get_prelim_flatness(L, F, PCTL, t, p):
|
|---|
| 293 | """Transform elevation percentile to a local lowness value using
|
|---|
| 294 | equation (1) and combined with flatness F to produce the preliminary
|
|---|
| 295 | valley flatness index (PVF) for the first step.
|
|---|
| 296 |
|
|---|
| 297 | Equation 3 (Gallant and Dowling, 2003)"""
|
|---|
| 298 |
|
|---|
| 299 | PVF = "tmp_PVF{L}" + \
|
|---|
| 300 | ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
|
|---|
| 301 | TMP_RAST[L].append(PVF)
|
|---|
| 302 |
|
|---|
| 303 | grass.mapcalc("$g = $a * (1.0 / (1.0 + pow(($x / $t), $p)))",
|
|---|
| 304 | g=PVF, a=F, x=PCTL, t=t, p=p)
|
|---|
| 305 |
|
|---|
| 306 | return PVF
|
|---|
| 307 |
|
|---|
| 308 |
|
|---|
| 309 | def get_prelim_flatness_rf(L, F, PCTL, t, p):
|
|---|
| 310 | """Transform elevation percentile to a local upness value using
|
|---|
| 311 | equation (1) and combined with flatness to produce the preliminary
|
|---|
| 312 | valley flatness index (PVF) for the first step
|
|---|
| 313 |
|
|---|
| 314 | Equation 3 (Gallant and Dowling, 2003)"""
|
|---|
| 315 |
|
|---|
| 316 | PVF = "tmp_PVF{L}".format(L=L+1) + \
|
|---|
| 317 | ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
|
|---|
| 318 | TMP_RAST[L].append(PVF)
|
|---|
| 319 |
|
|---|
| 320 | grass.mapcalc("$g = $a * (1.0 / (1.0 + pow(((1-$x) / $t), $p)))",
|
|---|
| 321 | g=PVF, a=F, x=PCTL, t=t, p=p)
|
|---|
| 322 |
|
|---|
| 323 | return PVF
|
|---|
| 324 |
|
|---|
| 325 |
|
|---|
| 326 | def get_valley_flatness(L, PVF, t, p):
|
|---|
| 327 | """Calculation of the valley flatness step VF
|
|---|
| 328 | Larger values of VF1 indicate increasing valley bottom character
|
|---|
| 329 | with values less than 0.5 considered not to be in valley bottoms
|
|---|
| 330 |
|
|---|
| 331 | Equation 4 (Gallant and Dowling, 2003)"""
|
|---|
| 332 |
|
|---|
| 333 | VF = "tmp_VF{L}".format(L=L+1) + \
|
|---|
| 334 | ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
|
|---|
| 335 | TMP_RAST[L].append(VF)
|
|---|
| 336 | grass.mapcalc("$g = 1 - (1.0 / (1.0 + pow(($x / $t), $p)))",
|
|---|
| 337 | g=VF, x=PVF, t=t, p=p)
|
|---|
| 338 |
|
|---|
| 339 | return VF
|
|---|
| 340 |
|
|---|
| 341 |
|
|---|
| 342 | def get_mrvbf(L, VF_Lminus1, VF_L, t):
|
|---|
| 343 | """Calculation of the MRVBF index
|
|---|
| 344 | Requires that L>1"""
|
|---|
| 345 |
|
|---|
| 346 | W = "tmp_W{L}".format(L=L+1) + \
|
|---|
| 347 | ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
|
|---|
| 348 | MRVBF = "MRVBF{L}".format(L=L+1) + \
|
|---|
| 349 | ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
|
|---|
| 350 |
|
|---|
| 351 | # Calculation of weight W2 (Equation 9)
|
|---|
| 352 | TMP_RAST[L].append(W)
|
|---|
| 353 | p = (math.log10(((L+1) - 0.5) / 0.1)) / math.log10(1.5)
|
|---|
| 354 | grass.mapcalc("$g = 1 - (1.0 / (1.0 + pow(($x / $t), $p)))",
|
|---|
| 355 | g=W, x=VF_L, t=t, p=p)
|
|---|
| 356 |
|
|---|
| 357 | # Calculation of MRVBF2 (Equation 8)
|
|---|
| 358 | TMP_RAST[L].append(MRVBF)
|
|---|
| 359 | grass.mapcalc("$MBF = ($W * ($L + $VF)) + ((1 - $W) * $VF1)",
|
|---|
| 360 | MBF=MRVBF, L=L, W=W, VF=VF_L, VF1=VF_Lminus1)
|
|---|
| 361 |
|
|---|
| 362 | return MRVBF
|
|---|
| 363 |
|
|---|
| 364 |
|
|---|
| 365 | def get_combined_flatness(L, F1, F2):
|
|---|
| 366 | """Calculates the combined flatness index
|
|---|
| 367 |
|
|---|
| 368 | Equation 13 (Gallant and Dowling, 2003)"""
|
|---|
| 369 |
|
|---|
| 370 | CF = "tmp_CF{L}".format(L=L+1) + \
|
|---|
| 371 | ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
|
|---|
| 372 | TMP_RAST[L].append(CF)
|
|---|
| 373 | grass.mapcalc("$CF = $F1 * $F2", CF=CF, F1=F1, F2=F2)
|
|---|
| 374 |
|
|---|
| 375 | return CF
|
|---|
| 376 |
|
|---|
| 377 |
|
|---|
| 378 | def get_smoothed_dem(L, DEM):
|
|---|
| 379 | """Smooth the DEM using an 11 cell averaging filter with gauss weighting of 3 radius"""
|
|---|
| 380 |
|
|---|
| 381 | smoothed = "tmp_DEM_smoothed_step{L}".format(L=L+1) + \
|
|---|
| 382 | ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
|
|---|
| 383 | TMP_RAST[L].append(smoothed)
|
|---|
| 384 |
|
|---|
| 385 | # g = 4.3565 * math.exp(-(3 / 3.0))
|
|---|
| 386 | r.neighbors(input=DEM, output=smoothed, size=11, gauss=3)
|
|---|
| 387 |
|
|---|
| 388 | return smoothed
|
|---|
| 389 |
|
|---|
| 390 |
|
|---|
| 391 | def refine(L, input, region, method='bilinear'):
|
|---|
| 392 | """change resolution back to base resolution and resample a raster"""
|
|---|
| 393 |
|
|---|
| 394 | input_padded = '_'.join(['tmp', input, '_padded']) + \
|
|---|
| 395 | ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
|
|---|
| 396 | TMP_RAST[L].append(input_padded)
|
|---|
| 397 | cell_padding(input=input, output=input_padded, radius=2)
|
|---|
| 398 |
|
|---|
| 399 | Region.write(region)
|
|---|
| 400 | input = '_'.join(['tmp', input, 'refined_base_resolution']) + \
|
|---|
| 401 | ''.join([random.choice(string.ascii_letters + string.digits) for n in range(4)])
|
|---|
| 402 | TMP_RAST[L].append(input)
|
|---|
| 403 |
|
|---|
| 404 | if method == 'bilinear':
|
|---|
| 405 | r.resamp_interp(input=input_padded,
|
|---|
| 406 | output=input,
|
|---|
| 407 | method="bilinear")
|
|---|
| 408 |
|
|---|
| 409 | if method == 'average':
|
|---|
| 410 | r.resamp_stats(input=input_padded,
|
|---|
| 411 | output=input,
|
|---|
| 412 | method='average',
|
|---|
| 413 | flags='w')
|
|---|
| 414 |
|
|---|
| 415 | return input
|
|---|
| 416 |
|
|---|
| 417 |
|
|---|
| 418 | def main():
|
|---|
| 419 | r_elevation = options['elevation']
|
|---|
| 420 | mrvbf = options['mrvbf'].split('@')[0]
|
|---|
| 421 | mrrtf = options['mrrtf'].split('@')[0]
|
|---|
| 422 | t_slope = float(options['t_slope'])
|
|---|
| 423 | t_pctl_v = float(options['t_pctl_v'])
|
|---|
| 424 | t_pctl_r = float(options['t_pctl_r'])
|
|---|
| 425 | t_vf = float(options['t_rf'])
|
|---|
| 426 | t_rf = float(options['t_rf'])
|
|---|
| 427 | p_slope = float(options['p_slope'])
|
|---|
| 428 | p_pctl = float(options['p_pctl'])
|
|---|
| 429 | moving_window_square = flags['s']
|
|---|
| 430 | min_cells = int(options['min_cells'])
|
|---|
| 431 |
|
|---|
| 432 | global current_region
|
|---|
| 433 | global TMP_RAST
|
|---|
| 434 | TMP_RAST = {}
|
|---|
| 435 | current_region = Region()
|
|---|
| 436 |
|
|---|
| 437 | # Some checks
|
|---|
| 438 | if t_slope <= 0 or t_pctl_v <= 0 or t_pctl_r <= 0 or t_vf <= 0 or t_rf <= 0 or \
|
|---|
| 439 | p_slope <= 0 or p_pctl <= 0:
|
|---|
| 440 | grass.fatal('Parameter values cannot be <= 0')
|
|---|
| 441 |
|
|---|
| 442 | if min_cells < 1:
|
|---|
| 443 | grass.fatal('Minimum number of cells in generalized DEM cannot be less than 1')
|
|---|
| 444 |
|
|---|
| 445 | if min_cells > current_region.cells:
|
|---|
| 446 | grass.fatal('Minimum number of cells in the generalized DEM cannot exceed the ungeneralized number of cells')
|
|---|
| 447 |
|
|---|
| 448 | ###########################################################################
|
|---|
| 449 | # Calculate number of levels
|
|---|
| 450 |
|
|---|
| 451 | levels = math.ceil(-math.log(float(min_cells)/current_region.cells) / math.log(3) - 2)
|
|---|
| 452 | levels = int(levels)
|
|---|
| 453 |
|
|---|
| 454 | if levels < 3:
|
|---|
| 455 | grass.fatal('MRVBF algorithm requires a greater level of generalization. Reduce number of min_cells')
|
|---|
| 456 |
|
|---|
| 457 | grass.message('Parameter Settings')
|
|---|
| 458 | grass.message('------------------')
|
|---|
| 459 | grass.message('min_cells = %d will result in %d generalization steps' % (min_cells, levels))
|
|---|
| 460 |
|
|---|
| 461 | TMP_RAST = {k: [] for k in range(levels)}
|
|---|
| 462 |
|
|---|
| 463 | ###########################################################################
|
|---|
| 464 | # Intermediate outputs
|
|---|
| 465 | Xres_step, Yres_step, DEM = [], [], []
|
|---|
| 466 | slope, F, PCTL, PVF, PVF_RF = [0]*levels, [0]*levels, [0]*levels, [0]*levels, [0]*levels
|
|---|
| 467 | VF, VF_RF, MRVBF, MRRTF = [0]*levels, [0]*levels, [0]*levels, [0]*levels
|
|---|
| 468 |
|
|---|
| 469 | ###########################################################################
|
|---|
| 470 | # Step 1 (L=0)
|
|---|
| 471 | # Base scale resolution
|
|---|
| 472 | L = 0
|
|---|
| 473 | Xres_step.append(current_region.ewres)
|
|---|
| 474 | Yres_step.append(current_region.nsres)
|
|---|
| 475 | DEM.append(r_elevation)
|
|---|
| 476 | radi = 3
|
|---|
| 477 |
|
|---|
| 478 | g.message(os.linesep)
|
|---|
| 479 | g.message("Step {L}".format(L=L+1))
|
|---|
| 480 | g.message("------")
|
|---|
| 481 |
|
|---|
| 482 | # Calculation of slope (S1) and calculation of flatness (F1) (Equation 2)
|
|---|
| 483 | grass.message("Calculation of slope and transformation to flatness F{L}...".format(L=L+1))
|
|---|
| 484 | slope[L] = get_slope(L, DEM[L])
|
|---|
| 485 | F[L] = get_flatness(L, slope[L], t_slope, p_slope)
|
|---|
| 486 |
|
|---|
| 487 | # Calculation of elevation percentile PCTL for step 1
|
|---|
| 488 | grass.message("Calculation of elevation percentile PCTL{L}...".format(L=L+1))
|
|---|
| 489 | PCTL[L] = get_percentile(L, DEM[L], radi, moving_window_square)
|
|---|
| 490 |
|
|---|
| 491 | # Transform elevation percentile to local lowness for step 1 (Equation 3)
|
|---|
| 492 | grass.message("Calculation of preliminary valley flatness index PVF{L}...".format(L=L+1))
|
|---|
| 493 | PVF[L] = get_prelim_flatness(L, F[L], PCTL[L], t_pctl_v, p_pctl)
|
|---|
| 494 | if mrrtf != '':
|
|---|
| 495 | grass.message("Calculation of preliminary ridge top flatness index PRF{L}...".format(L=L+1))
|
|---|
| 496 | PVF_RF[L] = get_prelim_flatness_rf(L, F[L], PCTL[L], t_pctl_r, p_pctl)
|
|---|
| 497 |
|
|---|
| 498 | # Calculation of the valley flatness step 1 VF1 (Equation 4)
|
|---|
| 499 | grass.message("Calculation of valley flatness VF{L}...".format(L=L+1))
|
|---|
| 500 | VF[L] = get_valley_flatness(L, PVF[L], t_vf, p_slope)
|
|---|
| 501 | if mrrtf != '':
|
|---|
| 502 | grass.message("Calculation of ridge top flatness RF{L}...".format(L=L+1))
|
|---|
| 503 | VF_RF[L] = get_valley_flatness(L, PVF_RF[L], t_rf, p_slope)
|
|---|
| 504 |
|
|---|
| 505 | ##################################################################################
|
|---|
| 506 | # Step 2 (L=1)
|
|---|
| 507 | # Base scale resolution
|
|---|
| 508 | L = 1
|
|---|
| 509 | Xres_step.append(current_region.ewres)
|
|---|
| 510 | Yres_step.append(current_region.nsres)
|
|---|
| 511 | DEM.append(r_elevation)
|
|---|
| 512 | t_slope /= 2.0
|
|---|
| 513 | radi = 6
|
|---|
| 514 |
|
|---|
| 515 | grass.message(os.linesep)
|
|---|
| 516 | grass.message("Step {L}".format(L=L+1))
|
|---|
| 517 | grass.message("------")
|
|---|
| 518 |
|
|---|
| 519 | # Calculation of flatness for step 2 (Equation 5)
|
|---|
| 520 | # The second step commences the same way with the original DEM at its base resolution,
|
|---|
| 521 | # using a slope threshold ts,2 half of ts,1:
|
|---|
| 522 | grass.message("Calculation of flatness F{L}...".format(L=L+1))
|
|---|
| 523 | F[L] = get_flatness(L, slope[L-1], t_slope, p_slope)
|
|---|
| 524 |
|
|---|
| 525 | # Calculation of elevation percentile PCTL for step 2 (radius of 6 cells)
|
|---|
| 526 | grass.message("Calculation of elevation percentile PCTL{L}...".format(L=L+1))
|
|---|
| 527 | PCTL[L] = get_percentile(L, r_elevation, radi, moving_window_square)
|
|---|
| 528 |
|
|---|
| 529 | # PVF for step 2 (Equation 6)
|
|---|
| 530 | grass.message("Calculation of preliminary valley flatness index PVF{L}...".format(L=L+1))
|
|---|
| 531 | PVF[L] = get_prelim_flatness(L, F[L], PCTL[L], t_pctl_v, p_pctl)
|
|---|
| 532 | if mrrtf != '':
|
|---|
| 533 | grass.message("Calculation of preliminary ridge top flatness index PRF{L}...".format(L=L+1))
|
|---|
| 534 | PVF_RF[L] = get_prelim_flatness_rf(L, F[L], PCTL[L], t_pctl_r, p_pctl)
|
|---|
| 535 |
|
|---|
| 536 | # Calculation of the valley flatness VF for step 2 (Equation 7)
|
|---|
| 537 | grass.message("Calculation of valley flatness VF{L}...".format(L=L+1))
|
|---|
| 538 | VF[L] = get_valley_flatness(L, PVF[L], t_vf, p_slope)
|
|---|
| 539 | if mrrtf != '':
|
|---|
| 540 | grass.message("Calculation of ridge top flatness RF{L}...".format(L=L+1))
|
|---|
| 541 | VF_RF[L] = get_valley_flatness(L, PVF_RF[L], t_rf, p_slope)
|
|---|
| 542 |
|
|---|
| 543 | # Calculation of MRVBF for step 2
|
|---|
| 544 | grass.message("Calculation of MRVBF{L}...".format(L=L+1))
|
|---|
| 545 | MRVBF[L] = get_mrvbf(L, VF_Lminus1=VF[L-1], VF_L=VF[L], t=t_pctl_v)
|
|---|
| 546 | if mrrtf != '':
|
|---|
| 547 | grass.message("Calculation of MRRTF{L}...".format(L=L+1))
|
|---|
| 548 | MRRTF[L] = get_mrvbf(L, VF_Lminus1=VF_RF[L-1], VF_L=VF_RF[L], t=t_pctl_r)
|
|---|
| 549 |
|
|---|
| 550 | # Update flatness for step 2 with combined flatness from F1 and F2 (Equation 10)
|
|---|
| 551 | grass.message("Calculation of combined flatness index CF{L}...".format(L=L+1))
|
|---|
| 552 | F[L] = get_combined_flatness(L, F[L-1], F[L])
|
|---|
| 553 |
|
|---|
| 554 | ##################################################################################
|
|---|
| 555 | # Remaining steps
|
|---|
| 556 | # DEM_1_1 refers to scale (smoothing) and resolution (cell size)
|
|---|
| 557 | # so that DEM_L1_L-1 refers to smoothing of current step,
|
|---|
| 558 | # but resolution of previous step
|
|---|
| 559 |
|
|---|
| 560 | for L in range(2, levels):
|
|---|
| 561 |
|
|---|
| 562 | t_slope /= 2.0
|
|---|
| 563 | Xres_step.append(Xres_step[L-1] * 3)
|
|---|
| 564 | Yres_step.append(Yres_step[L-1] * 3)
|
|---|
| 565 | radi = 6
|
|---|
| 566 |
|
|---|
| 567 | # delete temporary maps from L-2
|
|---|
| 568 | for tmap in TMP_RAST[L-2]:
|
|---|
| 569 | g.remove(type='raster', name=tmap, flags='f', quiet=True)
|
|---|
| 570 |
|
|---|
| 571 | grass.message(os.linesep)
|
|---|
| 572 | grass.message("Step {L}".format(L=L+1))
|
|---|
| 573 | grass.message("------")
|
|---|
| 574 |
|
|---|
| 575 | # Coarsen resolution to resolution of prevous step (step L-1) and smooth DEM
|
|---|
| 576 | if L >= 3:
|
|---|
| 577 | grass.run_command('g.region', ewres = Xres_step[L-1], nsres = Yres_step[L-1])
|
|---|
| 578 | grass.message('Coarsening resolution to ew_res={e} and ns_res={n}...'.format(
|
|---|
| 579 | e=Xres_step[L-1], n=Yres_step[L-1]))
|
|---|
| 580 |
|
|---|
| 581 | grass.message("DEM smoothing 11 x 11 windows with Gaussian smoothing kernel (sigma) 3...")
|
|---|
| 582 | DEM.append(get_smoothed_dem(L, DEM[L-1]))
|
|---|
| 583 |
|
|---|
| 584 | # Calculate slope
|
|---|
| 585 | grass.message("Calculation of slope...")
|
|---|
| 586 | slope[L] = get_slope(L, DEM[L])
|
|---|
| 587 |
|
|---|
| 588 | # Refine slope to base resolution
|
|---|
| 589 | if L >= 3:
|
|---|
| 590 | grass.message('Resampling slope back to base resolution...')
|
|---|
| 591 | slope[L] = refine(L, slope[L], current_region, method='bilinear')
|
|---|
| 592 |
|
|---|
| 593 | # Coarsen resolution to current step L and calculate PCTL
|
|---|
| 594 | grass.run_command('g.region', ewres=Xres_step[L], nsres=Yres_step[L])
|
|---|
| 595 | DEM[L] = refine(L, DEM[L], Region(), method = 'average')
|
|---|
| 596 | grass.message("Calculation of elevation percentile PCTL{L}...".format(L=L+1))
|
|---|
| 597 | PCTL[L] = get_percentile(L, DEM[L], radi, moving_window_square)
|
|---|
| 598 |
|
|---|
| 599 | # Refine PCTL to base resolution
|
|---|
| 600 | grass.message("Resampling PCTL{L} to base resolution...".format(L=L+1))
|
|---|
| 601 | PCTL[L] = refine(L, PCTL[L], current_region, method='bilinear')
|
|---|
| 602 |
|
|---|
| 603 | # Calculate flatness F at the base resolution
|
|---|
| 604 | grass.message("Calculate F{L} at base resolution...".format(L=L+1))
|
|---|
| 605 | F[L] = get_flatness(L, slope[L], t_slope, p_slope)
|
|---|
| 606 |
|
|---|
| 607 | # Update flatness with combined flatness CF from the previous step
|
|---|
| 608 | grass.message("Calculate combined flatness CF{L} at base resolution...".format(L=L+1))
|
|---|
| 609 | F[L] = get_combined_flatness(L, F1=F[L-1], F2=F[L])
|
|---|
| 610 |
|
|---|
| 611 | # Calculate preliminary valley flatness index PVF at the base resolution
|
|---|
| 612 | grass.message("Calculate preliminary valley flatness index PVF{L} at base resolution...".format(L=L+1))
|
|---|
| 613 | PVF[L] = get_prelim_flatness(L, F[L], PCTL[L], t_pctl_v, p_pctl)
|
|---|
| 614 | if mrrtf != '':
|
|---|
| 615 | grass.message("Calculate preliminary ridge top flatness index PRF{L} at base resolution...".format(L=L+1))
|
|---|
| 616 | PVF_RF[L] = get_prelim_flatness_rf(L, F[L], PCTL[L], t_pctl_r, p_pctl)
|
|---|
| 617 |
|
|---|
| 618 | # Calculate valley flatness index VF
|
|---|
| 619 | grass.message("Calculate valley flatness index VF{L} at base resolution...".format(L=L+1))
|
|---|
| 620 | VF[L] = get_valley_flatness(L, PVF[L], t_vf, p_slope)
|
|---|
| 621 | if mrrtf != '':
|
|---|
| 622 | grass.message("Calculate ridge top flatness index RF{L} at base resolution...".format(L=L+1))
|
|---|
| 623 | VF_RF[L] = get_valley_flatness(L, PVF_RF[L], t_rf, p_slope)
|
|---|
| 624 |
|
|---|
| 625 | # Calculation of MRVBF
|
|---|
| 626 | grass.message("Calculation of MRVBF{L}...".format(L=L+1))
|
|---|
| 627 | MRVBF[L] = get_mrvbf(L, VF_Lminus1=MRVBF[L-1], VF_L=VF[L], t=t_pctl_v)
|
|---|
| 628 | if mrrtf != '':
|
|---|
| 629 | grass.message("Calculation of MRRTF{L}...".format(L=L+1))
|
|---|
| 630 | MRRTF[L] = get_mrvbf(L, VF_Lminus1=MRRTF[L-1], VF_L=VF_RF[L], t=t_pctl_r)
|
|---|
| 631 |
|
|---|
| 632 | # Output final MRVBF
|
|---|
| 633 | grass.mapcalc("$x = $y", x = mrvbf, y=MRVBF[L])
|
|---|
| 634 |
|
|---|
| 635 | if mrrtf != '':
|
|---|
| 636 | grass.mapcalc("$x = $y", x = mrrtf, y=MRRTF[L])
|
|---|
| 637 |
|
|---|
| 638 | if __name__ == "__main__":
|
|---|
| 639 | options, flags = grass.parser()
|
|---|
| 640 | atexit.register(cleanup)
|
|---|
| 641 | sys.exit(main())
|
|---|
