Skip to content

add practical examples #7

@espinielli

Description

@espinielli

From a request via email from Michael Henry from Australia about great circle intersection of
[point, bearing]_1 with [point, bearing]_2:

library(nvctr)
#> Warning: package 'nvctr' was built under R version 4.0.3

# find point on great circle at distance, from Example 8
point_at_distance <- function(n_EA_E, distance, azimuth, r_Earth = 6371e3) {
  k_east_E <- unit(pracma::cross(base::t(R_Ee()) %*%
                                   c(1, 0, 0) %>%
                                   as.vector(), n_EA_E))
  k_north_E <- pracma::cross(n_EA_E, k_east_E)
  d_E <- k_north_E * cos(azimuth) + k_east_E * sin(azimuth)
  n_EB_E <- n_EA_E * cos(distance / r_Earth) + d_E * sin(distance / r_Earth)
  n_E2lat_lon(n_EB_E)
}

# test 1
A <- c(-10, 0)
n_EA_E <- lat_lon2n_E(rad(A[1]), rad(A[2]))
azimuth <- rad(0)
s_AB <- 100000 # distance (m)
n_EB_E <- point_at_distance(n_EA_E, s_AB, azimuth)
(B <-  n_EB_E %>% deg()) # reassuringly longitude remains 0
#> [1] -9.100678  0.000000

# test 2
A <- c(0, -10)
n_EA_E <- lat_lon2n_E(rad(A[1]), rad(A[2]))
azimuth <- rad(90)
s_AB <- 100000 # distance (m)
n_EB_E <- point_at_distance(n_EA_E, s_AB, azimuth)
(B <-  n_EB_E %>% deg()) # reassuringly latitude remains 0 (well, close enough, 5.5e-17 !)
#> [1]  5.506349e-17 -9.100678e+00


# great circle intersection, from Example 9
great_circle_intersection <- function(n_EA1_E, n_EA2_E, n_EB1_E, n_EB2_E) {
  n_EC_E_tmp <- unit(pracma::cross(
    pracma::cross(n_EA1_E, n_EA2_E),
    pracma::cross(n_EB1_E, n_EB2_E)))
  n_EC_E <- sign(pracma::dot(n_EC_E_tmp, n_EA1_E)) * n_EC_E_tmp
  
  n_EC_E
}

# (too?) simple test
A1 <- c(-10, 0)
n_EA1_E <- lat_lon2n_E(rad(A1[1]), rad(A1[2]))
A2 <- c(10, 0)
n_EA2_E <- lat_lon2n_E(rad(A2[1]), rad(A2[2]))
B1 <- c(0, -10)
n_EB1_E <- lat_lon2n_E(rad(B1[1]), rad(B1[2]))
B2 <- c(0, 10)
n_EB2_E <- lat_lon2n_E(rad(B2[1]), rad(B2[2]))

n_EC_E <- great_circle_intersection(n_EA1_E, n_EA2_E, n_EB1_E, n_EB2_E)
(C <- n_E2lat_lon(n_EC_E) %>% deg()) # reassuringly gives (0, 0)
#> [1] 0 0

# ([less] too?) simple test
A1 <- c(-10, -10)
n_EA1_E <- lat_lon2n_E(rad(A1[1]), rad(A1[2]))
A2 <- c(10, 10)
n_EA2_E <- lat_lon2n_E(rad(A2[1]), rad(A2[2]))
B1 <- c(10, -10)
n_EB1_E <- lat_lon2n_E(rad(B1[1]), rad(B1[2]))
B2 <- c(-10, 10)
n_EB2_E <- lat_lon2n_E(rad(B2[1]), rad(B2[2]))

n_EC_E <- great_circle_intersection(n_EA1_E, n_EA2_E, n_EB1_E, n_EB2_E)
(C <- n_E2lat_lon(n_EC_E) %>% deg()) # reassuringly gives (0, 0)
#> [1] 0 0

Created on 2020-11-20 by the reprex package (v0.3.0)

Created on 2020-11-20 by the reprex package (v0.3.0)

Created on 2020-11-20 by the reprex package (v0.3.0)

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions