invgeod
dist, az1, az2 = invgeod(lonlat1::Vector{<:Real}, lonlat2::Vector{<:Real}; proj::String=““, s_srs::String=”“, epsg::Integer=0, backward=false)
Solve the inverse geodesic problem.
Args:
lonlat1: - coordinates of point 1 in the given projection (or a matrix with several points).lonlat2: - coordinates of point 2 in the given projection (or a matrix with same size aslonlat1).projors_srs: - the given projection whose ellipsoid we move along. Can be a proj4 string or an WKT.epsg: - Alternative way of specifying the projection [Default is WGS84].backward: - Iftrue, return backard azimuths.
Returns
dist - A scalar with the distance between point 1 and point 2 (meters). Or a vector when lonlat1|2 have more than one pair of points.
az1 - azimuth at point 1 (degrees) ∈ [-180, 180)
az2 - (forward) azimuth at point 2 (degrees) ∈ [-180, 180)
Remarks:
If either point is at a pole, the azimuth is defined by keeping the longitude fixed, writing lat = 90 +/- eps, and taking the limit as eps -> 0+. ““” invgeod(lonlat1::Vector{<:Real}, lonlat2::Vector{<:Real}; proj::String=““, s_srs::String=”“, epsg::Integer=0, backward::Bool=false) = invgeod([lonlat1[1] lonlat1[2]], [lonlat2[1] lonlat2[2]]; proj=proj, s_srs=s_srs, epsg=epsg, backward=backward) function invgeod(lonlat1::Matrix{<:Real}, lonlat2::Matrix{<:Real}; proj::String=”“, s_srs::String=”“, epsg::Integer=0, backward::Bool=false) @assert (size(lonlat1) == size(lonlat2) || size(lonlat2,1) == 1)”Both matrices must have same size or second have one row” proj_string, projPJ_ptr, isgeog = helper_geod(proj, s_srs, epsg) (!isgeog) && (lonlat1 = xy2lonlat(lonlat1, proj_string); lonlat2 = xy2lonlat(lonlat2, proj_string)) # Convert to geog first d::Vector{Float64} = Vector{Float64}(undef, size(lonlat1,1)::Int) az1::Vector{Float64} = Vector{Float64}(undef, size(lonlat1,1)::Int) az2::Vector{Float64} = Vector{Float64}(undef, size(lonlat1,1)::Int) dist, azi1, azi2 = Ref{Cdouble}(), Ref{Cdouble}(), Ref{Cdouble}() for k = 1:size(lonlat1,1)::Int kk::Int = (size(lonlat2,1)::Int == 1) ? 1 : k # To allow the “all against one” case ccall((:geod_inverse, libproj), Cvoid, (Ptr{Cvoid},Cdouble,Cdouble,Cdouble, Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}), pointer_from_objref(get_ellipsoid(projPJ_ptr)), lonlat1[k,2], lonlat1[k,1], lonlat2[kk,2], lonlat2[kk,1], dist, azi1, azi2) d[k], az1[k], az2[k] = dist[], azi1[], azi2[] end proj_destroy(projPJ_ptr) if (backward) if (isa(az1, Real)) az1 += 180; (az1 > 180) && (az1 -= 360) az2 += 180; (az2 > 180) && (az2 -= 360) else az1 .+= 180; az2 .+= 180 [(az1[k] > 180) && (az1[k] -= 360) for k = 1:lastindex(az1)] [(az2[k] > 180) && (az2[k] -= 360) for k = 1:lastindex(az2)] end end return size(d,1) == 1 ? (d[1], az1[1], az2[1]) : (d, az1, az2) end