Oblique Mercator for straight line between points as map's equator

I’m trying to create a map with the line described here (the longest straight line only in the oceans) as the “equator” in an Oblique Mercator projection.

I’ve tried different combinations of the two points given, but I’m obviously not understanding something. What is the input supposed to be?

import pygmt

fig = pygmt.Figure()
# Using the origin projection pole
    region=[0, 360, -80, 80],
1 Like

Very cool.

Cant help you with PyGMT, but here’s a way to plot it in ShellGMT:

cat <<eof > points
66:40 25:17
162:14 58:37

gmt begin longest-sailable-line-path png
gmt set GMT_THEME=minimal
gmt pscoast -Gblack -JM20c -R-180/180/-70/70 -Dc -A10000
gmt plot points -Sc0.1c -Gred -W
gmt project -C66:40/25:17 -E162:14/58:37 -G10 -Q -L-32090/0 | gmt plot -Wblue -Bxafg -Byafg -B+t"Longest Sailable Straight Line Path on Earth"

gmt end show

rm points

1 Like

Good catch on that -L option. By reading its description I couldn’t tell that it what’s need to plot the longest part of the great circle. Though there is a bit of un-mentioned work as you had to compute the longest distance on a separate step.

Yes, I was few movements away from creating a post ala “How do I draw the longest part of a great circle”, but by looking at some project examples I plunged into something that looked at least looked right. The necessary distance I just read from the article (cheating…).

By the way, what is the best way of finding the longest part of a great crircle? :slight_smile:

I think we have to compute the total length and subtract the shortest path.
At least, that’s what I intend to do when I write the GMT.jl gcirc() function that will work for the ellipsoid as well (based on proj4, project is spherical only).

I normaly have these ideas in the shower but this time it was with a glass of wine.

We don’t need -L. We can plot from point to it’s antipode and from there to destination.

Thank you for the responses. I’m fine with using the CLI (though I haven’t had the chance to try it yet), and plotting the line will help diagnose, but this didn’t really answer my question: I’m trying to project an Oblique Mercator where that line is flat in the middle.

The best I can do is

gmt begin map
	gmt coast -R-180/180/-40/40 -JOb66:40/25:17/162:14/58:37/20c -Bafg -Gred
	gmt project -C66:40/25:17 -E162:14/58:37 -G10 -Q -L-32090/0 | gmt plot -W1p,blue
gmt end show

I am not able to fiddle with -R to get the line centered in the plot though. Perhaps it is possible but gave up after 5 minutes.

I think we need to use the equivalent -JOa or -JOb version instead. So compute the oblique pole and use a point mid-way long the path as the origin and pole.

Taking the cross-product of your two points with gmt vector -Tx gives the pole, and the half-way point on @Andreas line results in a working projection:

-R-180/180/-40/40 -JOc-82.1244285454/-52.0295503563/-38.9210614494/29.6353788676/20c

1 Like

I was able to find an ellipsoidal (geodesic) alternative. But with the points found for the great circle case a geodesic won’t work. See details at https://www.generic-mapping-tools.org/GMTjl_doc/tutorials/longest_sail/longestsail/ @WHFSmith please see if I’m not saying nonsense. This geodesics story is quite confusing.


@Joaquim I do not think you are saying nonsense, but one of the ways in which I now feel my brain getting old is that it takes me more time than I can give this problem to understand issues in differential geometry. I am intrigued by the problem, and particularly by your statement that swapping the roles of the start and end points gives a different answer (should this be true or is it a bug?). The following remarks are completely useless and do not present a solution, but maybe they will help someone smarter than me (or more handy with differential geometry) help us think it through.

First, I note that the term “longest straight line path” in https://arxiv.org/pdf/1804.07389.pdf is not really properly defined, as any path that is contained in the surface of a sphere or an ellipsoid and which is not of length zero is not contained in a straight line. On a sphere, the shortest path between two distinct points is a great circle and so is contained in a plane; what the problem as originally posed seemed to have in mind is the portion of the great circle defined by two distinct points which is not the shortest path but the remainder of the great circle after the shortest path is removed. I will call this the “longer way around” path. Because it is in the same plane and the same great circle as the shortest path, it is unambiguously defined (assuming that the two points are not antipodal): as a point moves along the shortest path in the sphere the antipode of that point traces the “longer way around” path on the sphere.

