[GRASS-SVN] r72229 - grass/trunk/raster/r.slope.aspect
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Feb 11 13:32:57 PST 2018
Author: mmetz
Date: 2018-02-11 13:32:57 -0800 (Sun, 11 Feb 2018)
New Revision: 72229
Modified:
grass/trunk/raster/r.slope.aspect/main.c
grass/trunk/raster/r.slope.aspect/r.slope.aspect.html
Log:
r.slope.aspect: add -n flag to produce aspect CW from North with flat = -9999 (like gdaldem)
Modified: grass/trunk/raster/r.slope.aspect/main.c
===================================================================
--- grass/trunk/raster/r.slope.aspect/main.c 2018-02-11 21:29:12 UTC (rev 72228)
+++ grass/trunk/raster/r.slope.aspect/main.c 2018-02-11 21:32:57 UTC (rev 72229)
@@ -57,6 +57,22 @@
* the sign bug for the second order derivatives (jh) */
+/* convert aspect from CCW from East to CW from North
+ * aspect for flat areas is set to -9999 */
+static double aspect_cw_n(double aspect)
+{
+ /* aspect == 0: flat */
+ if (aspect == 0)
+ return -9999;
+
+ /* no modulus because of fp values */
+ aspect = (450 - aspect);
+ if (aspect >= 360)
+ aspect -= 360;
+
+ return aspect;
+}
+
int main(int argc, char *argv[])
{
struct Categories cats;
@@ -142,7 +158,7 @@
} parm;
struct
{
- struct Flag *a;
+ struct Flag *a, *n;
} flag;
G_gisinit(argv[0]);
@@ -261,6 +277,13 @@
_("Do not align the current region to the raster elevation map");
flag.a->guisection = _("Settings");
+ flag.n = G_define_flag();
+ flag.n->key = 'n';
+ flag.n->label =
+ _("Create aspect as degrees clockwise from North (azimuth), with flat = -9999");
+ flag.n->description =
+ _("Default: degrees counter-clockwise from East, with flat = 0");
+ flag.n->guisection = _("Settings");
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -395,7 +418,7 @@
west = Rast_col_to_easting(0.5, &window);
V = G_distance(east, north, east, south) * 4 / (factor * zfactor);
H = G_distance(east, ns_med, west, ns_med) * 4 / (factor * zfactor);
- /* ____________________________
+ /* ____________________________
|c1 |c2 |c3 |
| | | |
| | north | |
@@ -792,7 +815,9 @@
} /* computing slope */
if (aspect_fd > 0) {
- if (key == 0.)
+ double aspect_flat = 0.;
+
+ if (slp_in_perc == 0.)
aspect = 0.;
else if (dx == 0) {
if (dy > 0)
@@ -802,28 +827,28 @@
}
else {
aspect = (atan2(dy, dx) / degrees_to_radians);
- if ((aspect <= 0.5) && (aspect > 0) &&
- out_type == CELL_TYPE)
- aspect = 360.;
if (aspect <= 0.)
aspect = 360. + aspect;
}
- /* if it's not the case that the slope for this cell
- is below specified minimum */
- if (!((slope_fd > 0) && (slp_in_perc < min_slope))) {
- if (out_type == CELL_TYPE)
- *((CELL *) asp_ptr) = (CELL) (aspect + .5);
- else
- Rast_set_d_value(asp_ptr,
- (DCELL) aspect, data_type);
+ if (flag.n->answer) {
+ aspect_flat = -9999;
+ aspect = aspect_cw_n(aspect);
}
+
+ if (out_type == CELL_TYPE) {
+ if (aspect > 0 && aspect < 0.5)
+ aspect = 360;
+ *((CELL *) asp_ptr) = (CELL) (aspect + .5);
+ }
else
- Rast_set_null_value(asp_ptr, 1, data_type);
+ Rast_set_d_value(asp_ptr,
+ (DCELL) aspect, data_type);
+
asp_ptr = G_incr_void_ptr(asp_ptr, Rast_cell_size(data_type));
/* now update min and max */
- if (min_asp > aspect)
+ if (aspect > aspect_flat && min_asp > aspect)
min_asp = aspect;
if (max_asp < aspect)
max_asp = aspect;
@@ -990,44 +1015,81 @@
we are using reverse order so that the label looked up
for i-.5 is not the one defined for i-.5, i+.5 interval, but
the one defile for i-1.5, i-.5 interval which is added later */
- for (i = ceil(max_asp); i >= 1; i--) {
- if (i == 360)
- sprintf(buf, "east");
- else if (i == 360)
- sprintf(buf, "east");
- else if (i == 45)
- sprintf(buf, "north ccw of east");
- else if (i == 90)
- sprintf(buf, "north");
- else if (i == 135)
- sprintf(buf, "north ccw of west");
- else if (i == 180)
- sprintf(buf, "west");
- else if (i == 225)
- sprintf(buf, "south ccw of west");
- else if (i == 270)
- sprintf(buf, "south");
- else if (i == 315)
- sprintf(buf, "south ccw of east");
- else
- sprintf(buf, "%d degree%s ccw from east", i,
- i == 1 ? "" : "s");
+ if (!flag.n->answer) {
+ for (i = ceil(max_asp); i >= 1; i--) {
+ if (i == 360)
+ sprintf(buf, "east");
+ else if (i == 45)
+ sprintf(buf, "north ccw of east");
+ else if (i == 90)
+ sprintf(buf, "north");
+ else if (i == 135)
+ sprintf(buf, "north ccw of west");
+ else if (i == 180)
+ sprintf(buf, "west");
+ else if (i == 225)
+ sprintf(buf, "south ccw of west");
+ else if (i == 270)
+ sprintf(buf, "south");
+ else if (i == 315)
+ sprintf(buf, "south ccw of east");
+ else
+ sprintf(buf, "%d degree%s ccw from east", i,
+ i == 1 ? "" : "s");
+ if (data_type == CELL_TYPE) {
+ Rast_set_c_cat((CELL *) &i, (CELL *) &i, buf, &cats);
+ continue;
+ }
+ tmp1 = (double)i - .5;
+ tmp2 = (double)i + .5;
+ Rast_set_d_cat(&tmp1, &tmp2, buf, &cats);
+ }
if (data_type == CELL_TYPE) {
- Rast_set_c_cat((CELL *) &i, (CELL *) &i, buf, &cats);
- continue;
+ cat = 0;
+ Rast_set_c_cat(&cat, &cat, "no aspect", &cats);
}
- tmp1 = (double)i - .5;
- tmp2 = (double)i + .5;
- Rast_set_d_cat(&tmp1, &tmp2, buf, &cats);
+ else {
+ tmp1 = 0;
+ Rast_set_d_cat(&tmp1, &tmp1, "no aspect", &cats);
+ }
}
- if (data_type == CELL_TYPE) {
- cat = 0;
- Rast_set_c_cat(&cat, &cat, "no aspect", &cats);
- }
else {
- tmp1 = 0.;
- tmp2 = .5;
- Rast_set_d_cat(&tmp1, &tmp2, "no aspect", &cats);
+ for (i = ceil(max_asp); i >= 1; i--) {
+ if (i == 0 && i == 360)
+ sprintf(buf, "north");
+ else if (i == 45)
+ sprintf(buf, "north-east");
+ else if (i == 90)
+ sprintf(buf, "east");
+ else if (i == 135)
+ sprintf(buf, "south-east");
+ else if (i == 180)
+ sprintf(buf, "south");
+ else if (i == 225)
+ sprintf(buf, "south-west");
+ else if (i == 270)
+ sprintf(buf, "west");
+ else if (i == 315)
+ sprintf(buf, "north-west");
+ else
+ sprintf(buf, "%d degree%s cw from north", i,
+ i == 1 ? "" : "s");
+ if (data_type == CELL_TYPE) {
+ Rast_set_c_cat((CELL *) &i, (CELL *) &i, buf, &cats);
+ continue;
+ }
+ tmp1 = (double)i - .5;
+ tmp2 = (double)i + .5;
+ Rast_set_d_cat(&tmp1, &tmp2, buf, &cats);
+ }
+ if (data_type == CELL_TYPE) {
+ cat = -9999;
+ Rast_set_c_cat(&cat, &cat, "no aspect", &cats);
+ }
+ else {
+ tmp1 = -9999;
+ Rast_set_d_cat(&tmp1, &tmp1, "no aspect", &cats);
+ }
}
Rast_write_cats(aspect_name, &cats);
Rast_free_cats(&cats);
Modified: grass/trunk/raster/r.slope.aspect/r.slope.aspect.html
===================================================================
--- grass/trunk/raster/r.slope.aspect/r.slope.aspect.html 2018-02-11 21:29:12 UTC (rev 72228)
+++ grass/trunk/raster/r.slope.aspect/r.slope.aspect.html 2018-02-11 21:32:57 UTC (rev 72229)
@@ -18,28 +18,35 @@
parameter <b>zscale</b> must not be used</em>.
<p>
-The <b>aspect</b> output raster map indicates the direction that slopes are
-facing. The aspect categories represent the number degrees of east. Category
-and color table files are also generated for the aspect raster map. The aspect
-categories represent the number degrees of east and they increase
-counterclockwise: 90 degrees is North, 180 is West, 270 is South 360 is East.<br>
-Note: These values can be transformed to azimuth (0 is North, 90 is East, etc)
-values using <a href="r.mapcalc.html">r.mapcalc</a>:
+The <b>aspect</b> output raster map indicates the direction that slopes
+are facing counterclockwise from East: 90 degrees is North, 180 is
+West, 270 is South, 360 is East. Zero aspect indicates flat areas with
+zero slope. Category and color table files are also generated for the
+aspect raster map. <br> Note: These values can be transformed to
+azimuth values (90 is East, 180 is South, 270 is West, 360 is North)
+using <a href="r.mapcalc.html">r.mapcalc</a>:
<div class="code"><pre>
-# convert angles from CCW to north up
-r.mapcalc "azimuth_aspect = (450 - ccw_aspect) % 360"
+# convert angles from CCW from East to CW from North
+# modulus (%) can not be used with floating point aspect values
+r.mapcalc "azimuth_aspect = if(ccw_aspect == 0, 0, \
+ if(ccw_aspect < 90, 90 - ccw_aspect, \
+ 450 - ccw_aspect)))"
</pre></div>
+Alternatively, the <b>-n</b> flag can be used to produce aspect as
+degrees CW from North. Aspect for flat areas is then set to -9999
+(default: 0).
+
<p>
-The aspect is not defined for slope equal to zero.
-Thus, most cells with a very small slope end up having category 0,
-45, ..., 360 in <b>aspect</b> output.
-It is possible to reduce the bias in these directions
-by filtering out the aspect in areas where the terrain is almost flat.
-A option <b>min_slope</b> can be used to specify the minimum slope
-for which aspect is computed. The aspect for all cells with
-slope < <b>min_slope</b> is set to <i>null</i> (no-data).
+The aspect for slope equal to zero (flat areas) is set to zero (-9999
+with <b>-n</b> flag). Thus, most cells with a very small slope end up
+having category 0, 45, ..., 360 in <b>aspect</b> output. It is possible
+to reduce the bias in these directions by filtering out the aspect in
+areas where the terrain is almost flat. A option <b>min_slope</b> can
+be used to specify the minimum slope for which aspect is computed. For
+all cells with slope < <b>min_slope</b>, both slope and
+aspect are set to zero.
<center>
<img src="aspect_diagram.png" border="0">
@@ -152,11 +159,11 @@
The current mask is ignored.
<p>
-The algorithm used to determine slope and aspect uses a 3x3 neighborhood
-around each cell in the raster elevation map. Thus, it is not possible to determine
-slope and aspect for the cells adjacent to the edges in the elevation map
-layer. These cells are assigned a "zero slope" value (category 0) in both
-the slope and aspect raster maps.
+The algorithm used to determine slope and aspect uses a 3x3
+neighborhood around each cell in the raster elevation map. Thus, it is
+not possible to determine slope and aspect for the cells adjacent to
+the edges and NULL cells in the elevation map layer. These cells are
+set to nodata in both the slope and aspect raster maps.
<p>
Horn's formula is used to find the first order derivatives in x and y directions.
@@ -206,8 +213,8 @@
<div class="code"><pre>
g.region raster=elevation -p
-# generate aspect map with CCW orientation
-r.slope.aspect elevation=elevation aspect=myaspect
+# generate integer aspect map with degrees CCW from East
+r.slope.aspect elevation=elevation aspect=myaspect precision=CELL
# generate compass orientation and classify four major directions (N, E, S, W)
r.mapcalc "aspect_4_directions = eval( \\
More information about the grass-commit
mailing list