Completing graph edges “around the world” to GeoPlot an itinerary

Question

My central goal is to plot a complex international inter-city itinerary on a world map, which might include repetitions of a single segment (e.g., New York to London or London to New York). The original difficulty stemmed from the fact that GeoPath plotted every such "equivalent" segment atop each other, obscuring the number of such passages and confusing the viewer.

@Jose's clever solution to this problem was to use the internal routines of Graph to, in essence, plot a MultiGraph of the route atop a map, where the edges for an equivalent segment would be automatically displaced.

This works well (at least for spatially localized itineraries), but reveals or introduces two new problems when plotting a world itinerary:

• The graph path doesn't "know" the shortest route between cities, and thus sometimes renders an incredibly long segment when the actual shorter segment went "on the other side of the world."
• The edge rendering isn't smart enough to complete an edge that is broken because it goes "on the other side of the world."

See the example map below. Indeed, there are separate edges for the Cape Town to Beijing (as desired), but alas only the portion of the edge near Cape Town is shown--not the portion near Beijing. Moreover, the code doesn't "know" that the proper (shortest) route from San Francisco to Hong Kong is across the Pacific Ocean.

Here is my modification of @Jose's code, where I've used DirectedEdge (to show direction of travel) and where I've colored the segments sequentially by Hue (to help the viewer understand the sequence of the route) and used Partition to speed the entry of cities on an itinerary:

ents = Partition[
      {Entity["City", {"SanFrancisco", "California", "UnitedStates"}], 
       Entity["City", {"HongKong", "HongKong", "HongKong"}], 
       Entity["City", {"Beijing", "Beijing", "China"}], 
       Entity["City", {"CapeTown", "WesternCape", "SouthAfrica"}], 
       Entity["City", {"Beijing", "Beijing", "China"}]}, 2, 1];

map = GeoGraphics[ents, 
         GeoProjection -> "Robinson", 
         GeoBackground -> "CountryBorders"];

proj = GeoProjection /. Options[map, GeoProjection];

uents = Union[Flatten[ents]];

coords = GeoGridPosition[GeoPosition[uents], proj][[1]];

graph = Graph[uents,
   MapThread[
    Style, {(DirectedEdge @@@ ents), 
     Hue /@ (Range[Length[ents]]/Length[ents])}],
   VertexCoordinates -> coords,
   EdgeStyle -> Thick];

Show[map[[1]], graph, ImageSize -> 800]

Question

In short, how can @Jose's code be modified to solve the two problems listed above: 1) produce "split edges" that go the shortest distance "around the world," and 2) show both portions of the edges?

Incidentally, I wouldn't mind using a different GeoProjection if it would solve these problems. However, I've tried every such projection and none quite solve the full problem. Perhaps, too, there is a way to embed a Graph on a sphere (the earth) and then convert that through a GeoProjection to a GeoMap.

Here is Mercator and a truncated range, which avoids breaking paths, but retains the problem of an incorrect path (from San Francisco to Hong Kong):

I think this functionality will be of great use to the Mathematica community.

Answer 1

Borrowing my ideas from your linked post:

Options[multiGeoPath] = {"Directed" -> True, "Offset" -> 1000};

multiGeoPath[spec:{{_, _}..}, OptionsPattern[]] :=
  Block[{locs, gpos, spec2, seen, dirQ, offset, cnt, path},
    locs = Union @@ spec;
    gpos = GeoPosition /@ locs;
    spec2 = spec /. Dispatch[Thread[locs -> gpos]];

    Scan[(seen[#] = 0)&, Sort /@ spec2];
    dirQ = TrueQ[OptionValue["Directed"]];
    offset = OptionValue["Offset"];

    Table[
      cnt = (seen[Sort[s]] += 1);
      path = offsetGeoPath[s, (-1)^cnt * Quotient[cnt, 2] * offset];
      If[dirQ, 
        Arrow[path],
        path
      ],
      {s, spec2}
    ]
]

offsetGeoPath[spec_, offset_] :=
  Block[{pts, degs, len},
    pts = Reverse[GeoGraphics`GeoEvaluate[GeoPath[spec]][[1]], {2}];
    degs = Prepend[GeoDirection @@@ Partition[pts, 2, 1] - Quantity[90, "Degrees"], Quantity[0, "Degrees"]];

    len = Length[pts];

    Line @ Table[
      GeoDestination[
        GeoPosition[pts[[i]]], 
        GeoDisplacement[{1000offset Sin[\[Pi] (i - 1)/(len - 1)], degs[[i]]}]
      ],
      {i, Length[pts]}
    ]
  ]

Your example:

ents = Partition[
  {Entity["City", {"SanFrancisco", "California", "UnitedStates"}], 
   Entity["City", {"HongKong", "HongKong", "HongKong"}], 
   Entity["City", {"Beijing", "Beijing", "China"}], 
   Entity["City", {"CapeTown", "WesternCape", "SouthAfrica"}], 
   Entity["City", {"Beijing", "Beijing", "China"}]}, 2, 1];

GeoGraphics[{Thick, Arrowheads[Medium], 
  MapIndexed[{ColorData[112][First[#2]], #1} &, multiGeoPath[ents]]}, 
  GeoRange -> "World", ImageSize -> Large]

Here's the full circuit. We can see here that the offset choices are not optimal, but it's a start:

ents2 = Join[Most[ents], Reverse[Rest[Reverse[ents]], {2}]];

GeoGraphics[{Thick, Arrowheads[Medium], 
  MapIndexed[{ColorData[112][First[#2]], #1} &, multiGeoPath[ents2]]},
  GeoRange -> "World", ImageSize -> Large]