Introduction
Introduction
This project involves determining the least distance path between two coastal points using the Wolfram Language. Even for a human, this optimal path is not necessarily trivial as several questions arise when trying to evaluate this path, such as which direction(s) to head in, when to hug the coastline, when to cut across the coastline, etc. Additionally, a shortest path on a map (a 2D projection of the globe) might not represent the actual shortest path on a 3D globe.
To solve this problem, I will first go over how to determine this shortest path, and how I obtained this path using Wolfram. Then, I will apply this logic to a relatively small island nation before demonstrating it on a global scale.
To solve this problem, I will first go over how to determine this shortest path, and how I obtained this path using Wolfram. Then, I will apply this logic to a relatively small island nation before demonstrating it on a global scale.
Methods
Methods
Jamaica
Jamaica
To see how to find the shortest path between two coastal points, lets look at the coastline of Jamaica below, where the two green disks are the points we wish to find the shortest distance between:
In[]:=
JamaicaPoints=LastSortCountryData,{"Shape","Mercator"}[[1,3,1]];SeedRandom[9];randompoints=RandomSample[JamaicaPoints,2];disks=Sequence[Green,Disk[#,.02]&/@randompoints];JamaicaCoastline=Graphics[{Line[JamaicaPoints],disks}]
Out[]=
To find the shortest distance, we can visualize a string at the bottom point going through the top point. If we pull on the string as tight as we can, the string will then become the shortest path between the two points that does not go through land. This path will consist of several straight lines.
Originally, when I tried to generate this distance, I thought of designing an iterative algorithm to check all possible slopes of lines coming from a point between two bounds. These bounds would be the lines formed by the adjacent points on the coastline. In the diagram below the red lines formed by the blue and yellow points are the bounds to check for the green point.
Originally, when I tried to generate this distance, I thought of designing an iterative algorithm to check all possible slopes of lines coming from a point between two bounds. These bounds would be the lines formed by the adjacent points on the coastline. In the diagram below the red lines formed by the blue and yellow points are the bounds to check for the green point.
In[]:=
SeedRandom[10];shortcoastline=JamaicaPoints[[1;;40]];randompoints=RandomSample[shortcoastline,2];
In[]:=
possibleSailDirections[startpoint_,endpoint_,coastline_](*given a startpoint, endpoint, and list of points which is the coastline*):=Module[(*left and right direction limits are the points on the coastline which are just before and after the startpoint*){leftdirectionlimit=coastline[[Position[coastline,startpoint][[1,1]]-1]],rightdirectionlimit=coastline[[Position[coastline,startpoint][[1,1]]+1]]},(*make a graphic showing halflines of the right and left direction limit*)Graphics[{Line[coastline],Red,HalfLine[{startpoint,rightdirectionlimit}],Red,HalfLine[{startpoint,leftdirectionlimit}],Green,Disk[startpoint,.001],Blue,Disk[leftdirectionlimit,.001],Yellow,Disk[rightdirectionlimit,.001],Red,Disk[endpoint,.001]},ImageSize->Large]]
In[]:=
possibleSailDirections[randompoints[[1]],randompoints[[2]],shortcoastline]
Out[]=
However, I found out that there was a built-in Resource Function in Wolfram called RegionFindShortestPath, which will determine the shortest path between two points in a Mesh Region (this is the ocean in our case). I won't go in depth about the function in this essay, but essentially it will find all possible paths between a given point and inward facing points (from the perspective of the region) in the region. Going back to the string analogy, these inward facing points (the points projecting into the ocean) are where the string gets "stuck," and the shortest path will consist of lines between such points.
My mentor and I were able to create a Mesh Region of the ocean surrounding Jamaica with the following function:
My mentor and I were able to create a Mesh Region of the ocean surrounding Jamaica with the following function:
In[]:=
JamaicaMesh=DiscretizeRegion[RegionDifference[Rectangle[{Min[Flatten[JamaicaPoints,{1}][[All,1]]]-.0001,Min[Flatten[JamaicaPoints,{1}][[All,2]]]-.0001},{Max[Flatten[JamaicaPoints,{1}][[All,1]]]+.0001,Max[Flatten[JamaicaPoints,{1}][[All,2]]]+.0001}],Polygon[JamaicaPoints]]]
Out[]=
Then, I was able to create an image that dynamically finds the shortest path between two points in this mesh region.
func=Manipulate[Show[JamaicaMesh,Graphics[{Thick,Red,Dynamic@Line@func[p1,p2]}]],{{p1,randompoints[[1]]},Locator},{{p2,randompoints[[2]]},Locator},SaveDefinitions->True]
Out[]=
However, this path is not necessarily the shortest path between the two points on a 3D globe, just on a 2D projection. In order to minimize this distance on a 3D projection, I used Wolfram's GeoPath function, which can find the shortest path between two points on the globe with their coordinates as input (the example below uses different points than the example above). I had to extract the coastline of Jamaica using "Coordinates" because Wolfram's "Shape" function does not use the standard latitude-longitude system, which the GeoPath function needs.
In[]:=
coords=FlattenReverseCountryData,"Coordinates",2,1;randcoords=RandomSample[coords,2];RandomSeed[4];coordsmesh=DiscretizeRegion[RegionDifference[Rectangle[{Min[Flatten[coords,{1}][[All,1]]]-.0001,Min[Flatten[coords,{1}][[All,2]]]-.0001},{Max[Flatten[coords,{1}][[All,1]]]+.0001,Max[Flatten[coords,{1}][[All,2]]]+.0001}],Polygon[coords]]];
In[]:=
func2=ResourceFunction["RegionFindShortestPath"][coordsmesh];Iconize[func2,"func2"]GeoGraphics[GeoPath[Reverse[func2[randcoords[[1]]][randcoords[[2]]],3],"Geodesic"],GeoRange->{{17.5,19},{-78.5,-76}}]
Out[]=
Out[]=
Applying to the World
Applying to the World
When trying to generalize this algorithm to the globe, I had several problems. Wolfram's built-in RegionFindShortestPath function is relatively inefficient since it calculates all paths from a given point to inner points in a mesh region. Although it worked for Jamaica, it took too long to run on a detailed map of the world. We tried to fix this by taking the Schematic Polygon (basically a less precise version of a country's border) of the world:
In[]:=
img=Graphics[CountryData["World","SchematicPolygon"]]
Out[]=
However, when running the previous function to generate a mesh region of the ocean, the world still had too many points, and my kernel crashed when trying to run it. Then, I realized that some of the Schematic Polygons for certain small archipelagos have several thousands of points in their Schematic Polygon, and they were the ones accounting for the majority of the points in the above map. I fixed this by manually removing these smaller archipelagos (as they are not relevant for this project). Additionally, we found another, more simpler, way to generate the Mesh Region using Image Processing:
With this, we were able to apply the dynamic shortest path on a global scale:
There were still a couple things wrong with this map though. First, the distance between Spain and Morocco is so small, that the Triangulate Mesh Region did not have any points in between them, so the Strait of Gibraltar was treated as a connected piece of land. Second, the Suez Canal was not treated as a body of water either. Third, optimal paths that travel across the Pacific Ocean were not recognized since this was a 2D map (for example, the path from Los Angeles to Shanghai would travel through the Indian and Atlantic Oceans rather than the Pacific).
To fix the first two problems, I modified the SchematicCoordinates of Spain, Morocco, and Egypt to make the body of water between them larger and more prevalent, allowing the function to calculate paths through the Strait of Gibraltar and the Suez Canal.
To fix the first two problems, I modified the SchematicCoordinates of Spain, Morocco, and Egypt to make the body of water between them larger and more prevalent, allowing the function to calculate paths through the Strait of Gibraltar and the Suez Canal.
For reference, this is what Spain and Morocco looked like before modification:
Whereas this is what they looked like after:
To solve the Pacific Ocean problem, I stitched two copies of the map together, and then selected a certain region of it to form a new map as shown:
Then, I ran the RegionFindShortestPath on this map, and got this result:
Conclusion
Conclusion
This project went in a direction that I did not expect it would before starting. I didn't think I would be able to get a dynamic map that finds the shortest distance between coastal points. I also thought that I would have to take other geographical factors into account (currents, winds, etc.) when evaluating the shortest path.
Future Directions
Future Directions
As aforementioned, I think it would be interesting to see how much (if, at all) this path changes taking these factors into consideration. Also, with a more powerful computer, I would like to see this shortest distance calculated on a more precise map of the globe. I wasn't able to incorporate the GeoPath into a dynamically updated map, so that would be a good extension as well.
Acknowledgements
Acknowledgements
I would like to thank my mentor, Stephen Schroeder for helping me with this project. It was very challenging, long process and I was inexperienced with the language, but he helped me significantly. I would also like to thank TAs Drake Hayes and Anika Patel for helping me throughout camp, and my mentor group, and the rest of the staff and students for making camp enjoyable!