I visited the furthest Pole and placed a register there: http://sangabrielmnts.myfreeforum.org/about7040.html. Do visit and sign in!
So I was out hiking with a friend, and a question came up about where the least accessible point of the San Gabriel Mountains was, with "accessible" defined as having a road or trail nearby. I decided to answer this question. The resulting code is in a fresh repository: https://github.com/dkogan/inaccessibility.
This is called the Pole of Inaccessibility: a point that is as far away as possible from a given set of objects. Locations of such poles are known for the most landlocked spot on earth or most far away from land. Here we limit ourselves to the San Gabriel Mountains, and try to stay away from roads and trails.
Input data processing
OpenStreetMap has open data I can use to map out all the roads and trails. This is the input dataset.
For 2D geometry, the best approach to compute the Pole of Inaccessibility appears to be to construct a Voronoi diagram of the geometry we're trying to stay away from, and to find the Voronoi vertex corresponding to the furthest-away point.
Our world is not 2D. Instead, it has varying elevation sitting on top of an ellipsoid. The grand purpose here is to compute a location that hardy people can visit and to tell everybody they did it, so extreme accuracy is not required. Thus I claim that assuming the world is locally-flat and using the Voronoi-diagram-based method is sufficient. So I construct a plane that best describes my query area and project all my input points to this plane. I use a plane that is tangent to the Earth's surface at the center of the query area. This clearly wouldn't work if trying to find the pole of inaccessibility of something as large as an ocean, for instance, but it works here.
To compute the tangent plane, I assume the Earth is spherical. As I move along the tangent plane away from the point of tangency, the elevation error grows:
E = sqrt(Rearth2 + d2) - Rearth
The San Gabriels are about 80km across, and the tangent plane sits in the middle, so at worst d = 40km and the error is about 125m. That's plenty good enough. Plot (source):
I ignore the ellipsoid shape of the Earth outright. I ignore the topography as well, since including it in my distance metrics would require a fancier algorithm than making a Voronoi diagram, and it would make the notion of "inaccessibility" more ambiguous.
Pole of Inaccessibility computation
I want to use the most basic Voronoi algorithm, so I represent my input as a set of points only; no line segments. To get reasonable accuracy, I make sure to sample each road at least every 100m.
Now that I have my set of dense-enough points in 2D, I construct the Voronoi diagram. Without constraints the furthest-away point would be infinitely far off to one side, so generally people constrain the solution to lie within the convex hull of the input points. This means that the Pole of Inaccessibility lies either on a Voronoi vertex or at an intersection of a Voronoi edge and the convex hull of the input. In my case there are generally more roads at the edges of my query area that in the interior (less stuff in the mountains than in the flats), so I simply assume that the Pole of Inaccessibility is not on the convex hull. This simplifies my implementation since I simply ignore all the Voronoi vertices that are outside of the query region.
So I need to look at every Voronoi vertex, check the distance between it and an adjacent input point, and return the vertex with the largest such distance.
Each step in the process lives in its own program. This simplifies implementation and makes it easy to work on each piece separately.
First we query OSM. This is done with the
query.sh script. It takes in corners
of the query area, constructs the query, sends it off to the server, and stores
query.sh takes 4 arguments; lat0, lon0, lat1, lon1, and stores its
output in a file called
query_$lat0_$lon0_$lat1_$lon1.json. The query uses the
OSM Overpass query language. By default I simply look at all the roads, trails
(everything with a
[out:json]; way ["highway"] ($lat0,$lon0,$lat1,$lon1); (._;>;); out;
If I want to only consider roads in my computation (allow trails), then I can exclude trails from the query:
[out:json]; way ["highway"] ["highway" != "footway" ] ["highway" != "path" ] ($lat0,$lon0,$lat1,$lon1); (._;>;); out;
$ ./query.sh 34.1390884 -118.4944153 34.5020298 -117.5852966 % Total % Received % Xferd Average Speed Time Time Time Current Dload Upload Total Spent Left Speed 100 28.3M 0 28.3M 0 138 72803 0 --:--:-- 0:06:48 --:--:-- 138k $ ls -lh query* -rw-r--r-- 1 dkogan dkogan 29M May 6 04:12 query_34.1390884_-118.4944153_34.5020298_-117.5852966.json
Next, I take the lat/lon pairs, map them to the tangent plane and make sure the
data is sufficiently dense. This is done by the
massage_input.pl script. It
takes in the
query_....json file we just obtained, and generates a
points_$lat0_$lon0_$lat1_$lon1.dat file that is a list of (x,y) tuples in my
plane. There's a small header of 4 values, representing the bounds of my data so
that I can reject outlying vertices, as described earlier.
$ ./massage_input.pl query_34.1390884_-118.4944153_34.5020298_-117.5852966.json $ ls -lh points* -rw-r--r-- 1 dkogan dkogan 4.3M May 6 04:20 points_34.1390884_-118.4944153_34.5020298_-117.5852966.dat
Pole of Inaccessibility computation
Now we can compute the Voronoi diagram. I use boost::polygon to do this. I had concerns that this step would be prohibitively slow, but the algorithm and this implementation are quick-enough such that this "just works".
points_....dat file is inputs on standard input. Note that this is
different from the other tools that read a file on the commandline instead.
For each Voronoi vertex I get an arbitrary neighboring edge, and an arbitrary neighboring cell. The distance between the vertex and the cell center is identical for any such edge, cell by definition of a Voronoi vertex. I keep track of the cell with the largest distance between the vertex and the cell center, and I report the vertex with the largest such distance as my Pole of Inaccessibility.
$ ./voronoi < points_34.1390884_-118.4944153_34.5020298_-117.5852966.dat furthest point center, surrounding points: 25541 -78 25308 4259 21223 -543 26931 -4192 distance: 4342.873206
Bam! So the Pole of Inaccessibility is about 4.3 km from the nearest trail/road. The coordinates here are in my 2D tangent plane, which isn't super useful. Now I convert them to lat/lon and I'm done.
Unmapping the planar coordinates
I do this with the massaging script as before simply by passing the coords in on the commandline:
$ ./massage_input.pl query_34.1390884_-118.4944153_34.5020298_-117.5852966.json 25541 -78 25308 4259 21223 -543 26931 -4192 34.3206972918426,-117.761740671673 34.3597096140465,-117.764149633831 34.3165155855012,-117.808770212857 34.2837076589351,-117.746734473897
I did this three times
- avoiding all roads, trails
- avoiding all roads
- avoiding all paved, driveable roads
All Poles of Inaccessibility are above the East Fork of the San Gabriel River, by Ross Mountain and Iron Mountain:
|pole||lat,lon||distance to nearest (m)|
|roads only||34.3108, -117.7474||5677|
|paved roads||34.2992, -117.7167||8157|
This is shown nicely on the map:
Looks like the bounding spots for the roads,trails point are the road up to South Mt Hawkins, the PCT on top of Mt. Baden-Powell and the trail at the Bridge to Nowhere.
The bounding spots for the roads-only point is the same road up to South Mt Hawkins, the Cabin Flat Campground and Shoemaker Canyon Road.