Converting atan2 Excel code to R, to produce flow direction in Cartesian Coordinates

Question

I have earth flow displacements in two directions, converted to point data, which I am trying to combine for an overall offset direction (point to where the flow is going). I have magnitude in a separate raster file.

I'm trying to convert an old excel code into R so I can work with larger datasets. I have an IF statement in excel which will convert 2 displacements into direction in Cartesian co-ordinates using atan2. I know that Excel and R have opposite argument placements for atan2, so I've swapped them.

Where it gets messy is after atan2, when I need to convert to degrees and "bin" the results in a correct quadrant. Ultimately, I'd like the values in Cartesian coordinates (0 to 360), rather than -180 to 180 (radians).

I've tried running through these already asked questions, but am still lost:

How to map atan2() to degrees 0-360
Reshaping EPA wind speed & direction data with dcast in R

Here is my excel code:

(N2 = the displacement in north direction)
(E2 = the displacement in east direction)

=IF(
    (ATAN2(N2,E2)*180/PI())<0, 360+(ATAN2(N2,E2)*180/PI()), 
    ATAN2(N2,E2)*180/PI()
)

I tried running this in R like this:

theta_rad <- atan2(E,N)
if(theta_rad < 0)
    theta_rad +2*pi

theta_deg <- theta_rad*(180/pi)`  # Convert to degrees

I've also tried the following, which just crashes R

dir <- (theta_rad > 0 ? x : (2*pi + theta_rad)) * 360 / (2*pi)

Warnings: The condition has length >1 and only the first element will be used.

Answer 1

Because atan2 returns on the domain of (-pi,pi] but you want it starting at 0 instead ((0,2pi] in radians), then modulus works for you here.

cbind(a=seq(-pi, pi, len=11), b=(2*pi+seq(-pi, pi, len=11))%%(2*pi))
#                a         b
#  [1,] -3.1415927 3.1415927
#  [2,] -2.5132741 3.7699112
#  [3,] -1.8849556 4.3982297
#  [4,] -1.2566371 5.0265482
#  [5,] -0.6283185 5.6548668
#  [6,]  0.0000000 0.0000000
#  [7,]  0.6283185 0.6283185
#  [8,]  1.2566371 1.2566371
#  [9,]  1.8849556 1.8849556
# [10,]  2.5132741 2.5132741
# [11,]  3.1415927 3.1415927

Converting this to degrees is just doing the standard conversion. I think your code should be

(360 + atan2(E, N) * 180 / pi) %% 360

Answer 2

thank you for your response. Based from your code, I was able to come up with an azimuth function as shown below. But I do not understand why I have to use either 450 or 90 on the add.degree variable:

add.degree <- ifelse(delta.xy[1] < 0 & delta.xy[2]>0, 450, 90)
return(add.degree-deg(atan2(delta.xy[2],delta.xy[1])))

Notheless, I found out later that the st_azimuth is now available as: nngeo::st_azimuth(st_point(.start),st_point(.end))

# quadrant 1
.start <- c(100.805, 79.600)
.end <- c(123.389, 123.400)
azimuth(.start, .end)
# [1] 27.27638

# quadrant 2
.start <- c(100.805, 123.400)
.end <- c(123.389, 79.600)
azimuth(.start, .end)
# [1] 152.7236

# quadrant 3
.start <- c(123.389, 123.400)
.end <- c(100.805, 79.600)
azimuth(.start, .end)
# [1] 207.2764

# quadrant 4
.start <- c(123.389, 79.600)
.end <- c(100.805, 123.400)
azimuth(.start, .end)
#[1] 332.7236

azimuth <- function(v.origin, v.ending) {
  delta.xy <- v.ending - v.origin
  add.degree <- ifelse(delta.xy[1] < 0 & delta.xy[2]>0, 450, 90)
  return(add.degree-deg(atan2(delta.xy[2],delta.xy[1])))
}