Given two (non-antipodal and distinct) points on the surface of an ellipsoid, there is only one unique path between them that is the shortest possible path. But is that sufficient to determine uniquely the ellipsoidal analog of the spherical “longer way around” path? How? Is the shortest ellipsoidal surface path contained in the intersection of the ellipsoid with a plane? (I don’t think so, but these days I have to drink decaffeinated coffee…) As a point moves along the shortest path on an ellipsoidal surface, does its antipodal point trace out the “longer way around” path on the ellipsoid? (I think it is clear what we mean by antipode on an ellipsoid.)

So how do we define the “long way around” on an ellipsoid? Note that it is not enough to say that we want the longest path between two points contained in an ellipsoid. To see this, imagine starting at the North Pole and taking a spiral path so that each time you have moved through 360 degrees of longitude you have moved a small amount southward. If you move an infinitesimally small amount southward with each rotation of longitude then you will eventually arrive at the South Pole but by an infinitely long path.

I am particularly intrigued by @Joaquim 's comment that he can find a solution that starts at point A and ends at point B, but if he swaps the positions of A and B then the path fails to be a solution. It seems intuitive (though differential geometry on strange surfaces is not intuitive to me; my intellect is infinitely smaller than Gauss’s) that the shortest path from A to B and the shortest path from B to A must be the same path, apart from a reversal of the direction of traverse. So, should we expect that our definition of longest “longer way around” path should also have that property (that swapping A and B doesn’t change the path)? If yes, then is there a bug in the code @Joaquim is using? Or is the problem that the longest path is so ill-posed that even an infinitely perfect coder cannot write an algorithm to compute it?

Is it possible that someone has defined the “longer way around” path as follows:
“Given two distinct and non-antipodal points A and B on an ellipsoid,
Find the shortest path from A to B, and the azimuth of this path at A; call this azimuth Alpha.
Then define the “longer way around” path as a path that departs from A with azimuth Alpha + Pi.”
Is this definition even workable? Would this definition arrive at B? Would it explain Joaquim’s paradox?

Maybe this is another way to define it:
Given two distinct and non-antipodal points, A and B, and with longitude defined to be periodic as usual, consider that A and B divide the 2 pi period of longitude into a larger and a smaller interval. The shortest path from A to B goes through each longitude in the smaller interval only once; that is, for each longitude in the smaller interval there is one and only one latitude at that longitude such that that lat,lon pair is a point on the shortest path. Now define the “long way around” path so that it goes through each longitude in the larger interval, but otherwise behaves like a geodesic path (in the surface of the ellipsoid, continuous, differentiable, with one and only one latitude occupied at each longitude in the set).
Is this the definition we have in mind? If so, is the swap(A,B) paradox correct and expected, or not?

Sorry I cannot offer solutions. I am eager to hear more about this from someone who is good at differential geometry on ellipsoids. (And please explain in enough detail so that non-experts can understand.)

Sorry. My last conjectured definition fails if points A and B are on the same meridian of longitude. I was assuming A and B had distinctly different longitudes. Obviously there are a lot of special cases that a smart algorithm would need to handle correctly.

Hi Walter,

I’m glad that those doubts on the ‘going around’ case are not only mine. When I fist started to do some calculations I was persuaded that I must be doing some mistake (which is not yet proven otherwise). The case that one cannot reach the starting point when going around the Earth following the shortest path between A and B, for the particular case that A = B, still make me feel that that there stuff in here that I don’t understand. There might be bug(s) in the code but I used the geod and invgeod functions from geographiclib (via PROJ) and those are considered top quality for these type of calculations. For interested readers I wrapped them in this Julia small function that can be read as a kind of pseudo-code.


Regarding this question,

“Is the shortest ellipsoidal surface path contained in the intersection of the ellipsoid with a plane? (I don’t think so, but these days I have to drink decaffeinated coffee…)”

