Every aircraft has an elevation ceiling, which is the maximum height at which it can fly. Although most aircraft have an elevation ceiling high enough that their travel paths are not affected by geographical features, some, such as drones, fly low enough that they need to look out for mountains and other elevations. This project developed a microsite which takes two geographical locations and an elevation ceiling as input, and then outputs a map with an approximation of the shortest route around any mountains that may be in the way.

Suppose we want the shortest path from LA to NYC with a ceiling of 7000 ft:

loc1 = Entity["City", {"LosAngeles", "California", "UnitedStates"}]["Coordinates"];
loc2 = Entity["City", {"NewYork", "NewYork", "UnitedStates"}]["Coordinates"];
ceiling = 7000;

The program first creates a "box" with the two locations at opposite corners. We then create an array of the flyable area using GeoElevationData and our ceiling as a threshold:

elevation[x1_, y1_, x2_, y2_, ceiling_] := 
 UnitStep[QuantityMagnitude[
    GeoElevationData[{{x1, y1}, {x2, y2}}, GeoZoomLevel -> 3, 
     UnitSystem -> "Imperial"]] - ceiling]

m = elevation[loc1[[1]], loc1[[2]], loc2[[1]], loc2[[2]], ceiling];

{w, l} = Dimensions[m];

From this, a network is created by connecting 'neighboring' points in our flyable area. We do this by starting out with a GridGraph, connect diagonally adjacent vertices, and then delete the invalid ones:

coordTrans[w_, l_, m_] := 
 Sort[(#2 - 1) w + w - #1 + 1 & @@@ Position[m, 1]]

remove = coordTrans[w, l, m];

networkCreate[w_, l_, g_] := 
 EdgeAdd[g, 
  Join[Table[
    If[Mod[n, w] == 1, Nothing, UndirectedEdge[n, n + w - 1]], {n, 1, 
     w (l - 1)}], 
   DeleteCases[
    Table[UndirectedEdge[n, n + w + 1], {n, 1, w (l - 1) - 1}], 
    Alternatives @@ 
     Table[UndirectedEdge[n, n + w + 1], {n, w , w (l - 2), w}]]]]

network = networkCreate[w, l, GridGraph[{w, l}]];

flyable = 
 Subgraph[network, Complement[Range[l w], remove], 
  VertexCoordinates -> 
   AbsoluteOptions[network, VertexCoordinates][[1, 2, 
    Complement[Range[l w], remove]]], ImageSize -> 1200]

We then add edge weights and find the shortest path from the bottom left corner to the upper right corner:

networkWeighting[UndirectedEdge[v1_, v2_]] := 
 If[Abs[v2 - v1] == 1 || Abs[v2 - v1] == w, 1, Sqrt[2]]

weights = networkWeighting /@ EdgeList[flyable];

flyable = Graph[EdgeList[flyable], EdgeWeight -> weights];

p = FindShortestPath[flyable, 1, l w];

Viola! Here's our path:

coordTransBack[p_, w_, 
  l_] := {Ceiling[#/w], l - ((Ceiling[#/w] - 1) w + l - #)} & /@ p

pathCoords = coordTransBack[p, w, l];

Graphics[{Red, Line[pathCoords]}]

Let's look at this atop a map:

geoPathCoords[loc1_, loc2_, w_, l_, 
  pathCoords_] := {(#2 - 1) (loc2[[1]] - loc1[[1]])/w + 
     loc1[[1]], (#1 - 1) (loc2[[2]] - loc1[[2]])/l + loc1[[2]]} & @@@ 
  pathCoords

geoPath = geoPathCoords[loc1, loc2, w, l, pathCoords];

GeoGraphics[{Thick, Red, GeoPath[geoPath]}]

There's a clear problem here. Once we're past the Rocky Mountains, we should be taking the greater circle path to our destination. Clearly we're not doing this here - this is because we were working with a discrete set of points over lat-lon space. To fix this, we can take a greedy approach and merge 2 consecutive line segments when the merge avoids the invalid regions. To do this, we express our invalid area and path as regions and use RegionDisjoint:

bad = ImageMesh[Image[m], Method -> "Exact", 
   DataRange -> Reverse[Sort /@ Transpose[{loc1, loc2}]]];

linePoints[p1_, p2_] := Table[(1 - t) p1 + t p2, {t, 0, 1, 0.1`}]

checkDisjoint[bad_, geoPath_, i_, j_] := 
 RegionDisjoint[bad, 
  Line[Reverse /@ linePoints[geoPath[[i]], geoPath[[j]]]]]

i = 0;
j = 1;
pathFinal = {};
While[i < Length[geoPath],
 While[checkDisjoint[bad, geoPath, i, j] && j < Length[geoPath],
  j++];
 AppendTo[pathFinal, geoPath[[j]]];
 i = j;
 j = i + 1;
 ]

The new path:

GeoGraphics[{Thick, Red, GeoPath[pathFinal]}]

Much better!

A potential issue is that our solution 'hugs' the mountains too much. In the real world, this would be too dangerous. A way around this is to pad our invalid region array, say using Dilation, before we build our connectivity network.

This algorithm produces a good approximation of the shortest route. However, it fails when our network is disconnected, even though there could be a potential path that goes outside the box. A further continuation of this project could be done to create larger boxes to look through.

Attachments: