Get direction raster in clockwise degrees starting from the North

[cgravelm]

I have been looking for a way to create a raster direction map using a DEM in GRASS (I am currently using the stable 7.0.0 version). So far, the 3 different scripts that do so (r.param.scale, r.shaded.aspect, and r.fill.dir) assume that 0 is at the East and count counterclockwise from there. The only way I can create a clockwise direction map with 0 at the North is to use the format "agnps" in r.fill.dir. However, this creates a map with directions ranging from 0 (equivalent to 0) to 8 (equivalent to 360 degrees), which can be easily transformed into a 0-360 scale, but lacks precision.

Would it be possible to add a flag in one of those scripts to create maps in degrees that have 0 to the North, 90 to the East and so forth? It would make it easier to integrate GRASS rasters to agent-based models.

[annakrat]

You can just use r.mapcalc "angle_cw = -angle_ccw + 450", this will give you what you need. It outputs angles from 90 to 450, but that's typically not a problem.

I guess it could be implemented in r.slope.aspect.

[hellik]

Replying to cgravelm:

I have been looking for a way to create a raster direction map using a DEM in GRASS (I am currently using the stable 7.0.0 version). So far, the 3 different scripts that do so (r.param.scale, r.shaded.aspect, and r.fill.dir) assume that 0 is at the East and count counterclockwise from there. The only way I can create a clockwise direction map with 0 at the North is to use the format "agnps" in r.fill.dir. However, this creates a map with directions ranging from 0 (equivalent to 0) to 8 (equivalent to 360 degrees), which can be easily transformed into a 0-360 scale, but lacks precision.

Would it be possible to add a flag in one of those scripts to create maps in degrees that have 0 to the North, 90 to the East and so forth? It would make it easier to integrate GRASS rasters to agent-based models.

have a look e.g. at ​aspect angles from cartesian (GRASS default) to compass angles for a r.mapcalc calculation

[neteler]

Replying to annakrat:

You can just use r.mapcalc "angle_cw = -angle_ccw + 450", this will give you what you need. It outputs angles from 90 to 450, but that's typically not a problem.

... or you solve it with an if() condition.

Shell script solution:

rotate_angle()
{
 is_negative=`echo "$1" | awk '{printf "%d\n", $1 < 0. ? 1 : 0}'`

 if [ $is_negative -eq 0 ] ; then
   tmp=`echo "$1"   | awk '{printf "%f\n", 360. - $1 + 90.}'`
   tmp=`echo "$tmp" | awk '{printf "%f\n", $1 >= 360. ? $1 - 360. : $1}'`
   echo "$tmp"
 else
   echo "NA"
 fi
}

rotate_angle 270
#[1] 180
rotate_angle 180
#[1] 270

Likewise you could use r.mapcalc and its eval() function.

I guess it could be implemented in r.slope.aspect.

Note the (old) patch:

raster/r.slope.aspect/r_sl_asp_northangle_diffs.tar.gz

Important: the output of r.slope.aspect can be used as input elsewhere, hence a flag would increase the risk to mess up subsequent calculations.

[hellik]

Replying to neteler:

Replying to annakrat:

You can just use r.mapcalc "angle_cw = -angle_ccw + 450", this will give you what you need. It outputs angles from 90 to 450, but that's typically not a problem.

... or you solve it with an if() condition.

taken from a python script

grass.mapcalc("$outmap = if( $cartesian == 0, 0, if( $cartesian < 90, 90 - $cartesian, 360 + 90 - $cartesian) )",
    outmap = r_aspect_compass, 
    cartesian = r_aspect)

[cmbarton]

Replying to neteler:

Replying to annakrat:

You can just use r.mapcalc "angle_cw = -angle_ccw + 450", this will give you what you need. It outputs angles from 90 to 450, but that's typically not a problem.

... or you solve it with an if() condition.

Shell script solution:

rotate_angle()
{
 is_negative=`echo "$1" | awk '{printf "%d\n", $1 < 0. ? 1 : 0}'`

 if [ $is_negative -eq 0 ] ; then
   tmp=`echo "$1"   | awk '{printf "%f\n", 360. - $1 + 90.}'`
   tmp=`echo "$tmp" | awk '{printf "%f\n", $1 >= 360. ? $1 - 360. : $1}'`
   echo "$tmp"
 else
   echo "NA"
 fi
}

rotate_angle 270
#[1] 180
rotate_angle 180
#[1] 270

