Get direction raster in clockwise degrees starting from the North

[landam]

Reported by cgravelm on 25 Mar 2015 22:54 UTC 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.

Operating system

MacOSX

GRASS GIS version and provenance

7.0.0

Migrated-From: https://trac.osgeo.org/grass/ticket/2637

[landam]

Attachment from hellik on 16 Feb 2018 17:27 UTC

https://trac.osgeo.org/grass/attachment/ticket/2637/aspect_nflag.png

[landam]

Attachment from hellik on 16 Feb 2018 17:27 UTC

https://trac.osgeo.org/grass/attachment/ticket/2637/aspect_nflag_color_color_manual_set.png

[landam]

Comment by annakrat on 26 Mar 2015 00:24 UTC 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.

[landam]

Comment by hellik on 26 Mar 2015 09:34 UTC Replying to [ticket:2637 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 http://trac.osgeo.org/grass/browser/grass-addons/grass7/raster/r.northerness.easterness/r.northerness.easterness.py#L61 aspect angles from cartesian (GRASS default) to compass angles for a r.mapcalc calculation

[landam]

Comment by neteler on 26 Mar 2015 09:48 UTC Replying to [comment:1 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.

[landam]

Comment by hellik on 26 Mar 2015 09:57 UTC Replying to [comment:3 neteler]:

Replying to [comment:1 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)

[landam]

Comment by cmbarton on 26 Mar 2015 19:47 UTC Replying to [comment:3 neteler]:

Replying to [comment:1 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

[landam]

Comment by glynn on 27 Mar 2015 16:03 UTC Replying to [comment:1 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.

[landam]

Comment by cmbarton on 27 Mar 2015 17:18 UTC 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.

[landam]

Comment by cmbarton on 27 Mar 2015 17:46 UTC Replying to [comment:7 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.

[landam]

Comment by neteler on 30 Jul 2015 21:19 UTC Ticket retargeted after 7.0.1 milestone closed

[landam]

Comment by neteler on 20 Nov 2015 17:08 UTC Ticket retargeted after milestone closed

[landam]

Comment by neteler on 28 Jan 2016 08:02 UTC Ticket retargeted after milestone closed

[landam]

Comment by neteler on 28 Jan 2016 08:06 UTC Ticket retargeted after 7.0.3 milestone closed

[landam]

Modified by @landam on 12 May 2016 06:36 UTC

[landam]

Modified by @landam on 25 Aug 2016 15:51 UTC

[landam]

Comment by @landam on 27 Aug 2016 13:42 UTC Milestone renamed

[landam]

Comment by neteler on 26 Jan 2018 11:40 UTC Ticket retargeted after milestone closed

[landam]

Comment by @landam on 28 Jan 2018 11:05 UTC What is status of this issue?

[landam]

Comment by hellik on 28 Jan 2018 13:11 UTC Replying to [comment:17 martinl]:

What is status of this issue?

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

[landam]

Comment by neteler on 28 Jan 2018 14:28 UTC 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.

[landam]

Comment by hellik on 28 Jan 2018 15:26 UTC Replying to [comment:19 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.

[landam]

Comment by cmbarton on 29 Jan 2018 06:27 UTC 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.

[landam]

Comment by mmetz on 29 Jan 2018 09:22 UTC Replying to [comment:21 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.

[landam]

Comment by hellik on 29 Jan 2018 09:38 UTC Replying to [comment:22 mmetz]:

Replying to [comment:21 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.

[landam]

Comment by hellik on 29 Jan 2018 09:51 UTC Replying to [comment:22 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?

[landam]

Comment by hellik on 29 Jan 2018 09:54 UTC Replying to [comment:24 hellik]:

Replying to [comment:22 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.

[landam]

Comment by mmetz on 29 Jan 2018 11:43 UTC Replying to [comment:24 hellik]:

Replying to [comment:22 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.

[landam]

Comment by hellik on 29 Jan 2018 12:01 UTC Replying to [comment:26 mmetz]:

Replying to [comment:24 hellik]:

Replying to [comment:22 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

[landam]

Comment by mmetz on 29 Jan 2018 14:02 UTC Replying to [comment:27 hellik]:

Replying to [comment:26 mmetz]:

Replying to [comment:24 hellik]:

Replying to [comment:22 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

[landam]

Comment by cmbarton on 29 Jan 2018 16:16 UTC This makes sense, especially as azimuth is east of north. It would be a welcome addition.

[landam]

Comment by mmetz on 11 Feb 2018 21:34 UTC Replying to [comment:29 cmbarton]:

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

Added to trunk https://trac.osgeo.org/grass/changeset/72220. Aspect as degrees CW from North with flat=-9999 can be produced with the new -n flag.

[landam]

Comment by cmbarton on 12 Feb 2018 01:24 UTC Thanks! Great news.

[landam]

Comment by neteler on 16 Feb 2018 13:44 UTC Replying to [comment:31 cmbarton]:

Thanks! Great news.

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

[landam]

Comment by hellik on 16 Feb 2018 14:04 UTC Replying to [comment:32 neteler]:

Replying to [comment:31 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.

[landam]

Comment by mmetz on 16 Feb 2018 15:50 UTC Replying to [comment:33 hellik]:

Replying to [comment:32 neteler]:

Replying to [comment:31 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.

[landam]

Comment by hellik on 16 Feb 2018 17:26 UTC Replying to [comment:34 mmetz]:

Replying to [comment:33 hellik]:

Replying to [comment:32 neteler]:

Replying to [comment:31 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

[landam]

Comment by hellik on 16 Feb 2018 22:08 UTC Replying to [comment:35 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?

[landam]

Comment by hellik on 16 Feb 2018 22:16 UTC Replying to [comment:36 hellik]:

Replying to [comment:35 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 (?)....

[landam]

Comment by mmetz on 18 Feb 2018 16:32 UTC Replying to [comment:37 hellik]:

Replying to [comment:36 hellik]:

Replying to [comment:35 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.

[landam]

Modified by neteler on 12 Jun 2018 20:48 UTC

[landam]

Comment by @landam on 25 Sep 2018 16:52 UTC All enhancement tickets should be assigned to 7.6 milestone.

[landam]

Comment by @landam on 25 Jan 2019 21:07 UTC Ticket retargeted after milestone closed