Location data

The nearest places, all at once

"Which trailhead is closest to each member?" is distances between every member and every trailhead, then the smallest in each row. NumPy's broadcasting does the whole grid in one call: give the members' coordinates shape (members, 1) and the trailheads' shape (trailheads,), and haversine_miles returns a (members, trailheads) table.

import numpy as np

EARTH_MILES = 3959
def haversine_miles(lat1, lon1, lat2, lon2):
    lat1, lon1, lat2, lon2 = map(np.radians, (lat1, lon1, lat2, lon2))
    a = np.sin((lat2 - lat1) / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin((lon2 - lon1) / 2) ** 2
    return 2 * EARTH_MILES * np.arcsin(np.sqrt(a))

trailheads = {"Fern Hollow": (45.47, -122.20), "North Ridge": (45.70, -121.88), "Pine Gap": (44.98, -122.05)}
names = list(trailheads)
t_lat = np.array([trailheads[n][0] for n in names])
t_lon = np.array([trailheads[n][1] for n in names])

members = {"Ana": (45.52, -122.68), "Bo": (45.60, -121.95), "Cy": (44.94, -123.03)}
m_lat = np.array([p[0] for p in members.values()])[:, None]      # shape (3, 1)
m_lon = np.array([p[1] for p in members.values()])[:, None]

dist = haversine_miles(m_lat, m_lon, t_lat, t_lon)                # shape (3, 3)
print(np.round(dist, 1))
for member, row in zip(members, dist):
    print(member, "->", names[int(np.argmin(row))], round(float(row.min()), 1), "miles")

[:, None] adds a second axis, turning a row of three numbers into a column, so the subtraction inside the formula pairs every member with every trailhead. For thousands of points this is fast; for millions, real systems use a spatial index, a structure that skips far-away points without measuring them.

import numpy as np

EARTH_MILES = 3959
def haversine_miles(lat1, lon1, lat2, lon2):
    lat1, lon1, lat2, lon2 = map(np.radians, (lat1, lon1, lat2, lon2))
    a = np.sin((lat2 - lat1) / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin((lon2 - lon1) / 2) ** 2
    return 2 * EARTH_MILES * np.arcsin(np.sqrt(a))


def nearest(places, lat, lon, n=2):
    return list(places.items())[:n]


stores = {"Portland": (45.52, -122.68), "Bend": (44.06, -121.31), "Seattle": (47.61, -122.33), "Salem": (44.94, -123.03)}
print(nearest(stores, 45.0, -122.9))

GeoJSON, and the order trap

Maps and location APIs exchange shapes as GeoJSON: ordinary JSON with a few agreed keys. A Point is one position, a Polygon is a closed ring of positions (the last repeats the first), and a FeatureCollection is a list of shapes, each with its own properties such as a name.

The trap: GeoJSON writes every position as [longitude, latitude], x then y, the reverse of how people say coordinates out loud. Swap them and a store in New York (40.71, -74.01) lands at latitude -74, deep in the Southern Ocean, while Portland gets an impossible latitude of -122.68.

import json

text = """{"type": "FeatureCollection", "features": [
  {"type": "Feature", "properties": {"name": "Fern Hollow trailhead"},
   "geometry": {"type": "Point", "coordinates": [-122.20, 45.47]}},
  {"type": "Feature", "properties": {"name": "Gorge delivery zone"},
   "geometry": {"type": "Polygon", "coordinates": [[[-122.4, 45.4], [-121.8, 45.4], [-121.8, 45.8], [-122.4, 45.8], [-122.4, 45.4]]]}}
]}"""
data = json.loads(text)
for f in data["features"]:
    g = f["geometry"]
    if g["type"] == "Point":
        lon, lat = g["coordinates"]                 # longitude first
        print(f["properties"]["name"], "lat", lat, "lon", lon)
    else:
        ring = g["coordinates"][0]                  # the outer ring of the polygon
        print(f["properties"]["name"], len(ring) - 1, "corners, closed:", ring[0] == ring[-1])

A Polygon's coordinates is a list of rings: the first is the outside edge, and any others are holes cut out of it. Most zones have just one ring, which is why the code takes [0].

Is the point inside the zone?

Deciding whether a delivery address is inside a zone is the point-in-polygon test. The classic method casts a ray: draw a line from the point straight east, and count how many of the polygon's edges it crosses. An odd count means inside (you crossed in and never out), even means outside. For small zones you can treat longitude and latitude as plain x and y, because the shape distortion is the same for the point and the edges.

For each edge from corner (x1, y1) to (x2, y2), the ray at height y crosses it when y is between the two corners' heights (one strictly above, one at or below), and the crossing's x, found along the edge, is to the right of the point.

def inside(lon, lat, ring):
    xs = [p[0] for p in ring]
    ys = [p[1] for p in ring]
    return min(xs) <= lon <= max(xs) and min(ys) <= lat <= max(ys)


square = [[-122.4, 45.4], [-121.8, 45.4], [-121.8, 45.8], [-122.4, 45.8], [-122.4, 45.4]]
print(inside(-122.0, 45.6, square), inside(-121.5, 45.6, square))

import numpy as np
import matplotlib.pyplot as plt

zone = [[-122.4, 45.4], [-121.8, 45.4], [-121.8, 45.8], [-122.4, 45.8], [-122.4, 45.4]]
rng = np.random.default_rng(3)
lons = rng.uniform(-122.7, -121.5, 40)
lats = rng.uniform(45.2, 46.0, 40)

def inside(lon, lat, ring):
    crossings = 0
    for (x1, y1), (x2, y2) in zip(ring, ring[1:]):
        if (y1 > lat) != (y2 > lat) and x1 + (lat - y1) * (x2 - x1) / (y2 - y1) > lon:
            crossings += 1
    return crossings % 2 == 1

in_zone = np.array([inside(x, y, zone) for x, y in zip(lons, lats)])
fig, ax = plt.subplots(figsize=(6, 4.5))
ax.plot(*zip(*zone), color="black", linewidth=1)
ax.scatter(lons[in_zone], lats[in_zone], label=f"in zone ({in_zone.sum()})")
ax.scatter(lons[~in_zone], lats[~in_zone], marker="x", label=f"outside ({(~in_zone).sum()})")
ax.set_aspect(1 / np.cos(np.radians(lats.mean())))
ax.set_xlabel("longitude")
ax.set_ylabel("latitude")
ax.legend(loc="upper right")
plt.show()