Is there a simple / easy way to transfer the longitude and latitude into the x value and y value (cartesian coordinate system) of a point in a planar map?
If you want great circle distance, and are happy to think of the earth as a sphere, then it should just be arccos(dot(v,w)) where v,w are 3-vectors like v=[sin(lat), cos(lat)*sin(long), cos(lat)*cos(long)]… up to units. But perhaps a package called Geodesy will already have this somewhere?
@improbable22, I read the package description of Geodesy and the above example if directly from their git link. But I’m not sure I totally understand it.
Now that I look at the readme, it sounds like “distance” is probably a straight-line distance in 3D. Is this what you want?
They say "Future work may focus on geodesics and related calculations " which is what I was aiming for. But the package is more sophisticated than assuming the world is a sphere, and thus working out the accurate great-circle distance (i.e. geodesic length) accurately would be harder work.
To get the cartesian (2d-) distance on a map, have a look at Proj4.jl which is a Julia wrapper around a relatively widely used C library. As @improbable22 already wrote, this depends on the map projection.
@scelles, there are several options for calculating great circle paths between points in Julia, including bearings:
Straightforward Haversine calculations, as explained above—works for spheres only, but easy to implement oneself;
Vincenty’s algorithm—works for moderately elliptical ellipsoids but has problems with near-antipodal points. A versions is available in Geodesics.jl (code mostly originally from GreatCircle.jl).
Charles Karney’s algorithms—work for even very eccentric bodies and all combinations of points, with no slowdown compared to (2). Implemented in GeographicLib.inverse, or Proj4.geod_inverse (via the Proj4 external library).
For constant bearing it should only be a question of converting to Mercator (or other conformal projection) and compute the conjugate of the of the straight line angle.
compute the conjugate of the of the straight line angle
@joa-quim I’m not sure to understand what you mean about “conjugate”
In Mercator projection, bearing should simply be angle between segment North-Point1 and Point1-Point2.
Is there some Julia packages to show earth (for example) and see how trajectory between two points can be displayed on it and using different map projections?
Well, the angle between 2 points in a cartesian axis (theta) is the angle counterclockwise from the x-axis, but bearing (azimuth) is the angle clockwise from North. So azimuth = 90 - theta
(to plot an loxodrome one could interpolate linearly between two points in a Mercator, convert back to gepgraphic (using mapproject) and add that line to the plot below)
Thanks for your examples with GMT.jl
Being able to also plot loxodrome will be great!
I hope it will be possible to have such plots with GeoMakie.jl also. Pinging @asinghvi17
It should be - if you have the geodesic line, all you need to do is get the coordinates in Mercator space and construct a LinRange for lat and lon - then, you can reproject into any coordinate system.
I haven’t been able to create an example yet, but probably will be able to in the next couple of days.
Turns out that GMT can compute loxodromes directly with the module sample1d (sorry the syntax for this case was not yet expanded to the more comprehensible keywords)
Hey, do you know how to get the geodesic line data (x,y,z)[lon,lat,alt] between two points coordinates?, I would like to have the data and then plot it with whatever library I wish, as you mentioned before.
Here is a translation to Julia of a Matlab code I wrote years ago. It is approximate in that it assumes a spherical earth model, but I found it completely adequate for plotting purposes.
using LinearAlgebra: ⋅, norm, ×
"""
gcinterp(p1, p2, dθ)
Compute list of points on a great circle interpolation assuming a spherical earth model.
### Input Arguments
- `p1`, `p2`: Iterables of length 2 containing the longitude λ and lattitude L
in degrees of the initial and terminal points of the desired arc.
- `dθ`: The desired angular increment in degrees for point spacing along
the interpolated great circle arc.
### Return Value
A vector of ordered pairs `(λᵢ, Lᵢ)` (longitudes and lattitudes in degrees) of
points on the great circle path joining `p1` and `p2`. The arc length spacing of the
points will be equal to or slightly less than `dθ`.
"""
function gcinterp(p1, p2, dθ)
λ₁, L₁ = float.(p1[1:2])
λ₂, L₂ = float.(p2[1:2])
rg1 = ll2rect(λ₁, L₁)
rg2 = ll2rect(λ₂, L₂)
ct0 = rg1 ⋅ rg2
st0 = norm(rg1 × rg2)
rg1 /= st0
rg2 /= st0
θ0 = atand(st0,ct0) # great circle arc length
n = 1+ceil(Int, θ0/dθ)
points = [(λ₁,L₁) for _ in 1:n]
for (i,θ) in enumerate(range(0, θ0, length=n))
rvec = rg1 * sind(θ0-θ) + rg2 * sind(θ)
points[i] = rect2ll(rvec)
end
points[1] = p1[1], p1[2] # Eliminate rounding errors
points[n] = p2[1], p2[2]
points
end
"""
ll2rect(λ, L)
Convert longitude `λ` and lattitude `L` (both in degreees) to a 3-vector of rectangular coordinates.
"""
function ll2rect(λ, L)
sinλ, cosλ = sincosd(λ)
sinL, cosL = sincosd(L)
[cosL * cosλ, cosL * sinλ, sinL]
end
"""
rect2ll(rvec)
Convert rectangular coordinates (a 3-vector of unit length) to a (longitude,lattitude) tuple.
Angular units on output are degrees.
"""
function rect2ll(rvec)
λ = atand(rvec[2], rvec[1])
L = asind(rvec[3])
(λ, L)
end
Off course. I tried to see if there was an option to generate loxodromes that I forgot the obvious.
)