[GRASS-dev] [GRASS GIS] #2637: Get direction raster in clockwise degrees starting from the North

#2637: Get direction raster in clockwise degrees starting from the North
-------------------------+--------------------------------------------------
Reporter: cgravelm | Owner: grass-dev@…
     Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.1
Component: Raster | Version: 7.0.0
Keywords: | Platform: MacOSX
      Cpu: Unspecified |
-------------------------+--------------------------------------------------
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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/2637&gt;
GRASS GIS <http://grass.osgeo.org>

#2637: Get direction raster in clockwise degrees starting from the North
-------------------------+--------------------------------------------------
Reporter: cgravelm | Owner: grass-dev@…
     Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.1
Component: Raster | Version: 7.0.0
Keywords: | Platform: MacOSX
      Cpu: Unspecified |
-------------------------+--------------------------------------------------

Comment(by 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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/2637#comment:1&gt;
GRASS GIS <http://grass.osgeo.org>

#2637: Get direction raster in clockwise degrees starting from the North
-------------------------+--------------------------------------------------
Reporter: cgravelm | Owner: grass-dev@…
     Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.1
Component: Raster | Version: 7.0.0
Keywords: | Platform: MacOSX
      Cpu: Unspecified |
-------------------------+--------------------------------------------------

Comment(by hellik):

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

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/2637#comment:2&gt;
GRASS GIS <http://grass.osgeo.org>

#2637: Get direction raster in clockwise degrees starting from the North
-------------------------+--------------------------------------------------
Reporter: cgravelm | Owner: grass-dev@…
     Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.1
Component: Raster | Version: 7.0.0
Keywords: | Platform: MacOSX
      Cpu: Unspecified |
-------------------------+--------------------------------------------------

Comment(by 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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/2637#comment:3&gt;
GRASS GIS <http://grass.osgeo.org>

#2637: Get direction raster in clockwise degrees starting from the North
-------------------------+--------------------------------------------------
Reporter: cgravelm | Owner: grass-dev@…
     Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.1
Component: Raster | Version: 7.0.0
Keywords: | Platform: MacOSX
      Cpu: Unspecified |
-------------------------+--------------------------------------------------

Comment(by hellik):

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)
}}}

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/2637#comment:4&gt;
GRASS GIS <http://grass.osgeo.org>

#2637: Get direction raster in clockwise degrees starting from the North
-------------------------+--------------------------------------------------
Reporter: cgravelm | Owner: grass-dev@…
     Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.1
Component: Raster | Version: 7.0.0
Keywords: | Platform: MacOSX
      Cpu: Unspecified |
-------------------------+--------------------------------------------------

Comment(by cmbarton):

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

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/2637#comment:5&gt;
GRASS GIS <http://grass.osgeo.org>

#2637: Get direction raster in clockwise degrees starting from the North
-------------------------+--------------------------------------------------
Reporter: cgravelm | Owner: grass-dev@…
     Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.1
Component: Raster | Version: 7.0.0
Keywords: | Platform: MacOSX
      Cpu: Unspecified |
-------------------------+--------------------------------------------------

Comment(by glynn):

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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/2637#comment:6&gt;
GRASS GIS <http://grass.osgeo.org>

#2637: Get direction raster in clockwise degrees starting from the North
-------------------------+--------------------------------------------------
Reporter: cgravelm | Owner: grass-dev@…
     Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.1
Component: Raster | Version: 7.0.0
Keywords: | Platform: MacOSX
      Cpu: Unspecified |
-------------------------+--------------------------------------------------

Comment(by 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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/2637#comment:7&gt;
GRASS GIS <http://grass.osgeo.org>

#2637: Get direction raster in clockwise degrees starting from the North
-------------------------+--------------------------------------------------
Reporter: cgravelm | Owner: grass-dev@…
     Type: enhancement | Status: new
Priority: normal | Milestone: 7.0.1
Component: Raster | Version: 7.0.0
Keywords: | Platform: MacOSX
      Cpu: Unspecified |
-------------------------+--------------------------------------------------

Comment(by cmbarton):

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.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/2637#comment:8&gt;
GRASS GIS <http://grass.osgeo.org>