# Terrain aspect similarity

I have a point data set with terrain aspect at the location of the points, measured in degrees. I want to compute the similarity of aspect between the points, e.g. 70° is similar to 90°. This is straight forward, but the problem arises for aspects near the north direction. A aspect of 350° should be similar to 10°. How can I detect these two aspects to be similar?

I am looking for a transformation of the values that takes into account that 358° and 2° is almost the same, even though it is far apart in my 0-360 metric.

• You simply could check, whether the second aspect is either within a range or aspect 1 +/- your "similarity boundary" (e.g. 20°), or outside a range of +/- 339 degree? You do the second check only, if the first fails, so a nested if clause.
– Erik
Sep 24, 2020 at 12:15
• Ok, however, my aim is to use these attributes for a distance matrix (like the result of dist() in R) – so I have no direct control of the comparison. I need to transform my values before I calculate the distance matrix. Sep 24, 2020 at 12:17
• I guess we need some more background here, as well as the software you're working with.
– Erik
Sep 24, 2020 at 12:28
• My aim is to calculate a distance matrix of my locations based on several attributes, one of them is the aspect. I want to calculate the distance matrix with dist() in R. Points with an aspect near north (e.g. 355° and 5°) should have a close distance even though they are far apart on the scale. Sep 24, 2020 at 12:39

Compute the difference, and then take 360-the absolute difference if the absolute difference is greater than 180.

Here's one way of doing that in an R function:

``````d180 = function(a,b){
d = a-b
d = ifelse(abs(d)<180 , d , 360-abs(d))
d
}
``````

and some examples:

``````> d180(0,0)
[1] 0
> d180(0,1)
[1] -1
> d180(0,180)
[1] 180
> d180(0,181)
[1] 179
> d180(358,2)
[1] 4
> d180(2,358)
[1] 4
``````

Whether the distance is positive or negative is a bit odd, return the `abs(d)` if you don't care...

A bit of a different perspective, you could just linearize aspect. Here is pseudo-code pulled from the linear aspect tool in our ArcGIS Geomorphometry & Gradient Metrics toolbox. Note; we use `(mod (..) * 100)` and `(36000(360*100) / 100)` allow for two decimal points since `Fmod` is no longer available in ArcGIS.

``````a = aspect(elev) # degrees
t1 = focal(sin(a/57.296), window=3,"sum")
t2 = focal(cos(a/57.296), window=3,"sum")
linear.aspect = mod(((450-(atan2(t1, t2) * 57.296)) * 100), 36000) / 100
``````

Another nice trick is to look at the Stage (1978) slope-aspect transformation that provides an interaction. Assuming slope in degrees (Stage paper assumes percent) and aspect in degrees, you can calculate it thus:

``````s = (slope(elev) / 0.572957795786) * 0.01
a = aspect(elev) * (pi / 180)
scosa = s * cos(a)
scosa = s * sin(a)
``````

Example transformed values for 50% slope across 10 aspects.

``````Aspect    cosine    sine
N         0.500     0.000
N30E      0.433     0.250
N45E      0.345     0.345
N60E      0.250     0.433
E         0.000     0.500
ESE      -0.354     0.354
S        -0.500     0.000
SSW      -0.354    -0.354
W         0.000    -0.500
``````

Or, linearize with a recentered origin. This is done in vegetation modeling to provide a metric that is optimized for site prodictivity (eg. Roberts and Cooper 1989 TRASP transformation) Assuming degrees, `1 - cos( (pi / 180) (a - 30)) / 2`

``````a = aspect(elev)
trasp = (1 - cos( (3.142 / 180) * (a - 30)) ) / 2
``````

These transformations are available in ArcGIS in the Geomorphometry & Gradient Metrics toolbox and in the R package spatialEco (on CRAN).

Retaining sign (pseudo):

``````a = modulus((targetangle-sourceangle),6.283185307179586)
if (a>3.141592653589793)
return a-6.283185307179586
else
return a
``````

(Talking degrees:) If modulus(-45, 360) doesn't return 315 in your language, then use modulus((-45 + 360), 360). Your modulus fn must support decimal values.