How can I view the equations used in the GRASS GIS modules? I found the code here: 1) URL: https://trac.osgeo.org/…/re…/raster/r.watershed/front/main.c; 2) URL: https://github.com/OSGeo/grass/tree/master/raster/r.watershed?fbclid=IwAR24KBPYl8iQjGQMphtl54dri5T0hU8cU5Q9DYycgAjeFG7RYQ8SVvAMn0Q But not the expressions.
-
1Dig deeper? See github.com/OSGeo/grass/blob/master/raster/r.watershed/seg/… What exact equation are you looking for? – mankoff Jan 10 '20 at 19:25
-
1Hint: Out grassbook.org comes with an appendix, explaining in-depth numerous equations implemented in GRASS GIS. – markusN Jan 10 '20 at 20:40
-
@mankoff, for LS equations. a) NETELER & MITASOVA (2008, pp. 158-166, 382), display LS equation in SI units; b) RENARD et al. (1997, p. 104) display LS equation not SI units; c) MITASOVA et al. (1996), display LS equation in SI units. d) github.com/OSGeo/grass/blob/master/raster/r.watershed/seg/…, display LS equation not in SI units? I see some terms "feet". Note that d)(code) is different from expressions in a). Was expecting form something like: lsfac=1.6*exp(flowacc*resolution/22.1,0.6)*exp(sin(slope)/0.09,1.3). Any comments? The LS GRASS code uses SI units? – Valter Albino Jan 10 '20 at 23:37
Add a comment
|
Concerning to the code (https://github.com/OSGeo/grass/blob/master/raster/r.watershed/seg/sg_factor.c#L68), the GUI gets the MDE in SI units and converts them to USA feet units:
int sg_factor(void)
{
int r, c;
CELL low_elev, hih_elev;
double height, length, S, sin_theta;
WAT_ALT wa;
ASP_FLAG af;
G_message(_("SECTION 5: RUSLE LS and/or S factor determination."));
for (r = nrows - 1; r >= 0; r--) {
G_percent(nrows - r, nrows, 3);
for (c = ncols - 1; c >= 0; c--) {
seg_get(&aspflag, (char *)&af, r, c);
if (FLAG_GET(af.flag, NULLFLAG))
continue;
seg_get(&watalt, (char *)&wa, r, c);
low_elev = wa.ele;
cseg_get(&r_h, &hih_elev, r, c);
dseg_get(&s_l, &length, r, c);
height = 1.0 * (hih_elev - low_elev) / ele_scale;
if (length > max_length) {
height *= max_length / length;
length = max_length;
}
sin_theta = height / sqrt(height * height + length * length);
if (height / length < .09)
S = 10.8 * sin_theta + .03;
else
S = 16.8 * sin_theta - .50;
if (ls_flag) {
length *= METER_TO_FOOT;
len_slp_equ(length, sin_theta, S, r, c);
}
if (sg_flag) {
dseg_put(&s_g, &S, r, c);
}
}
}
G_percent(nrows, nrows, 1); /* finish it */
return 0;
}
int len_slp_equ(double slope_length, double sin_theta, double S, int r, int c)
{
double rill, s_l_exp, /* m */
rill_ratio, /* Beta */
LS;
rill_ratio = (sin_theta / 0.0896) / (3.0 * pow(sin_theta, 0.8) + 0.56);
if (ril_flag) {
dseg_get(&ril, &rill, r, c);
}
else if (ril_value >= 0.0) {
rill = ril_value;
}
else
rill = 0.0;
/* rill_ratio equation from Steve Warren */
rill_ratio *= .5 + .005 * rill + .0001 * rill * rill;
s_l_exp = rill_ratio / (1 + rill_ratio);
LS = S * pow((slope_length / 72.6), s_l_exp);
dseg_put(&l_s, &LS, r, c);
return 0;
}