No, the intersection of a plane and the ellipsoid is a great ellipse and that is not the shortest path on the ellipsoid (see this from the GeographicLib manual for a discussion on that regard)

Regarding the paradox,

define the “longer way around” path as a path that departs from A with azimuth Alpha + Pi.” Is this definition even workable? Would this definition arrive at B? Would it explain Joaquim’s paradox?

I’m afraid not. See the three figures below. Note that the firs two are exactly equal.

“I am eager to hear more about this from someone who is good at differential geometry on ellipsoids. (And please explain in enough detail so that non-experts can understand.)”

I’ll take the risk of being inconvenient by pinging Charles Karney @cffk that I’m sure he will have something to teach us.

Here’s my take on these issues:

@WHFSmith On the question of the definition of the “longest straight line path”: “Straight” here should be thought of in the 2d world (without regard to how the 2d surface is embedded in 3d) and it means having zero geodesic curvature. (Perhaps “straightest” would be a better term that “straight”.) It is the path a ship would take with rudder amidships, a plane in straight and level flight, a car with the steering wheel locked in the center position, a model train track without the curved pieces, etc.

@Joaquim On “reverting the order of these two points would not have worked”: I don’t understand this at all. You can reverse a geodesic path very easily. For the case in question we have:

$ GeodSolve -: -p 9 -b <<EOF
25:14:34.6N 066:42:5.8E 213:28:14.2596640340 32065757.768734273
58:27:27.6N 162:12:5.8E 107:51:16.3789997267 32065757.768734273
58:27:27.6000000000N 162:12:05.7999999998E 107:51:16.3789997267
25:14:34.6000000000N 066:42:05.8000000003E 213:28:14.2596640340

@Andreas On finding the geodesic “the long way round”: If the shortest geodesic between A and B has length f*C where C is the approximate circumference of the earth, then what is the (non-shortest) geodesic with approximate length (1-f)*C? This has a unique answer unless f is close to zero or half. In these latter cases, the complication is that there are 3 possible solutions. These paths can all be found relatively easily (generalizing the method used to find the shortest geodesic).

By the way, you can sail a long way from Plymouth, England, on heading 210. If you can fiddle the path to go between the N and S islands of New Zealand, you can make it to Japan.

@cffk First, thank you very much for having looked at this interesting and, for me, still confusing problem.

I am sorry but there is one thing that I don’t understand in your example. From my understanding, to go from A to B from the longest path you compute the backward azimuth from A to B at A and depart from there with that backward azimuth. Next you do the same thing by swapping A and B. But when I try to confirm this I can’t recover the angle 107:51:16.3789997267

Since we need to compute backward azims it’s simple if we use decimal degrees

echo 25:14:34.6N 066:42:5.8E 213:28:14.2596640340 | gmtconvert
25.2429444444 66.7016111111 -146.529372316 (= 213.36715905483763)

echo 58:27:27.6N 162:12:5.8E 107:51:16.3789997267 | gmtconvert
58.4576666667 162.201611111 107.854549722

Compute the backward azimuths at the departing point.

$ echo 25:14:34.6N 066:42:5.8E 58:27:27.6N 162:12:5.8E | GeodSolve.exe -p 9 -i
33.36715905486175 108.33445121260358 7958785.658131862

echo 58:27:27.6N 162:12:5.8E 25:14:34.6N 066:42:5.8E | GeodSolve.exe -p 9 -i
-71.66554878739642 -146.63284094513824 7958785.658131862

First azim pops up fine like in your example

33.36715905486175 + 180 = 213.36715905486176

But the second one no

-71.66554878739642 + 180 = 108.33445121260358
108.33445121260358 != 107.854549722

Trying to do the reverse exercise with the azimuths obtained here do not show a reversible path

echo 25:14:34.6N 066:42:5.8E 213.36715905483763 32065757.768734273 | GeodSolve -: -p 9 -b
58:33:02.2083714721N 162:15:31.8312587573E 107:52:16.5265922991

$ echo 58:27:27.6N 162:12:5.8E 108.33445121260358 32065757.768734273 | GeodSolve -: -p 9 -b
25:29:38.4826548377N 066:17:00.3222169120E 213:26:42.8930980065