Likewise you could use r.mapcalc and its eval() function.

I guess it could be implemented in r.slope.aspect.

Note the (old) patch:

raster/r.slope.aspect/r_sl_asp_northangle_diffs.tar.gz

Important: the output of r.slope.aspect can be used as input elsewhere, hence a flag would increase the risk to mess up subsequent calculations.

I don't see how a flag would mess up other calculations. Different other routines use different ways of representing aspect. So you need to know what is needed and select that even now. There are also potential uses that could benefit by having aspect in the cardinal directions. It is a pain to always have to convert the output to make it readable for humans in normal ways too. A simple flag seems like a nice improvement here.

Michael

[glynn]

Replying to annakrat:

You can just use r.mapcalc "angle_cw = -angle_ccw + 450", this will give you what you need. It outputs angles from 90 to 450, but that's typically not a problem.

Or, if you want to limit the range to 0..360, use:

r.mapcalc "angle_cw = (450 - angle_ccw) % 360"

Similarly, if you want -180..180 (again, clockwise from North):

r.mapcalc "angle_cw = (630 - angle_ccw) % 360 - 180"

Similar formulae can be used for any other convention. The main caveat is that the left-hand operand to the % (modulo) operator needs to be non-negative in order to obtain a non-negative result.

[cmbarton]

Thanks for the simple mapcalc approach. I think, however, the point is not that this cannot be calculated in mapcalc. It is that it would be convenient to have the primary aspect module for GRASS have an option to calculate aspect in cardinal degrees from north, what most users would expect from an aspect calculation in a GIS. The counter clockwise from E. is a path dependent legacy from the early days of GIS when rasters were treated like a 2D graph. Other modules have come to expect that, which is OK. Perhaps it even makes some math easier in some uses. But we no longer calculate raster cell position from the lower left. So it is odd that r.slope.aspect does not at least have an option to treat aspect like a geographic value instead of only like a vector angle on a graph.

[cmbarton]

Replying to cmbarton:

Thanks for the simple mapcalc approach. I think, however, the point is not that this cannot be calculated in mapcalc. It is that it would be convenient to have the primary aspect module for GRASS have an option to calculate aspect in cardinal degrees from north, what most users would expect from an aspect calculation in a GIS. The counter clockwise from E. is a path dependent legacy from the early days of GIS when rasters were treated like a 2D graph. Other modules have come to expect that, which is OK. Perhaps it even makes some math easier in some uses. But we no longer calculate raster cell position from the lower left. So it is odd that r.slope.aspect does not at least have an option to treat aspect like a geographic value instead of only like a vector angle on a graph.

The other thing is that all the mapcalc solutions require 2 passes through a map to get aspect as degrees from north. If this were a flag in r.slope.aspect, it could be done in 1 pass. This is important if you are doing a lot of big maps.

[neteler]

Ticket retargeted after 7.0.1 milestone closed

[neteler]

Ticket retargeted after milestone closed

[neteler]

Ticket retargeted after milestone closed

[neteler]

Ticket retargeted after 7.0.3 milestone closed

[martinl]

Milestone renamed

[neteler]

Ticket retargeted after milestone closed

[martinl]

What is status of this issue?

[hellik]

Replying to martinl:

What is status of this issue?

a flag to calculate direction raster in clockwise degrees starting from the North isn't implemented yet

[neteler]

I am not sure at all that a flag would be useful, esp. since other modules depend on the current implementation.

See also comment:6 above on how to calculate it easily with r.mapcalc.

[hellik]

Replying to neteler:

I am not sure at all that a flag would be useful, esp. since other modules depend on the current implementation.

See also comment:6 above on how to calculate it easily with r.mapcalc.

I can't see that a flag interferes with other modules when the actual behavior is still default.

If not already there, adding the r.mapcalc examples may help to clarify for other users.

[cmbarton]

Such a flag would still be useful for a variety of applications and should not interfere with other modules if the default behavior remains the same. Note that in r.slope.aspect there is the choice of outputting slope as degrees or percent, with degrees being the default.

Just because outputting direction as in a compass bearing is not useful for some people does not mean that it would not be useful for other. For anyone interested in knowing what compass direction a location is facing (and that is a question that people have), the current output is completely misleading.

The manual states "The aspect output raster map indicates the direction that slopes are facing". To most people, the word "direction" refers to a compass direction = azimuth. If they keep reading, it gets increasingly baffling from this point of view. "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."

