The Story

While flying a plane over the forests of Wisconsin a couple years ago, I thought to myself "If I lost an engine, I would be awfully long way from an emergency landing airport. I wonder where the worst place to lose an engine in Wisconsin would be, if we define worst as "furthest from any airport"?

To scratch my itch, I mapped out Wisconsin's borders, downloaded a data file of airport coordinates, and wrote a program that brute-forces the answer. It took a while to run, naturally, and I've always wondered how I could have done better.

The Goal

I'm mostly interested in improvements to (or replacements for) the list comprehension that determines the furthest point from an airport. This is simply for my own information/education.

The Changes

I've made the following changes to make the code more compact for posting:

  • To simplify the problem of determining borders I've changed the state from Wisconsin to Colorado
  • I've pruned the number of airports in Colorado down to 4
  • I'm only considering coordinates located at 1 degree intervals
  • This same code is also posted on GitHub.

The Code

from geopy.distance import great_circle
from shapely.geometry import Polygon, Point
import numpy, itertools

def buildGrid():
    '''Return a grid of all latitude and longitude coordinates in the state,
       at whatever resolution. One of these will be our solution point.

    # Use shapely to build a polygon representing the land area of the state
    # For states like Colorado this is trivial but for non-square states we
    # would import the coordinate list from a separate file
    stateBorders = [
        (-109.0489, 40.9996), 
        (-109.0448, 37.0004)]
    borders = Polygon( stateBorders )

    # This represents the bounding latitude and longitude of the state. With
    # a state like Colorado this approximates the actual borders, but with 
    # a state such as Florida this would also contain much ocean.

    # Step by a full degree for a coarse grid or less (such as a tenth of 
    # a degree) for a finer solution
    latstep = 1.0
    lonstep = -1.0
    latitudeBox = numpy.arange(35,43,latstep)
    longitudeBox = numpy.arange(-101,-110,lonstep)

    # Return all points that meet the criteria for being within the 
    # actual state borders.  For example, remove all ocean from 
    # Florida, leaving only land area.
    return (i for i in itertools.product(latitudeBox,longitudeBox) if Point(i[1],i[0]).within(borders) )

def distanceToNearestAirport(point, allAirports):
    '''Given a point, determine the distance to the nearest airport
    Point is a tuple of the following format: (latitude, longitude)
    allAirports is an OrderedDict containing at least value['lat'] and value['long']
    great_circle is provided by geopy, nautical is to make the answer easy
    to verify using standard aeronautical charts'''
    return min( [great_circle(point, (eachAirport['lat'],eachAirport['long'])).nautical for eachAirport in allAirports] )

def main():

    '''Figure out the place in a state that is furthest from the nearest airport'''

    # Airport data would be stored on disk as a pickled OrderedDict
    # in the working directory. This is a small subset.
    allAirports = [
        {'long': -104.8493156, 'lat': 39.57013424, 'code': 'APA'},
        {'long': -104.6731767, 'lat': 39.86167312, 'code': 'DEN'},
        {'long': -104.5376322, 'lat': 39.78420646, 'code': 'FTG'},
        {'long': -105.1172046, 'lat': 39.90881199, 'code': 'BJC'}

    # this list comprehension walks through the grid of points and creates a
    # list of tuples: (the distance from a provided point to the nearest
    # airport, and the point itself)
    # the max function takes the highest-distance value from the list
    # of tuples, which is our answer

    howfar, furthestpoint = (max(
        ((distanceToNearestAirport(eachpoint, allAirports), eachpoint )
        for eachpoint in buildGrid() ) ) ) 

    # Since i also want to know which airport is nearest to the point, we run
    # our known furthest point back through the distance algorithm again. 
    # It outputs a tuple of the following format: (distance, airportname)
    result = min( [
        (great_circle(furthestpoint, (eachAirport['lat'], eachAirport['long']) ).nautical, eachAirport['code'])
        for eachAirport in allAirports] )

    # Output!
    print("Latitude of most desolate point: {} \n"
            "Longitude of furthest point: {} \n"
            "Nearest airport to point: {}\n"
            "Distance from point to airport: {}\n"
            .format( furthestpoint[0] , furthestpoint[1] , result[1] , result[0]) )

  • \$\begingroup\$ What do you use as your data source? \$\endgroup\$ – 200_success Oct 6 '18 at 14:41
  • \$\begingroup\$ @200_success - Data files from the FAA and (i think) the US Postal service that I've manually edited to make them suitable for import. You can see them on the "python" branch of my GitHub. \$\endgroup\$ – Steve V. Oct 6 '18 at 16:40

If this were a problem on the Euclidean plane, and not on the surface of the Earth, then it could be solved exactly by computing the Voronoi diagram for the set of airport locations, and intersecting the Voronoi regions with the bounding polygon for the state boundary. The most distant point will be one of the vertices of the intersected polygons.

A Voronoi diagram is a division of the plane into convex regions, each region containing one of the seed points (in the case in the post, one of the airports) and all other points that are closer to that seed point than any of the others.

For example, let's start with the set of Colorado airports with at least 10,000 enplanements per year:

# (longitude, latitude, code) for selected airports in Colorado.
    (-106.868889, 39.223056, 'ASE'),
    (-104.700833, 38.805833, 'COS'),
    (-104.673056, 39.861667, 'DEN'),
    (-107.753889, 37.151389, 'DRO'),
    (-106.917778, 39.642500, 'EGE'),
    (-108.526667, 39.122500, 'GJT'),
    (-106.933056, 38.533889, 'GUC'),
    (-107.217778, 40.481111, 'HDN'),
    (-107.894242, 38.509794, 'MTJ'),

# (longitude, latitude) for boundary of state of Colorado.
    (-109.0448, 37.0004),
    (-102.0424, 36.9949),
    (-102.0534, 41.0006),
    (-109.0489, 40.9996),
    (-109.0448, 37.0004),

Boundary of Colorado as a green dashed polygon. Selected airports as blue circles.

We can compute the Voronoi diagram for this set of points using scipy.spatial.Voronoi, and plot it using Matplotlib and scipy.spatial.voronoi_plot_2d, like this:

from geopy.distance import great_circle
import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial import Voronoi, voronoi_plot_2d

_, ax = plt.subplots()
points = np.array([(lon, lat) for lon, lat, code in AIRPORTS])
voronoi_plot_2d(Voronoi(points), ax=ax, show_vertices=False)
plt.plot(*points.T, 'b.')

boundary = np.array(BOUNDARY)
x, y = boundary.T
plt.plot(x, y, 'g--')
plt.xlim(round(x.min() - 1), round(x.max() + 1))
plt.ylim(round(y.min() - 1), round(y.max() + 1))

Boundary of Colorado as a green dashed polygon. Selected airports as blue circles. Finite Voronoi regions for the airports as solid black lines. Infinite Voronoi ridges as black dashed lines.

The solid black polygons are the finite Voronoi regions, and the dashed black lines are the boundaries of the Voronoi regions extending to infinity. Note that some of the finite regions (and all of the infinite regions) extend outside the state boundary, so we need to intersect these regions with the state boundary. This can be done using the Shapely geometry library and the function voronoi_polygons from this answer at Stack Overflow.

diameter = np.linalg.norm(boundary.ptp(axis=0))
boundary_polygon = Polygon(boundary)
for poly in voronoi_polygons(Voronoi(points), diameter):
    coords = np.array(poly.intersection(boundary_polygon).exterior.coords)
    plt.plot(*coords.T, 'r-')

Voronoi regions for each airport intersected with boundary of Colorado and drawn as red polygons. Airports as blue circles.

Now the most distant point is one of the vertices of the polygons resulting from intersecting the Voronoi regions with the boundary, and the Voronoi property means we only have to check each vertex against the seed point for its region:

farthest_point = None
farthest_airport = None
farthest_dist = 0
polygons = voronoi_polygons(Voronoi(points), diameter)
for airport, poly in zip(AIRPORTS, polygons):
    lon, lat, _ = airport
    coords = np.array(poly.intersection(boundary_polygon).exterior.coords)
    plt.plot(*coords.T, 'r-')
    for coord in coords:
        dist = great_circle(coord[::-1], (lat, lon))
        if dist > farthest_dist:
            farthest_point = coord
            farthest_airport = airport
            farthest_dist = dist
plt.plot(*farthest_point, 'go')
lon, lat, code = farthest_airport
plt.annotate(xy=(lon, lat), xytext=(lon + 0.1, lat + 0.05), s=code)

Voronoi regions for each airport intersected with boundary of Colorado and drawn as red polygons. Airports as blue circles. Furthest point from any airport as green circle in southeast corner of the state. Its closest airport, Colorado Springs, labelled.

This solves the planar version of the problem exactly in time \$O(\max(m, n \log n))\$, where \$m\$ is the number of vertices in the bounding polygon, and \$n\$ is the number of airports. For small areas on the surface of the Earth this will be a better approximation than the one in the post (which is only as accurate as allowed by latstep and lonstep).

  • \$\begingroup\$ The first paragraph is unnecessarily restrictive: it is equally as true in spherical geometry as in Euclidean geometry, since the clip imposed by the state boundaries will prevent the same region from containing antipodal points. The only other caveat is that the library support is less likely to be available. \$\endgroup\$ – Peter Taylor Oct 11 '18 at 11:21
  • \$\begingroup\$ It's the lack of library support that I am thinking of. There's scipy.spatial.SphericalVoronoi, but I'm not aware of a spherical geometry equivalent of Shapely's intersection method. (You can project to the plane, intersect there, and project back, but dealing with the poles and 180° longitude is tricky.) \$\endgroup\$ – Gareth Rees Oct 11 '18 at 14:44
  • \$\begingroup\$ I apologize for how long it took to accept this answer, I thought I had marked it as accepted already. \$\endgroup\$ – Steve V. Dec 15 '18 at 13:52
    # Step by a full degree for a coarse grid or less (such as a tenth of 
    # a degree) for a finer solution
    latstep = 1.0
    lonstep = -1.0

If you want to do it by sampling, there's no need to stick to a single resolution. Given a "square" (I realise that in spherical geometry it's not actually a square) and the distances of the four corners, you can bound the maximum obtainable in that region. Filter out those regions whose upper bound is definitely not good enough to win, subdivide the regions which remain, and repeat.

However, there's a far more geometrical way to do it. Compute the Voronoi diagram. The point will be either a vertex of the Voronoi diagram, an intersection between the Voronoi diagram and the boundary of the state, or a vertex of the boundary of the state.

A quick and dirty way of doing it would just be to calculate the lines of bisection between each pair of airports, calculate the points of intersections of pairs of those lines and/or the boundary of the state, and then measure the distance for each of these points. For a handful of airports this would finish running before you finished coding and debugging the Voronoi calculation.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.