This is exactly the situation represented above by the last two figures of my last example of three.

So, my question is: how was the azimuth of 107:51:16.3789997267 obtained?

@Joaquim My azimuths for the long way round were

at A 213:28:14.2596640340 = 213.47062768445389 = -146.52937231554611
at B 107:51:16.3789997267 = 107.85454972214630

(You mistakenly convert my azimuth at A to 213.36715905483763.) The azimuths for the short way round are

at A 33.36715905486175
at B 108.33445121260358

The discrepancies between the short way and long way azimuths are then

at A 33.36715905486175+180-213.47062768445389 = -0.103
at B 108.33445121260358-107.85454972214630 = 0.480

(If you use your incorrectly converted azimuth, the discrepancy at A is indeed tiny: 33.36715905486175+180-213.36715905483763 = 0.00000000002.) Admittedly the discrepancy at A is smaller that at B; but I wouldn’t read too much into this. Basically, finding the long way round path by reversing the azimuths for the short way round path is always going to be approximate.

I found my long way round path by doing adjusting the distance and the azimuth for a geodesic starting at A so that the endpoint is at B. I did this 2d optimization problem graphically in Octave (taking about 2 iterations for each added digit). (I used the my native Octave implementation of the geodesic problem.)

The right way to do this is as follows: Start at the most poleward point (in Kamchatka) and a starting guess for the azimuth (e.g., the reversed short way round azimuth). Solve the hybrid problem finding the longitude at the second intersection of the path with the circle of latitude for the other point. This is a straightforward (non-iterative) calculation. Now apply Newton’s method to adjust the initial azimuth so that the longitude at the second point matches the given one. There’s a simple formula giving the derivative dlon2/dazi1 in terms of the reduced length m12. After a 2-3 iterations you’ll have full accuracy. (This just mimics the solution for the shortest geodesic for which you would pick the first intersection with the circle of latitude.)

Sorry for the mistake in first azim (and have done those conversion so many times).

A thanks for the further explanation. I think I’ll be able to come out with something similar to what you described. I’m using geod_direct and geod_inverse via PROJ (it exports the C wraps of them) but that doesn’t give access to m12.

So the take out points of this are:

  • Finding the long way round path by reversing the azimuths for the short way round path is always going to be approximate.

  • If one do that, we fall on the paradox of the non-reversibility.

  • But there is a solution

As a note on another different subject, and taking advantage of your attention. It would be lovely if more GeographicLib function were included in PROJ. I’m thinking in particular the functions that are used in the gravity program (then I could use them in GMT.jl :slight_smile: ). And it looks to me that would be within the domain of interest of PROJ users.

And once more, tanks for your contribute on this. (I’ll try to look at the Bristol route too).

As a follow up and hopefully settle it, I was able to devise an approach that obtains the reversible geodesic with a single step, no iteration needed (at least for this case). Details for this are better explained in the code itself

With it we can get a reversible geodesic up to the decimeter level and the angles are very close to those reported in Charles’ post.

A Julia example

ptPak = [66.7016111111 25.2429444444];
ptKam = [162.201611111 58.4576666667];

D = geodesic(ptPak, ptKam, longest=true);
println(D[end, 1], "\t", D[end, 2], "\t", D[1,3])
162.2016143757807       58.4576723686032        -146.52937862182463

D = geodesic(ptKam, ptPak, longest=true);
println(D[end, 1], "\t", D[end, 2], "\t", D[1,3])
66.70161540794047       25.24294132647047       107.85454447496345

I’ve updated the GMT.jl docs example and show here the new ~longest route.

Getting the Gravity stuff into PROJ will be a heavy lift, I fear. There’s a fair amount of work on my part required to get the classes extricated from the rest of GeographicLib. Then there’s the question of whether this mass of code really belongs in PROJ. My instinct would be to say that it doesn’t; everything else in PROJ is 2d (maybe with 1 or 2 trivial exceptions).

By the way, Plymouth is a better starting point for a long sailing route compared to Portsmouth (or Bristol). However, I don’t see how to avoid hitting the Falklands or New Zealand… But I didn’t try very hard.