I understand the history of this (from before GRASS rasters were georegistered) and realize that because of this long legacy effect, other raster modules have been written to expect this kind of non-cartographic direction metric. But a simple flag to let users change this to represent what a map direction is commonly expected to mean seems like a useful enhancement.

[mmetz]

Replying to cmbarton:

Such a flag would still be useful for a variety of applications and should not interfere with other modules if the default behavior remains the same.

Such a flag would be easy to implement, but a clear definition of aspect with regard to zero or small slope is currently missing. Generally, for zero slope, zero aspect is produced, i.e. a value of zero does not mean East, but zero slope, while East is encoded as 360 degrees. However, if min_slope is > 0 and slope < min_slope, slope is set to zero, but aspect is set to NULL. This is not consistent. Therefore I suggest to produce zero aspect for zero slope or if slope is < min_slope. NULL is then reserved for those cells where slope and aspect can not be computed. This also means that when converting CCW from East to CW from North, North must be encoded as 360 degrees, not 0 degrees.

[hellik]

Replying to mmetz:

Replying to cmbarton:

Such a flag would still be useful for a variety of applications and should not interfere with other modules if the default behavior remains the same.

Such a flag would be easy to implement, but a clear definition of aspect with regard to zero or small slope is currently missing. Generally, for zero slope, zero aspect is produced, i.e. a value of zero does not mean East, but zero slope, while East is encoded as 360 degrees. However, if min_slope is > 0 and slope < min_slope, slope is set to zero, but aspect is set to NULL. This is not consistent. Therefore I suggest to produce zero aspect for zero slope or if slope is < min_slope. NULL is then reserved for those cells where slope and aspect can not be computed. This also means that when converting CCW from East to CW from North, North must be encoded as 360 degrees, not 0 degrees.

sounds reasonable.

[hellik]

Replying to mmetz:

Therefore I suggest to produce zero aspect for zero slope or if slope is < min_slope. NULL is then reserved for those cells where slope and aspect can not be computed. This also means that when converting CCW from East to CW from North, North must be encoded as 360 degrees, not 0 degrees.

what about to assign -1 for cells with zero slope or slope is < min_slope?

[hellik]

Replying to hellik:

Replying to mmetz:

Therefore I suggest to produce zero aspect for zero slope or if slope is < min_slope. NULL is then reserved for those cells where slope and aspect can not be computed. This also means that when converting CCW from East to CW from North, North must be encoded as 360 degrees, not 0 degrees.

what about to assign -1 for cells with zero slope or slope is < min_slope?

then North can start with 0 degrees; which is the most common use case.

[mmetz]

Replying to hellik:

Replying to mmetz:

Therefore I suggest to produce zero aspect for zero slope or if slope is < min_slope. NULL is then reserved for those cells where slope and aspect can not be computed. This also means that when converting CCW from East to CW from North, North must be encoded as 360 degrees, not 0 degrees.

what about to assign -1 for cells with zero slope or slope is < min_slope?

This would change the current default behaviour with min_slope = 0 and aspect = 0 if slope = 0. Direction = 0 has an equivalent meaning in r.cost, r.walk, r.drain, r.path, r.watershed, r.stream.extract: it is a valid cell, there is no way away from here, stop routing here.

Also, -1 can be interpreted as 359, just like 0 can be interpreted as 360, no improvement there. IMHO, the corresponding aspect value for zero slope must be clearly documented, then it can be easily changed if need be.

[hellik]

Replying to mmetz:

Replying to hellik:

Replying to mmetz:

Therefore I suggest to produce zero aspect for zero slope or if slope is < min_slope. NULL is then reserved for those cells where slope and aspect can not be computed. This also means that when converting CCW from East to CW from North, North must be encoded as 360 degrees, not 0 degrees.

what about to assign -1 for cells with zero slope or slope is < min_slope?

This would change the current default behaviour with min_slope = 0 and aspect = 0 if slope = 0. Direction = 0 has an equivalent meaning in r.cost, r.walk, r.drain, r.path, r.watershed, r.stream.extract: it is a valid cell, there is no way away from here, stop routing here.

Also, -1 can be interpreted as 359, just like 0 can be interpreted as 360, no improvement there. IMHO, the corresponding aspect value for zero slope must be clearly documented, then it can be easily changed if need be.

a short note how other GIS calculate aspect:

​https://gis.stackexchange.com/questions/238999/understanding-aspect-units-in-qgis

[mmetz]

Replying to hellik:

Replying to mmetz:

Replying to hellik:

Replying to mmetz:

Therefore I suggest to produce zero aspect for zero slope or if slope is < min_slope. NULL is then reserved for those cells where slope and aspect can not be computed. This also means that when converting CCW from East to CW from North, North must be encoded as 360 degrees, not 0 degrees.

what about to assign -1 for cells with zero slope or slope is < min_slope?

This would change the current default behaviour with min_slope = 0 and aspect = 0 if slope = 0. Direction = 0 has an equivalent meaning in r.cost, r.walk, r.drain, r.path, r.watershed, r.stream.extract: it is a valid cell, there is no way away from here, stop routing here.

Also, -1 can be interpreted as 359, just like 0 can be interpreted as 360, no improvement there. IMHO, the corresponding aspect value for zero slope must be clearly documented, then it can be easily changed if need be.

a short note how other GIS calculate aspect:

​https://gis.stackexchange.com/questions/238999/understanding-aspect-units-in-qgis

no mentioning of aspect for flat areas.

gdaldem uses -9999 or 0 with -zero_for_flat

[cmbarton]

This makes sense, especially as azimuth is east of north. It would be a welcome addition.

[mmetz]

Replying to cmbarton:

This makes sense, especially as azimuth is east of north. It would be a welcome addition.

Added to trunk r72220. Aspect as degrees CW from North with flat=-9999 can be produced with the new -n flag.

[cmbarton]

Thanks! Great news.

[neteler]

Replying to cmbarton:

Thanks! Great news.

Please report if it works ok (It may be a backport candidate).

[hellik]

Replying to neteler:

Replying to cmbarton:

Thanks! Great news.

Please report if it works ok (It may be a backport candidate).

tested it a bit; seems to work so far.

the only thing for improving seems to assign a working raster color scheme for this new aspect type map.

as far as what I've seen, the result raster map is black.

[mmetz]

Replying to hellik:

Replying to neteler:

Replying to cmbarton:

Thanks! Great news.

Please report if it works ok (It may be a backport candidate).

tested it a bit; seems to work so far.

the only thing for improving seems to assign a working raster color scheme for this new aspect type map.

as far as what I've seen, the result raster map is black.

The current aspect color scheme is

0% black
50% white
100% black

assuming that aspect is in degrees (CCW from East or CW from North), a color scheme like

0 black
180 white
360 black

should work.

[hellik]

Replying to mmetz:

Replying to hellik:

Replying to neteler:

Replying to cmbarton:

Thanks! Great news.

Please report if it works ok (It may be a backport candidate).

tested it a bit; seems to work so far.

the only thing for improving seems to assign a working raster color scheme for this new aspect type map.

as far as what I've seen, the result raster map is black.

I mean the result map where the color scheme is automatically set by r.slope.aspect (see attached: aspect_nflag.png

The current aspect color scheme is

0% black
50% white
100% black

assuming that aspect is in degrees (CCW from East or CW from North), a color scheme like

0 black
180 white
360 black

should work.

this scheme looks good, see attached aspect_nflag_color_color_manual_set.png

[hellik]

Replying to hellik:

0 black
180 white
360 black

should work.

this scheme looks good, see attached aspect_nflag_color_color_manual_set.png

looking at aspect_nflag_color_color_manual_set.png, it seems to be just the other way around?

[hellik]

Replying to hellik:

Replying to hellik:

0 black
180 white
360 black

should work.

this scheme looks good, see attached aspect_nflag_color_color_manual_set.png

looking at aspect_nflag_color_color_manual_set.png, it seems to be just the other way around?

it seems to be more a visual effect as the human eye is used to get an illumination in maps of top-left (?)....

[mmetz]

Replying to hellik:

Replying to hellik:

Replying to hellik:

0 black
180 white
360 black

should work.

this scheme looks good, see attached aspect_nflag_color_color_manual_set.png

looking at aspect_nflag_color_color_manual_set.png, it seems to be just the other way around?

try

0 white
180 black
360 white

IMHO it's a matter of taste if you make north white and south black or vice versa.

it seems to be more a visual effect as the human eye is used to get an illumination in maps of top-left (?)....

Such an illumination applies to shaded reliefs, see e.g. the azimuth option of r.relief.

[martinl]

All enhancement tickets should be assigned to 7.6 milestone.

[martinl]

Ticket retargeted after milestone closed

[martinl]

Ticket retargeted after milestone closed